diff --git a/modules/util/gfs_reader.py b/modules/util/gfs_reader.py
index 1e4a872fb089d347033ef38652a8390af0c770fd..0bc0785ce8f42c10373d19e16debda62027fde8c 100644
--- a/modules/util/gfs_reader.py
+++ b/modules/util/gfs_reader.py
@@ -1,13 +1,10 @@
import datetime
from datetime import timezone
import glob
+import os
import numpy as np
import xarray as xr
from util.util import value_to_index
-from metpy import *
-# import cartopy
-# from cartopy import *
-# import cartopy.crs as ccrs
# gfs_directory = '/apollo/cloud/Ancil_Data/clavrx_ancil_data/dynamic/gfs/'
gfs_directory = '/Users/tomrink/data/gfs/'
@@ -16,8 +13,15 @@ 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])
+NZ = plevs.shape[0]
class MyGenericException(Exception):
@@ -25,6 +29,13 @@ class MyGenericException(Exception):
self.message = message
+def get_timestamp(filename):
+ toks = filename.split('.')
+ tstr = toks[1].split('_')[0]
+ dto = datetime.datetime.strptime(tstr, gfs_date_format + '%H').replace(tzinfo=timezone.utc)
+ return dto.timestamp()
+
+
def get_time_tuple_utc(timestamp):
dt_obj = datetime.datetime.fromtimestamp(timestamp, timezone.utc)
return dt_obj, dt_obj.timetuple()
@@ -33,15 +44,18 @@ def get_time_tuple_utc(timestamp):
def get_bounding_gfs_files(timestamp):
dt_obj, time_tup = get_time_tuple_utc(timestamp)
date_str = dt_obj.strftime(gfs_date_format)
- dt_obj0 = datetime.datetime.strptime(date_str, gfs_date_format).replace(tzinfo=timezone.utc)
- dt_obj1 = dt_obj0 + datetime.timedelta(days=1)
- date_str_1 = dt_obj1.strftime(gfs_date_format)
-
- flist_0 = glob.glob(gfs_directory+'gfs.'+date_str+'??_F012.h5')
- flist_1 = glob.glob(gfs_directory+'gfs.'+date_str_1+'00_F012.h5')
- if (len(flist_0) == 0) or (len(flist_1) == 0):
+ dt_obj = datetime.datetime.strptime(date_str, gfs_date_format).replace(tzinfo=timezone.utc)
+ dt_obj_r = dt_obj + datetime.timedelta(days=1)
+ date_str_r = dt_obj_r.strftime(gfs_date_format)
+ dt_obj_l = dt_obj - datetime.timedelta(days=1)
+ date_str_l = dt_obj_l.strftime(gfs_date_format)
+
+ flist_l = glob.glob(gfs_directory+'gfs.'+date_str_l+'??_F012.h5')
+ flist = glob.glob(gfs_directory+'gfs.'+date_str+'??_F012.h5')
+ flist_r = glob.glob(gfs_directory+'gfs.'+date_str_r+'??_F012.h5')
+ filelist = flist_l + flist + flist_r
+ if len(filelist) == 0:
return None, None, None, None
- filelist = flist_0 + flist_1
ftimes = []
for pname in filelist: # TODO: make better with regular expressions (someday)
@@ -77,14 +91,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]
@@ -111,13 +124,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]
@@ -147,26 +159,28 @@ def get_horz_layer_s(xr_dataset, fld_names, press, lon_range=None, lat_range=Non
return sub_fld
-def get_time_interpolated_layer(ds_0, ds_1, time_0, time_1, time, fld_name, press, lon_range=None, lat_range=None):
- lyr_0 = get_horz_layer(ds_0, fld_name, press, lon_range=lon_range, lat_range=lat_range)
- lyr_1 = get_horz_layer(ds_1, fld_name, press, lon_range=lon_range, lat_range=lat_range)
+def get_time_interpolated_layer(xr_dataset_s, time_s, time, fld_name, press, lon_range=None, lat_range=None, method='linear'):
+ layer_s = []
+ for ds in xr_dataset_s:
+ lyr = get_horz_layer(ds, fld_name, press, lon_range=lon_range, lat_range=lat_range)
+ layer_s.append(lyr)
- lyr = xr.concat([lyr_0, lyr_1], 'time')
- lyr = lyr.assign_coords(time=[time_0, time_1])
+ lyr = xr.concat(layer_s, 'time')
+ lyr = lyr.assign_coords(time=time_s)
- intrp_lyr = lyr.interp(time=time, method='linear')
- intrp_lyr = intrp_lyr.expand_dims('channel')
- intrp_lyr = intrp_lyr.assign_coords(channel=[fld_name])
+ intrp_lyr = lyr.interp(time=time, method=method)
return intrp_lyr
-def get_time_interpolated_layer_s(ds_0, ds_1, time_0, time_1, time, fld_name_s, press, lon_range=None, lat_range=None):
- lyr_s_0 = get_horz_layer_s(ds_0, fld_name_s, press, lon_range=lon_range, lat_range=lat_range)
- lyr_s_1 = get_horz_layer_s(ds_1, fld_name_s, press, lon_range=lon_range, lat_range=lat_range)
+def get_time_interpolated_layer_s(xr_dataset_s, time_s, time, fld_name_s, press, lon_range=None, lat_range=None, method='linear'):
+ layer_s = []
+ for ds in xr_dataset_s:
+ lyr = get_horz_layer_s(ds, fld_name_s, press, lon_range=lon_range, lat_range=lat_range)
+ layer_s.append(lyr)
- lyr = xr.concat([lyr_s_0, lyr_s_1], 'time')
- lyr = lyr.assign_coords(time=[time_0, time_1])
+ lyr = xr.concat(layer_s, 'time')
+ lyr = lyr.assign_coords(time=time_s)
intrp_lyr = lyr.interp(time=time, method='linear')
@@ -177,7 +191,6 @@ def get_vert_profile(xr_dataset, fld_name, lons, lats):
if lon360: # convert -180/+180 to 0,360
lons = np.where(lons < 0, lons + 360, lons)
- plevs = xr_dataset['pressure levels']
fld = xr_dataset[fld_name]
fld = fld.assign_coords(fakeDim2=lon_coords, fakeDim1=lat_coords, fakeDim0=plevs)
@@ -194,7 +207,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:
@@ -230,31 +242,31 @@ def get_point(xr_dataset, fld_name, lons, lats):
return intrp_fld
-def get_time_interpolated_vert_profile(ds_0, ds_1, time0, time1, fld_name, time, lons, lats):
-
- prof_0 = get_vert_profile(ds_0, fld_name, lons, lats)
- prof_1 = get_vert_profile(ds_1, fld_name, lons, lats)
+def get_time_interpolated_vert_profile(xr_dataset_s, time_s, fld_name, time, lons, lats, method='linear'):
+ prof_s = []
+ for ds in xr_dataset_s:
+ vp = get_vert_profile(ds, fld_name, lons, lats)
+ prof_s.append(vp)
- prof = xr.concat([prof_0, prof_1], 'time')
- prof = prof.assign_coords(time=[time0, time1])
-
- intrp_prof = prof.interp(time=time, method='linear')
+ prof = xr.concat(prof_s, 'time')
+ prof = prof.assign_coords(time=time_s)
+ intrp_prof = prof.interp(time=time, method=method)
intrp_prof = intrp_prof.values
return intrp_prof
-def get_time_interpolated_vert_profile_s(ds_0, ds_1, time0, time1, fld_name_s, time, lons, lats):
-
- prof_0 = get_vert_profile_s(ds_0, fld_name_s, lons, lats)
- prof_1 = get_vert_profile_s(ds_1, fld_name_s, lons, lats)
-
- prof = xr.concat([prof_0, prof_1], 'time')
- prof = prof.assign_coords(time=[time0, time1])
+def get_time_interpolated_vert_profile_s(xr_dataset_s, time_s, fld_name_s, time, lons, lats, method='linear'):
+ prof_s = []
+ for ds in xr_dataset_s:
+ vp = get_vert_profile_s(ds, fld_name_s, lons, lats)
+ prof_s.append(vp)
- intrp_prof = prof.interp(time=time, method='linear')
+ prof = xr.concat(prof_s, 'time')
+ prof = prof.assign_coords(time=time_s)
+ intrp_prof = prof.interp(time=time, method=method)
intrp_prof = intrp_prof.values
return intrp_prof
@@ -272,3 +284,103 @@ def get_time_interpolated_point(ds_0, ds_1, time0, time1, fld_name, time, lons,
intrp_vals = intrp_vals.values
return intrp_vals
+
+
+def get_time_interpolated_voxel(xr_dataset_s, time_s, time, fld_name, lon, lat, press, x_width=5, y_width=5, z_width=3, method='linear'):
+ vox_s = []
+ for ds in xr_dataset_s:
+ vox = get_voxel(ds, fld_name, lon, lat, press, x_width=x_width, y_width=y_width, z_width=z_width)
+ vox_s.append(vox)
+
+ vox = xr.concat(vox_s, 'time')
+ vox = vox.assign_coords(time=time_s)
+
+ intrp_vox = vox.interp(time=time, method=method)
+
+ return intrp_vox
+
+
+def get_voxel(xr_dataset, fld_name, lon, lat, press, x_width=5, y_width=5, z_width=3):
+ if lon360:
+ if lon < 0:
+ lon += 360
+
+ fld = xr_dataset[fld_name]
+
+ 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)
+ x_h = int(x_width / 2)
+ p_h = int(z_width / 2)
+
+ y_start = y_c - y_h
+ x_start = x_c - x_h
+ z_start = p_c - p_h
+ if y_start < 0 or x_start < 0 or z_start < 0:
+ return None
+
+ y_stop = y_c + y_h + 1
+ x_stop = x_c + x_h + 1
+ z_stop = p_c + p_h + 1
+ if y_stop > NY-1 or x_stop > NX-1 or z_stop > NZ-1:
+ return None
+
+ sub_fld = fld[y_start:y_stop, x_start:x_stop, z_start:z_stop]
+
+ sub_fld = sub_fld.expand_dims('channel')
+ sub_fld = sub_fld.assign_coords(channel=[fld_name], fakeDim2=lon_coords[x_start:x_stop],
+ fakeDim1=lat_coords[y_start:y_stop], fakeDim0=plevs[z_start:z_stop])
+
+ return sub_fld
+
+
+def get_time_interpolated_voxel_s(xr_dataset_s, time_s, time, fld_name_s, lon, lat, press, x_width=5, y_width=5, z_width=3, method='linear'):
+ vox_s = []
+ for ds in xr_dataset_s:
+ vox = get_voxel_s(ds, fld_name_s, lon, lat, press, x_width=x_width, y_width=y_width, z_width=z_width)
+ vox_s.append(vox)
+
+ vox = xr.concat(vox_s, 'time')
+ vox = vox.assign_coords(time=time_s)
+
+ intrp_vox = vox.interp(time=time, method=method)
+
+ return intrp_vox
+
+
+def get_voxel_s(xr_dataset, fld_name_s, lon, lat, press, x_width=5, y_width=5, z_width=3):
+ if lon360:
+ if lon < 0:
+ lon += 360
+
+ 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)
+ x_h = int(x_width / 2)
+ p_h = int(z_width / 2)
+
+ y_start = y_c - y_h
+ x_start = x_c - x_h
+ z_start = p_c - p_h
+ if y_start < 0 or x_start < 0 or z_start < 0:
+ return None
+
+ y_stop = y_c + y_h + 1
+ x_stop = x_c + x_h + 1
+ z_stop = p_c + p_h + 1
+ if y_stop > NY-1 or x_stop > NX-1 or z_stop > NZ-1:
+ return None
+
+ sub_fld_s = []
+ for name in fld_name_s:
+ fld = xr_dataset[name]
+ sub_fld = fld[y_start:y_stop, x_start:x_stop, z_start:z_stop]
+ sub_fld_s.append(sub_fld)
+ sub_fld = xr.concat(sub_fld_s, 'channel')
+
+ sub_fld = sub_fld.assign_coords(channel=fld_name_s, fakeDim2=lon_coords[x_start:x_stop],
+ fakeDim1=lat_coords[y_start:y_stop], fakeDim0=plevs[z_start:z_stop])
+
+ return sub_fld
\ No newline at end of file