diff --git a/modules/aeolus/aeolus_amv.py b/modules/aeolus/aeolus_amv.py index 55e3e8325577a4fc3a89853d56447390f1085ecf..c307d86d0864c402e616c9eddd82f467e3a5c134 100644 --- a/modules/aeolus/aeolus_amv.py +++ b/modules/aeolus/aeolus_amv.py @@ -5,6 +5,7 @@ import numpy as np import xarray as xr from netCDF4 import Dataset, Dimension, Variable from aeolus.geos_nav import GEOSNavigation +from aeolus.datasource import get_datasource from util.util import haversine_np, LatLonTuple, GenericException from amv.intercompare import best_fit, bin_data_by, get_press_bin_ranges, spd_dir_from_uv, uv_from_spd_dir, \ direction_difference, run_best_fit_gfs @@ -22,213 +23,213 @@ num_elems = 5424 num_lines = 5424 -class AMVFiles: - - def __init__(self, files_path, file_time_span, pattern, band='14'): - self.flist = glob.glob(files_path + pattern) - if len(self.flist) == 0: - raise GenericException('no matching files found in: ' + files_path) - self.band = band - self.ftimes = [] - self.span_seconds = datetime.timedelta(minutes=file_time_span).seconds - for pname in self.flist: - dto = self.get_datetime(pname) - dto_start = dto - dto_end = dto + datetime.timedelta(minutes=file_time_span) - self.ftimes.append((dto_start.timestamp(), dto_end.timestamp())) - self.ftimes = np.array(self.ftimes) - self.flist = np.array(self.flist) - sidxs = np.argsort(self.ftimes[:, 0]) - self.ftimes = self.ftimes[sidxs, :] - self.flist = self.flist[sidxs] - - def get_datetime(self, pathname): - pass - - def get_navigation(self): - pass - - def get_file_containing_time(self, timestamp): - k = -1 - for i in range(len(self.ftimes.shape[0])): - if (timestamp >= self.ftimes[i, 0]) and (timestamp < self.ftimes[i, 1]): - k = i - break - if k < 0: - return None, None, None - - return self.flist[k], self.ftimes[k, 0], k - - def get_file(self, timestamp): - diff = self.ftimes[:, 0] - timestamp - midx = np.argmin(np.abs(diff)) - if np.abs(self.ftimes[midx, 0] - timestamp) < self.span_seconds: - return self.flist[midx], self.ftimes[midx, 0], midx - else: - return None, None, None - - def get_parameters(self): - pass - - def get_out_parameters(self): - pass - - def get_meta_dict(self): - pass - - -class Framework(AMVFiles): - def __init__(self, files_path, file_time_span, band='14'): - super().__init__(files_path, file_time_span, '*WINDS_AMV_EN-'+band+'*.nc', band) - - self.elem_name = 'Element' - self.line_name = 'Line' - self.lon_name = 'Longitude' - self.lat_name = 'Latitude' - - self.out_params = ['Lon', 'Lat', 'Element', 'Line', 'pressure', 'wind_speed', 'wind_direction', 'BestFitPresLvl'] - self.params = ['MedianPress', 'Wind_Speed', 'Wind_Dir', 'BestFitPresLvl'] - self.meta_dict = {'Lon': ('degrees east', 'f4'), 'Lat': ('degrees north', 'f4'), 'Element': (None, 'i4'), 'Line': (None, 'i4'), - 'pressure': ('hPa', 'f4'), 'wind_speed': ('m s-1', 'f4'), 'wind_direction': ('degrees', 'f4'), - 'BestFitPresLvl': ('hPa', 'f4')} - - def get_parameters(self): - return self.params - - def get_out_parameters(self): - return self.out_params - - def get_meta_dict(self): - return self.meta_dict - - def get_navigation(self): - return GEOSNavigation(sub_lon=-75.0) - - def get_datetime(self, pathname): - fname = os.path.split(pathname)[1] - toks = fname.split('_') - dstr = toks[4] - tstr = toks[5] - dtstr = dstr + tstr - dto = datetime.datetime.strptime(dtstr, '%Y%j%H%M').replace(tzinfo=timezone.utc) - - return dto - - -class FrameworkCloudHeight(AMVFiles): - def __init__(self, files_path, file_time_span): - super().__init__(files_path, file_time_span, '*_CLOUD_HEIGHT_EN'+'*.nc', None) - - self.elem_name = 'Element' - self.line_name = 'Line' - self.lon_name = 'Longitude' - self.lat_name = 'Latitude' - - self.out_params = ['CldTopPres', 'CldTopHght', 'CldOptDpth'] - self.params = ['CldTopPres', 'CldTopHght', 'CldOptDpth'] - self.meta_dict = {'CldTopPres': ('hPa', 'f4'), 'CldTopHght': ('km', 'f4'), 'CldOptDpth': ('km', 'f4')} - - def get_parameters(self): - return self.params - - def get_out_parameters(self): - return self.out_params - - def get_meta_dict(self): - return self.meta_dict - - def get_navigation(self): - return GEOSNavigation(sub_lon=-75.0) - - def get_datetime(self, pathname): - fname = os.path.split(pathname)[1] - toks = fname.split('_') - dstr = toks[4] - tstr = toks[5] - dtstr = dstr + tstr - dto = datetime.datetime.strptime(dtstr, '%Y%j%H%M').replace(tzinfo=timezone.utc) - - return dto - - -class OPS(AMVFiles): - def __init__(self, files_path, file_time_span, band='14'): - super().__init__(files_path, file_time_span, 'OR_ABI-L2-DMWF*'+'C'+band+'*.nc', band) - - self.elem_name = None - self.line_name = None - self.lon_name = 'lon' - self.lat_name = 'lat' - - self.out_params = ['Lon', 'Lat', 'Element', 'Line', 'pressure', 'wind_speed', 'wind_direction'] - self.params = ['pressure', 'wind_speed', 'wind_direction'] - self.meta_dict = {'Lon': ('degrees east', 'f4'), 'Lat': ('degrees north', 'f4'), 'Element': (None, 'i4'), 'Line': (None, 'i4'), - 'pressure': ('hPa', 'f4'), 'wind_speed': ('m s-1', 'f4'), 'wind_direction': ('degrees', 'f4')} - - def get_navigation(self): - # return GEOSNavigation(sub_lon=-75.2) ? - return GEOSNavigation(sub_lon=-75.0) - - def get_datetime(self, pathname): - fname = os.path.split(pathname)[1] - toks = fname.split('_') - dtstr = toks[3] - dtstr = dtstr[:-3] - dto = datetime.datetime.strptime(dtstr, 's%Y%j%H%M').replace(tzinfo=timezone.utc) - - return dto - - def get_parameters(self): - return self.params - - def get_out_parameters(self): - return self.out_params - - def get_meta_dict(self): - return self.meta_dict - - -class CarrStereo(AMVFiles): - def __init__(self, files_path, file_time_span, band='14'): - super().__init__(files_path, file_time_span, 'tdw_qc_GOES*'+'ch_'+band+'.nc', band) - - self.elem_name = 'Element' - self.line_name = 'Line' - self.lon_name = 'Lon' - self.lat_name = 'Lat' - - self.out_params = ['Lon', 'Lat', 'Element', 'Line', 'V_3D_u', 'V_3D_v', 'H_3D', 'pres', 'Fcst_Spd', 'Fcst_Dir', 'SatZen', - 'InversionFlag', 'CloudPhase', 'CloudType'] - - self.params = ['V_3D', 'H_3D', 'pres', 'Fcst_Spd', 'Fcst_Dir', 'SatZen', - 'InversionFlag', 'CloudPhase', 'CloudType'] - - self.meta_dict = {'H_3D': ('m', 'f4'), 'pres': ('hPa', 'f4'), 'Fcst_Spd': ('m s-1', 'f4'), - 'Fcst_Dir': ('degree', 'f4'), - 'SatZen': ('degree', 'f4'), 'InversionFlag': (None, 'u1'), - 'CloudPhase': (None, 'u1'), 'CloudType': (None, 'u1'), - 'V_3D_u': ('m s-1', 'f4'), 'V_3D_v': ('m s-1', 'f4'), 'Lon': ('degrees east', 'f4'), - 'Lat': ('degrees north', 'f4'), 'Element': (None, 'i4'), 'Line': (None, 'i4')} - - def get_navigation(self): - return GEOSNavigation(sub_lon=-137.0) - - def get_datetime(self, pathname): - fname = os.path.split(pathname)[1] - toks = fname.split('_') - dtstr = toks[3] - dto = datetime.datetime.strptime(dtstr, '%Y%j.%H%M.ch').replace(tzinfo=timezone.utc) - - return dto - - def get_parameters(self): - return self.params - - def get_out_parameters(self): - return self.out_params - - def get_meta_dict(self): - return self.meta_dict +# class AMVFiles: +# +# def __init__(self, files_path, file_time_span, pattern, band='14'): +# self.flist = glob.glob(files_path + pattern) +# if len(self.flist) == 0: +# raise GenericException('no matching files found in: ' + files_path) +# self.band = band +# self.ftimes = [] +# self.span_seconds = datetime.timedelta(minutes=file_time_span).seconds +# for pname in self.flist: +# dto = self.get_datetime(pname) +# dto_start = dto +# dto_end = dto + datetime.timedelta(minutes=file_time_span) +# self.ftimes.append((dto_start.timestamp(), dto_end.timestamp())) +# self.ftimes = np.array(self.ftimes) +# self.flist = np.array(self.flist) +# sidxs = np.argsort(self.ftimes[:, 0]) +# self.ftimes = self.ftimes[sidxs, :] +# self.flist = self.flist[sidxs] +# +# def get_datetime(self, pathname): +# pass +# +# def get_navigation(self): +# pass +# +# def get_file_containing_time(self, timestamp): +# k = -1 +# for i in range(len(self.ftimes.shape[0])): +# if (timestamp >= self.ftimes[i, 0]) and (timestamp < self.ftimes[i, 1]): +# k = i +# break +# if k < 0: +# return None, None, None +# +# return self.flist[k], self.ftimes[k, 0], k +# +# def get_file(self, timestamp): +# diff = self.ftimes[:, 0] - timestamp +# midx = np.argmin(np.abs(diff)) +# if np.abs(self.ftimes[midx, 0] - timestamp) < self.span_seconds: +# return self.flist[midx], self.ftimes[midx, 0], midx +# else: +# return None, None, None +# +# def get_parameters(self): +# pass +# +# def get_out_parameters(self): +# pass +# +# def get_meta_dict(self): +# pass +# +# +# class Framework(AMVFiles): +# def __init__(self, files_path, file_time_span, band='14'): +# super().__init__(files_path, file_time_span, '*WINDS_AMV_EN-'+band+'*.nc', band) +# +# self.elem_name = 'Element' +# self.line_name = 'Line' +# self.lon_name = 'Longitude' +# self.lat_name = 'Latitude' +# +# self.out_params = ['Lon', 'Lat', 'Element', 'Line', 'pressure', 'wind_speed', 'wind_direction', 'BestFitPresLvl'] +# self.params = ['MedianPress', 'Wind_Speed', 'Wind_Dir', 'BestFitPresLvl'] +# self.meta_dict = {'Lon': ('degrees east', 'f4'), 'Lat': ('degrees north', 'f4'), 'Element': (None, 'i4'), 'Line': (None, 'i4'), +# 'pressure': ('hPa', 'f4'), 'wind_speed': ('m s-1', 'f4'), 'wind_direction': ('degrees', 'f4'), +# 'BestFitPresLvl': ('hPa', 'f4')} +# +# def get_parameters(self): +# return self.params +# +# def get_out_parameters(self): +# return self.out_params +# +# def get_meta_dict(self): +# return self.meta_dict +# +# def get_navigation(self): +# return GEOSNavigation(sub_lon=-75.0) +# +# def get_datetime(self, pathname): +# fname = os.path.split(pathname)[1] +# toks = fname.split('_') +# dstr = toks[4] +# tstr = toks[5] +# dtstr = dstr + tstr +# dto = datetime.datetime.strptime(dtstr, '%Y%j%H%M').replace(tzinfo=timezone.utc) +# +# return dto +# +# +# class FrameworkCloudHeight(AMVFiles): +# def __init__(self, files_path, file_time_span): +# super().__init__(files_path, file_time_span, '*_CLOUD_HEIGHT_EN'+'*.nc', None) +# +# self.elem_name = 'Element' +# self.line_name = 'Line' +# self.lon_name = 'Longitude' +# self.lat_name = 'Latitude' +# +# self.out_params = ['CldTopPres', 'CldTopHght', 'CldOptDpth'] +# self.params = ['CldTopPres', 'CldTopHght', 'CldOptDpth'] +# self.meta_dict = {'CldTopPres': ('hPa', 'f4'), 'CldTopHght': ('km', 'f4'), 'CldOptDpth': ('km', 'f4')} +# +# def get_parameters(self): +# return self.params +# +# def get_out_parameters(self): +# return self.out_params +# +# def get_meta_dict(self): +# return self.meta_dict +# +# def get_navigation(self): +# return GEOSNavigation(sub_lon=-75.0) +# +# def get_datetime(self, pathname): +# fname = os.path.split(pathname)[1] +# toks = fname.split('_') +# dstr = toks[4] +# tstr = toks[5] +# dtstr = dstr + tstr +# dto = datetime.datetime.strptime(dtstr, '%Y%j%H%M').replace(tzinfo=timezone.utc) +# +# return dto +# +# +# class OPS(AMVFiles): +# def __init__(self, files_path, file_time_span, band='14'): +# super().__init__(files_path, file_time_span, 'OR_ABI-L2-DMWF*'+'C'+band+'*.nc', band) +# +# self.elem_name = None +# self.line_name = None +# self.lon_name = 'lon' +# self.lat_name = 'lat' +# +# self.out_params = ['Lon', 'Lat', 'Element', 'Line', 'pressure', 'wind_speed', 'wind_direction'] +# self.params = ['pressure', 'wind_speed', 'wind_direction'] +# self.meta_dict = {'Lon': ('degrees east', 'f4'), 'Lat': ('degrees north', 'f4'), 'Element': (None, 'i4'), 'Line': (None, 'i4'), +# 'pressure': ('hPa', 'f4'), 'wind_speed': ('m s-1', 'f4'), 'wind_direction': ('degrees', 'f4')} +# +# def get_navigation(self): +# # return GEOSNavigation(sub_lon=-75.2) ? +# return GEOSNavigation(sub_lon=-75.0) +# +# def get_datetime(self, pathname): +# fname = os.path.split(pathname)[1] +# toks = fname.split('_') +# dtstr = toks[3] +# dtstr = dtstr[:-3] +# dto = datetime.datetime.strptime(dtstr, 's%Y%j%H%M').replace(tzinfo=timezone.utc) +# +# return dto +# +# def get_parameters(self): +# return self.params +# +# def get_out_parameters(self): +# return self.out_params +# +# def get_meta_dict(self): +# return self.meta_dict +# +# +# class CarrStereo(AMVFiles): +# def __init__(self, files_path, file_time_span, band='14'): +# super().__init__(files_path, file_time_span, 'tdw_qc_GOES*'+'ch_'+band+'.nc', band) +# +# self.elem_name = 'Element' +# self.line_name = 'Line' +# self.lon_name = 'Lon' +# self.lat_name = 'Lat' +# +# self.out_params = ['Lon', 'Lat', 'Element', 'Line', 'V_3D_u', 'V_3D_v', 'H_3D', 'pres', 'Fcst_Spd', 'Fcst_Dir', 'SatZen', +# 'InversionFlag', 'CloudPhase', 'CloudType'] +# +# self.params = ['V_3D', 'H_3D', 'pres', 'Fcst_Spd', 'Fcst_Dir', 'SatZen', +# 'InversionFlag', 'CloudPhase', 'CloudType'] +# +# self.meta_dict = {'H_3D': ('m', 'f4'), 'pres': ('hPa', 'f4'), 'Fcst_Spd': ('m s-1', 'f4'), +# 'Fcst_Dir': ('degree', 'f4'), +# 'SatZen': ('degree', 'f4'), 'InversionFlag': (None, 'u1'), +# 'CloudPhase': (None, 'u1'), 'CloudType': (None, 'u1'), +# 'V_3D_u': ('m s-1', 'f4'), 'V_3D_v': ('m s-1', 'f4'), 'Lon': ('degrees east', 'f4'), +# 'Lat': ('degrees north', 'f4'), 'Element': (None, 'i4'), 'Line': (None, 'i4')} +# +# def get_navigation(self): +# return GEOSNavigation(sub_lon=-137.0) +# +# def get_datetime(self, pathname): +# fname = os.path.split(pathname)[1] +# toks = fname.split('_') +# dtstr = toks[3] +# dto = datetime.datetime.strptime(dtstr, '%Y%j.%H%M.ch').replace(tzinfo=timezone.utc) +# +# return dto +# +# def get_parameters(self): +# return self.params +# +# def get_out_parameters(self): +# return self.out_params +# +# def get_meta_dict(self): +# return self.meta_dict def get_amvs(amv_files, timestamp, filepath=None): @@ -425,8 +426,8 @@ def create_file2(filename, raob_to_amv_dct, raob_dct, amv_files): def run_best_fit_all(): - amv_files = Framework('/apollo/cloud/scratch/AMV_BUST/FMWK_TEST3/', 10) - prd_files = FrameworkCloudHeight('/apollo/cloud/scratch/AMV_BUST/FMWK_TEST3/', 10) + amv_files = get_datasource('/apollo/cloud/scratch/AMV_BUST/FMWK_TEST3/', 10, 'FMWK') + prd_files = get_datasource('/apollo/cloud/scratch/AMV_BUST/FMWK_TEST3/', 10, 'FMWK_CLD_HGT') raob_dir = '/home/rink/data/raob/' raob_files = ['raob_soundings20191117_0000.cdf', @@ -464,17 +465,17 @@ def run_best_fit_all(): for k, file in enumerate(raob_files): raob_dct, ts = get_raob_dict_cdf(raob_dir+file) m_d = match_amvs_to_raobs(raob_dct, ts, amv_files=amv_files) - #bf_dct = run_best_fit(m_d, raob_dct, gfs_dir+gfs_files[k]) - #prd_dct = get_product_at_locs(m_d, ts, prd_files) - #out_list.append((bf_dct, prd_dct)) - - amvs = get_amvs(amv_files, ts) - bfs = run_best_fit_gfs(amvs, gfs_dir+gfs_files[k], amv_lat_idx=0, amv_lon_idx=1, amv_prs_idx=4, amv_spd_idx=5, amv_dir_idx=6) - alons = amvs[:, 0] - alats = amvs[:, 1] - prds = get_product_at_lat_lons(prd_files, ts, alons, alats, filepath=None) - - out_list.append((amvs, np.array(bfs), prds)) + bf_dct = run_best_fit(m_d, raob_dct, gfs_dir+gfs_files[k]) + prd_dct = get_product_at_locs(m_d, ts, prd_files) + out_list.append((bf_dct, prd_dct)) + + # amvs = get_amvs(amv_files, ts) + # bfs = run_best_fit_gfs(amvs, gfs_dir+gfs_files[k], amv_lat_idx=0, amv_lon_idx=1, amv_prs_idx=4, amv_spd_idx=5, amv_dir_idx=6) + # alons = amvs[:, 0] + # alats = amvs[:, 1] + # prds = get_product_at_lat_lons(prd_files, ts, alons, alats, filepath=None) + # + # out_list.append((amvs, np.array(bfs), prds)) return out_list @@ -628,9 +629,9 @@ def run_best_fit(raob_to_amv_dct, raob_dct, gfs_filename=None): bf = best_fit(amv_spd, amv_dir, amv_prs, amv_lat, amv_lon, raob_spd, raob_dir, raob_prs) bf_list.append(bf) - pdiff = amv_prs - raob_prs - lev_idx = np.argmin(np.abs(pdiff)) - if np.abs(pdiff[lev_idx]) > 100.0: + pdiff = np.abs(amv_prs - raob_prs) + lev_idx = np.argmin(pdiff) + if pdiff[lev_idx] > 60.0: tup = (raob_spd[lev_idx], raob_dir[lev_idx], raob_prs[lev_idx], -9) else: tup = (raob_spd[lev_idx], raob_dir[lev_idx], raob_prs[lev_idx], 0) @@ -1721,7 +1722,7 @@ def create_amv_to_aeolus_match_file(aeolus_files_dir, amv_files_dir, outfile=Non amv_files = None if amv_source == 'CARR': - amv_files = CarrStereo(amv_files_dir, 60, band) + amv_files = get_datasource(amv_files_dir, 60, 'CARR', band=band) m_d = match_amvs_to_aeolus(a_d, amv_files_dir, amv_source, band, amv_files)