snapshot...

modules/util/gfs_reader.py

@@ -12,8 +12,14 @@ gfs_date_format = '%y%m%d'
 # force incoming longitude to (0, 360) to match GFS
 lon360 = True
 # GFS half degree resolution
-lat_coords = np.linspace(-90, 90, 361)
-lon_coords = np.linspace(0, 359.5, 720)
+NX = 720
+NY = 361
+lat_coords = np.linspace(-90, 90, NY)
+lon_coords = np.linspace(0, 359.5, NX)
+
+plevs = np.array([10.0, 20.0, 30.0, 50.0, 70.0, 100.0, 150.0, 200.0, 250.0, 300.0,
+                  350.0, 400.0, 450.0, 500.0, 550.0, 600.0, 650.0, 700.0, 750.0, 800.0,
+                  850.0, 900.0, 925.0, 950.0, 975.0, 1000.0])
 
 
 class MyGenericException(Exception):
@@ -80,14 +86,13 @@ def get_bounding_gfs_files(timestamp):
 
 
 def get_horz_layer(xr_dataset, fld_name, press, lon_range=None, lat_range=None):
-    plevs = xr_dataset['pressure levels']
-    p_idx = value_to_index(plevs.values, press)
+    p_idx = value_to_index(plevs, press)
     fld = xr_dataset[fld_name]
 
     x_lo = 0
-    x_hi = xr_dataset.fakeDim2.shape[0]
+    x_hi = NX
     y_lo = 0
-    y_hi = xr_dataset.fakeDim1.shape[0]
+    y_hi = NY
 
     if lon_range is not None:
         lon_lo = lon_range[0]
@@ -114,13 +119,12 @@ def get_horz_layer(xr_dataset, fld_name, press, lon_range=None, lat_range=None):
 
 
 def get_horz_layer_s(xr_dataset, fld_names, press, lon_range=None, lat_range=None):
-    plevs = xr_dataset['pressure levels']
-    p_idx = value_to_index(plevs.values, press)
+    p_idx = value_to_index(plevs, press)
 
     x_lo = 0
-    x_hi = xr_dataset.fakeDim2.shape[0]
+    x_hi = NX
     y_lo = 0
-    y_hi = xr_dataset.fakeDim1.shape[0]
+    y_hi = NY
 
     if lon_range is not None:
         lon_lo = lon_range[0]
@@ -199,7 +203,6 @@ def get_vert_profile_s(xr_dataset, fld_name_s, lons, lats):
     if lon360:  # convert -180,+180 to 0,360
         lons = np.where(lons < 0, lons + 360, lons)
 
-    plevs = xr_dataset['pressure levels']
 
     fld_s = []
     for fld_name in fld_name_s:
@@ -297,10 +300,10 @@ def get_voxel(xr_dataset, fld_name, lon, lat, press, x_width=5, y_width=5, z_wid
     if lon360:
         if lon < 0:
             lon += 360
-    plevs = xr_dataset['pressure levels']
+
     fld = xr_dataset[fld_name]
 
-    p_c = value_to_index(plevs.values, press)
+    p_c = value_to_index(plevs, press)
     x_c = value_to_index(lon_coords, lon)
     y_c = value_to_index(lat_coords, lat)
     y_h = int(y_width / 2)
@@ -334,9 +337,8 @@ def get_voxel_s(xr_dataset, fld_name_s, lon, lat, press, x_width=5, y_width=5, z
     if lon360:
         if lon < 0:
             lon += 360
-    plevs = xr_dataset['pressure levels']
 
-    p_c = value_to_index(plevs.values, press)
+    p_c = value_to_index(plevs, press)
     x_c = value_to_index(lon_coords, lon)
     y_c = value_to_index(lat_coords, lat)
     y_h = int(y_width / 2)