From 9c85dbea189825db5365152108d3231062b60c0a Mon Sep 17 00:00:00 2001 From: tomrink <rink@ssec.wisc.edu> Date: Thu, 19 Nov 2020 21:00:35 -0600 Subject: [PATCH] snapshot... --- modules/util/gfs_reader.py | 32 +++++++++++++++++--------------- 1 file changed, 17 insertions(+), 15 deletions(-) diff --git a/modules/util/gfs_reader.py b/modules/util/gfs_reader.py index 48c43795..5c360a06 100644 --- a/modules/util/gfs_reader.py +++ b/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) -- GitLab