snapshot...

modules/contrail/util.py

@@ -13,6 +13,9 @@ from util.util import get_grid_values_all
 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.units import units
 
 gfs_files = GFSfiles('/Users/tomrink/data/contrail/gfs/')
 # GEOSNavigation needs to be updated to support GOES-18
@@ -93,7 +96,19 @@ def extract(mask_image, image_ts, clavrx_path):
     lon_range = [np.min(contrail_lons), np.max(contrail_lons)]
     lat_range = [np.min(contrail_lats), np.max(contrail_lats)]
 
-    wind_3d = get_volume(xr_dataset, ['u-wind','v-wind'], lon_range=[lon_range[0], lon_range[1]], lat_range=[lat_range[0], lat_range[1]])
+    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]])
+    uwind_3d = uwind_3d.transpose('Pressure', 'Latitude', 'Longitude')
+    vwind_3d = vwind_3d.transpose('Pressure', 'Latitude', 'Longitude')
+    temp_3d = temp_3d.transpose('Pressure', 'Latitude', 'Longitude')
+
+    # uwind_2d = get_horz_layer(xr_dataset, 'u-wind', 500.0, lon_range=[lon_range[0], lon_range[1]], lat_range=[lat_range[0], lat_range[1]])
+    # 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)
+    static_3d = static_stability(temp_3d.coords['Pressure'] * units.hPa, temp_3d)
 
     xr_dataset.close()