Skip to content
Snippets Groups Projects
Commit f321cc27 authored by Tom Rink's avatar Tom Rink
Browse files

Merge remote-tracking branch 'origin/master'

parents 4d6e2f73 1e474eb0
No related branches found
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
modules/util/gfs_reader.py
+ 165
− 53
View file @ f321cc27
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
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment