snapshot...

modules/contrail/util.py

@@ -14,7 +14,7 @@ from util.gfs_reader import *
 from util.geos_nav import GEOSNavigation
 from aeolus.datasource import GFSfiles
 import cartopy.crs as ccrs
-from metpy.calc import shearing_deformation, static_stability
+from metpy.calc import shearing_deformation, static_stability, wind_speed, first_derivative
 from metpy.units import units
 
 gfs_files = GFSfiles('/Users/tomrink/data/contrail/gfs/')
@@ -68,9 +68,8 @@ def extract(mask_image, image_ts, clavrx_path):
     contrail_lons = contrail_lons[keep]
     contrail_lats = contrail_lats[keep]
 
-    bins = np.arange(100, 1000, 100)
-
     # Indexes of contrail_press for individual bins
+    bins = np.arange(100, 1000, 100)
     binned_indexes = np.digitize(contrail_press, bins)
 
     # Store the indexes in a dictionary where the key is the bin number and value is the list of indexes
@@ -78,19 +77,19 @@ def extract(mask_image, image_ts, clavrx_path):
     for bin_num in np.unique(binned_indexes):
         bins_dict[bin_num] = np.where(binned_indexes == bin_num)[0]
 
-    # MetPy for shearing deformation, static stability, vertical shear?, others...
-    # This does point by point computation of model parameters for each contrail pixel
-    voxel_dict = {key: [] for key in bins_dict.keys()}
-    for key in bins_dict.keys():
-        print('working on pressure level: ', bins[key])
-        for c_idx in bins_dict[key]:
-            lon = contrail_lons[c_idx]
-            lat = contrail_lats[c_idx]
-            press = contrail_press[c_idx]
-
-            wind_3d = get_voxel_s(xr_dataset, ['u-wind','v-wind'], lon, lat, press)
-            if wind_3d is not None:
-                voxel_dict[key].append(wind_3d)
+    # # This does point by point computation of model parameters for each contrail pixel
+    # voxel_dict = {key: [] for key in bins_dict.keys()}
+    # for key in bins_dict.keys():
+    #     print('working on pressure level: ', bins[key])
+    #     for c_idx in bins_dict[key]:
+    #         lon = contrail_lons[c_idx]
+    #         lat = contrail_lats[c_idx]
+    #         press = contrail_press[c_idx]
+    #
+    #         wind_3d = get_voxel_s(xr_dataset, ['u-wind','v-wind'], lon, lat, press)
+    #         if wind_3d is not None:
+    #             voxel_dict[key].append(wind_3d)
+    # ------------------------------------------------------------------------------------
 
     # This section will compute model parameters in bulk for the region then pull for each contrail pixel
     lon_range = [np.min(contrail_lons), np.max(contrail_lons)]
@@ -99,6 +98,7 @@ def extract(mask_image, image_ts, clavrx_path):
     uwind_3d = get_volume(xr_dataset, 'u-wind', 'm s-1', lon_range=[lon_range[0], lon_range[1]], lat_range=[lat_range[0], lat_range[1]])
     vwind_3d = get_volume(xr_dataset, 'v-wind', 'm s-1', lon_range=[lon_range[0], lon_range[1]], lat_range=[lat_range[0], lat_range[1]])
     temp_3d = get_volume(xr_dataset, 'temperature', 'degK', lon_range=[lon_range[0], lon_range[1]], lat_range=[lat_range[0], lat_range[1]])
+    rh_3d = get_volume(xr_dataset, 'rh', '%', lon_range=[lon_range[0], lon_range[1]], lat_range=[lat_range[0], lat_range[1]])
 
     uwind_3d = uwind_3d.transpose('Pressure', 'Latitude', 'Longitude')
     vwind_3d = vwind_3d.transpose('Pressure', 'Latitude', 'Longitude')
@@ -108,8 +108,10 @@ def extract(mask_image, image_ts, clavrx_path):
     # vwind_2d = get_horz_layer(xr_dataset, 'v-wind', 500.0, lon_range=[lon_range[0], lon_range[1]], lat_range=[lat_range[0], lat_range[1]])
     # shear_2d = shearing_deformation(uwind_2d, vwind_2d)
 
-    shear_3d = shearing_deformation(uwind_3d, vwind_3d)
+    horz_shear_3d = shearing_deformation(uwind_3d, vwind_3d)
     static_3d = static_stability(temp_3d.coords['Pressure'] * units.hPa, temp_3d)
+    horz_wind_spd_3d = wind_speed(uwind_3d, vwind_3d)
+    vert_shear_3d = first_derivative(horz_wind_spd_3d, axis=0, x=temp_3d.coords['Pressure'] * units.hPa)
 
     xr_dataset.close()