From 3cb5e6504d30fb5f19d398a421cfe03b43239f10 Mon Sep 17 00:00:00 2001
From: tomrink <rink@ssec.wisc.edu>
Date: Tue, 31 Oct 2023 21:14:09 -0500
Subject: [PATCH] snapshot...
---
modules/util/util.py | 698 -------------------------------------------
1 file changed, 698 deletions(-)
diff --git a/modules/util/util.py b/modules/util/util.py
index c4b2b13d..d8db73d5 100644
--- a/modules/util/util.py
+++ b/modules/util/util.py
@@ -15,7 +15,6 @@ from scipy.interpolate import RectBivariateSpline, interp2d
from scipy.ndimage import gaussian_filter
from scipy.signal import medfilt2d
from cartopy.crs import Geostationary, Globe
-import geos_nav
LatLonTuple = namedtuple('LatLonTuple', ['lat', 'lon'])
@@ -667,19 +666,6 @@ def add_noise(data, noise_scale=0.01, seed=None, copy=True):
return data
-# Keep for reference. These pkl files are incompatible with latest version of Cartopy
-# f = open(ancillary_path+'geos_crs_goes16_FD.pkl', 'rb')
-# geos_goes16_fd = pickle.load(f)
-# f.close()
-#
-# f = open(ancillary_path+'geos_crs_goes16_CONUS.pkl', 'rb')
-# geos_goes16_conus = pickle.load(f)
-# f.close()
-#
-# f = open(ancillary_path+'geos_crs_H08_FD.pkl', 'rb')
-# geos_h08_fd = pickle.load(f)
-# f.close()
-
crs_goes16_fd = Geostationary(central_longitude=-75.0, satellite_height=35786023.0, sweep_axis='x',
globe=Globe(ellipse=None, semimajor_axis=6378137.0, semiminor_axis=6356752.31414,
inverse_flattening=298.2572221))
@@ -692,9 +678,6 @@ crs_h08_fd = Geostationary(central_longitude=140.7, satellite_height=35785.863,
globe=Globe(ellipse=None, semimajor_axis=6378.137, semiminor_axis=6356.7523,
inverse_flattening=298.25702))
-geos_goes16_fd = geos_nav.GEOSNavigation()
-geos_goes16_conus = geos_nav.GEOSNavigation(CFAC=5.6E-05, LFAC=-5.6E-05, COFF=-0.101332, LOFF=0.128212, num_elems=2500, num_lines=1500)
-
def get_cartopy_crs(satellite, domain):
if satellite == 'GOES16':
@@ -816,285 +799,6 @@ exmp_file_conus = '/Users/tomrink/data/OR_ABI-L1b-RadC-M6C14_G16_s20193140811215
# Full Disk
exmp_file_fd = '/Users/tomrink/data/OR_ABI-L1b-RadF-M6C16_G16_s20212521800223_e20212521809542_c20212521809596.nc'
-# keep for reference
-# if domain == 'CONUS':
-# exmpl_ds = xr.open_dataset(exmp_file_conus)
-# elif domain == 'FD':
-# exmpl_ds = xr.open_dataset(exmp_file_fd)
-# mdat = exmpl_ds.metpy.parse_cf('Rad')
-# geos = mdat.metpy.cartopy_crs
-# xlen = mdat.x.values.size
-# ylen = mdat.y.values.size
-# exmpl_ds.close()
-
-# Taiwan domain:
-# lon, lat = 120.955098, 23.834310
-# elem, line = (1789, 1505)
-# # UR from Taiwan
-# lon, lat = 135.0, 35.0
-# elem_ur, line_ur = (2499, 995)
-taiwan_i0 = 1079
-taiwan_j0 = 995
-taiwan_lenx = 1420
-taiwan_leny = 1020
-# geos.transform_point(135.0, 35.0, ccrs.PlateCarree(), False)
-# geos.transform_point(106.61, 13.97, ccrs.PlateCarree(), False)
-taiwain_extent = [-3342, -502, 1470, 3510] # GEOS coordinates, not line, elem
-
-
-# ------------ This code will not be needed when we implement a Fully Convolutional CNN -----------------------------------
-# Generate and return tiles of name_list parameters
-def make_for_full_domain_predict(h5f, name_list=None, satellite='GOES16', domain='FD', res_fac=1,
- data_extent=None):
- w_x = 16
- w_y = 16
- i_0 = 0
- j_0 = 0
- s_x = int(w_x / res_fac)
- s_y = int(w_y / res_fac)
-
- geos, xlen, xmin, xmax, ylen, ymin, ymax = get_cartopy_crs(satellite, domain)
- if satellite == 'GOES16':
- if data_extent is not None:
- ll_pt_a = data_extent[0]
- ll_pt_b = data_extent[1]
-
- if satellite == 'H08':
- xlen = taiwan_lenx
- ylen = taiwan_leny
- i_0 = taiwan_i0
- j_0 = taiwan_j0
- elif satellite == 'H09':
- xlen = taiwan_lenx
- ylen = taiwan_leny
- i_0 = taiwan_i0
- j_0 = taiwan_j0
-
- grd_dct = {name: None for name in name_list}
-
- cnt_a = 0
- for ds_name in name_list:
- fill_value, fill_value_name = get_fill_attrs(ds_name)
- gvals = get_grid_values(h5f, ds_name, j_0, i_0, None, num_j=ylen, num_i=xlen, fill_value_name=fill_value_name, fill_value=fill_value)
- if gvals is not None:
- grd_dct[ds_name] = gvals
- cnt_a += 1
-
- if cnt_a > 0 and cnt_a != len(name_list):
- raise GenericException('weirdness')
-
- grd_dct_n = {name: [] for name in name_list}
-
- n_x = int(xlen/s_x) - 1
- n_y = int(ylen/s_y) - 1
-
- r_x = xlen - (n_x * s_x)
- x_d = 0 if r_x >= w_x else int((w_x - r_x)/s_x)
- n_x -= x_d
-
- r_y = ylen - (n_y * s_y)
- y_d = 0 if r_y >= w_y else int((w_y - r_y)/s_y)
- n_y -= y_d
-
- ll = [j_0 + j*s_y for j in range(n_y)]
- cc = [i_0 + i*s_x for i in range(n_x)]
-
- for ds_name in name_list:
- for j in range(n_y):
- j_ul = j * s_y
- j_ul_b = j_ul + w_y
- for i in range(n_x):
- i_ul = i * s_x
- i_ul_b = i_ul + w_x
- grd_dct_n[ds_name].append(grd_dct[ds_name][j_ul:j_ul_b, i_ul:i_ul_b])
-
- grd_dct = {name: None for name in name_list}
- for ds_name in name_list:
- grd_dct[ds_name] = np.stack(grd_dct_n[ds_name])
-
- return grd_dct, ll, cc
-
-
-def make_for_full_domain_predict_viirs_clavrx(h5f, name_list=None, res_fac=1, day_night='DAY', use_dnb=False):
- w_x = 16
- w_y = 16
- i_0 = 0
- j_0 = 0
- s_x = int(w_x / res_fac)
- s_y = int(w_y / res_fac)
-
- ylen = h5f['scan_lines_along_track_direction'].shape[0]
- xlen = h5f['pixel_elements_along_scan_direction'].shape[0]
-
- use_nl_comp = False
- if (day_night == 'NIGHT' or day_night == 'AUTO') and use_dnb:
- use_nl_comp = True
-
- grd_dct = {name: None for name in name_list}
-
- cnt_a = 0
- for ds_name in name_list:
- name = ds_name
-
- if use_nl_comp:
- if ds_name == 'cld_reff_dcomp':
- name = 'cld_reff_nlcomp'
- elif ds_name == 'cld_opd_dcomp':
- name = 'cld_opd_nlcomp'
-
- fill_value, fill_value_name = get_fill_attrs(name)
- gvals = get_grid_values(h5f, name, j_0, i_0, None, num_j=ylen, num_i=xlen, fill_value_name=fill_value_name, fill_value=fill_value)
- if gvals is not None:
- grd_dct[ds_name] = gvals
- cnt_a += 1
-
- if cnt_a > 0 and cnt_a != len(name_list):
- raise GenericException('weirdness')
-
- # TODO: need to investigate discrepencies with compute_lwc_iwc
- # if use_nl_comp:
- # cld_phase = get_grid_values(h5f, 'cloud_phase', j_0, i_0, None, num_j=ylen, num_i=xlen)
- # cld_dz = get_grid_values(h5f, 'cld_geo_thick', j_0, i_0, None, num_j=ylen, num_i=xlen)
- # reff = grd_dct['cld_reff_dcomp']
- # opd = grd_dct['cld_opd_dcomp']
- #
- # lwc_nlcomp, iwc_nlcomp = compute_lwc_iwc(cld_phase, reff, opd, cld_dz)
- # grd_dct['iwc_dcomp'] = iwc_nlcomp
- # grd_dct['lwc_dcomp'] = lwc_nlcomp
-
- grd_dct_n = {name: [] for name in name_list}
-
- n_x = int(xlen/s_x) - 1
- n_y = int(ylen/s_y) - 1
-
- r_x = xlen - (n_x * s_x)
- x_d = 0 if r_x >= w_x else int((w_x - r_x)/s_x)
- n_x -= x_d
-
- r_y = ylen - (n_y * s_y)
- y_d = 0 if r_y >= w_y else int((w_y - r_y)/s_y)
- n_y -= y_d
-
- ll = [j_0 + j*s_y for j in range(n_y)]
- cc = [i_0 + i*s_x for i in range(n_x)]
-
- for ds_name in name_list:
- for j in range(n_y):
- j_ul = j * s_y
- j_ul_b = j_ul + w_y
- for i in range(n_x):
- i_ul = i * s_x
- i_ul_b = i_ul + w_x
- grd_dct_n[ds_name].append(grd_dct[ds_name][j_ul:j_ul_b, i_ul:i_ul_b])
-
- grd_dct = {name: None for name in name_list}
- for ds_name in name_list:
- grd_dct[ds_name] = np.stack(grd_dct_n[ds_name])
-
- lats = get_grid_values(h5f, 'latitude', j_0, i_0, None, num_j=ylen, num_i=xlen)
- lons = get_grid_values(h5f, 'longitude', j_0, i_0, None, num_j=ylen, num_i=xlen)
-
- ll_2d, cc_2d = np.meshgrid(ll, cc, indexing='ij')
-
- lats = lats[ll_2d, cc_2d]
- lons = lons[ll_2d, cc_2d]
-
- return grd_dct, ll, cc, lats, lons
-
-
-def make_for_full_domain_predict2(h5f, satellite='GOES16', domain='FD', res_fac=1):
- w_x = 16
- w_y = 16
- i_0 = 0
- j_0 = 0
- s_x = int(w_x / res_fac)
- s_y = int(w_y / res_fac)
-
- geos, xlen, xmin, xmax, ylen, ymin, ymax = get_cartopy_crs(satellite, domain)
- if satellite == 'H08':
- xlen = taiwan_lenx
- ylen = taiwan_leny
- i_0 = taiwan_i0
- j_0 = taiwan_j0
-
- n_x = int(xlen/s_x)
- n_y = int(ylen/s_y)
-
- solzen = get_grid_values(h5f, 'solar_zenith_angle', j_0, i_0, None, num_j=ylen, num_i=xlen)
- satzen = get_grid_values(h5f, 'sensor_zenith_angle', j_0, i_0, None, num_j=ylen, num_i=xlen)
- solzen = solzen[0:(n_y-1)*s_y:s_y, 0:(n_x-1)*s_x:s_x]
- satzen = satzen[0:(n_y-1)*s_y:s_y, 0:(n_x-1)*s_x:s_x]
-
- return solzen, satzen
-# -------------------------------------------------------------------------------------------
-
-
-def prepare_evaluate(h5f, name_list, satellite='GOES16', domain='FD', res_fac=1, offset=0):
- w_x = 16
- w_y = 16
- i_0 = 0
- j_0 = 0
- s_x = w_x // res_fac
- s_y = w_y // res_fac
-
- geos, xlen, xmin, xmax, ylen, ymin, ymax = get_cartopy_crs(satellite, domain)
- if satellite == 'H08':
- xlen = taiwan_lenx
- ylen = taiwan_leny
- i_0 = taiwan_i0
- j_0 = taiwan_j0
- elif satellite == 'H09':
- xlen = taiwan_lenx
- ylen = taiwan_leny
- i_0 = taiwan_i0
- j_0 = taiwan_j0
-
- n_x = (xlen // s_x)
- n_y = (ylen // s_y)
-
- # r_x = xlen - (n_x * s_x)
- # x_d = 0 if r_x >= w_x else int((w_x - r_x)/s_x)
- # n_x -= x_d
- #
- # r_y = ylen - (n_y * s_y)
- # y_d = 0 if r_y >= w_y else int((w_y - r_y)/s_y)
- # n_y -= y_d
-
- ll = [(offset+j_0) + j*s_y for j in range(n_y)]
- cc = [(offset+i_0) + i*s_x for i in range(n_x)]
-
- grd_dct_n = {name: [] for name in name_list}
-
- cnt_a = 0
- for ds_name in name_list:
- fill_value, fill_value_name = get_fill_attrs(ds_name)
- gvals = get_grid_values(h5f, ds_name, j_0, i_0, None, num_j=ylen, num_i=xlen, fill_value_name=fill_value_name, fill_value=fill_value)
- gvals = np.expand_dims(gvals, axis=0)
- if gvals is not None:
- grd_dct_n[ds_name] = gvals
- cnt_a += 1
-
- if cnt_a > 0 and cnt_a != len(name_list):
- raise GenericException('weirdness')
-
- solzen = get_grid_values(h5f, 'solar_zenith_angle', j_0, i_0, None, num_j=ylen, num_i=xlen)
- satzen = get_grid_values(h5f, 'sensor_zenith_angle', j_0, i_0, None, num_j=ylen, num_i=xlen)
- solzen = solzen[0:n_y*s_y:s_y, 0:n_x*s_x:s_x]
- satzen = satzen[0:n_y*s_y:s_y, 0:n_x*s_x:s_x]
-
- return grd_dct_n, solzen, satzen, ll, cc
-
-
-flt_level_ranges_str = {k: None for k in range(5)}
-flt_level_ranges_str[0] = '0_2000'
-flt_level_ranges_str[1] = '2000_4000'
-flt_level_ranges_str[2] = '4000_6000'
-flt_level_ranges_str[3] = '6000_8000'
-flt_level_ranges_str[4] = '8000_15000'
-
-# flt_level_ranges_str = {k: None for k in range(1)}
-# flt_level_ranges_str[0] = 'column'
-
def get_cf_nav_parameters(satellite='GOES16', domain='FD'):
param_dct = None
@@ -1152,371 +856,6 @@ def get_cf_nav_parameters(satellite='GOES16', domain='FD'):
return param_dct
-def write_icing_file(clvrx_str_time, output_dir, preds_dct, probs_dct, x, y, lons, lats, elems, lines):
- outfile_name = output_dir + 'icing_prediction_'+clvrx_str_time+'.h5'
- h5f_out = h5py.File(outfile_name, 'w')
-
- dim_0_name = 'x'
- dim_1_name = 'y'
-
- prob_s = []
- pred_s = []
-
- flt_lvls = list(preds_dct.keys())
- for flvl in flt_lvls:
- preds = preds_dct[flvl]
- pred_s.append(preds)
- icing_pred_ds = h5f_out.create_dataset('icing_prediction_level_'+flt_level_ranges_str[flvl], data=preds, dtype='i2')
- icing_pred_ds.attrs.create('coordinates', data='y x')
- icing_pred_ds.attrs.create('grid_mapping', data='Projection')
- icing_pred_ds.attrs.create('missing', data=-1)
- icing_pred_ds.dims[0].label = dim_0_name
- icing_pred_ds.dims[1].label = dim_1_name
-
- for flvl in flt_lvls:
- probs = probs_dct[flvl]
- prob_s.append(probs)
- icing_prob_ds = h5f_out.create_dataset('icing_probability_level_'+flt_level_ranges_str[flvl], data=probs, dtype='f4')
- icing_prob_ds.attrs.create('coordinates', data='y x')
- icing_prob_ds.attrs.create('grid_mapping', data='Projection')
- icing_prob_ds.attrs.create('missing', data=-1.0)
- icing_prob_ds.dims[0].label = dim_0_name
- icing_prob_ds.dims[1].label = dim_1_name
-
- prob_s = np.stack(prob_s, axis=-1)
- max_prob = np.max(prob_s, axis=2)
-
- icing_prob_ds = h5f_out.create_dataset('max_icing_probability_column', data=max_prob, dtype='f4')
- icing_prob_ds.attrs.create('coordinates', data='y x')
- icing_prob_ds.attrs.create('grid_mapping', data='Projection')
- icing_prob_ds.attrs.create('missing', data=-1.0)
- icing_prob_ds.dims[0].label = dim_0_name
- icing_prob_ds.dims[1].label = dim_1_name
-
- max_lvl = np.argmax(prob_s, axis=2)
-
- icing_pred_ds = h5f_out.create_dataset('max_icing_probability_level', data=max_lvl, dtype='i2')
- icing_pred_ds.attrs.create('coordinates', data='y x')
- icing_pred_ds.attrs.create('grid_mapping', data='Projection')
- icing_pred_ds.attrs.create('missing', data=-1)
- icing_pred_ds.dims[0].label = dim_0_name
- icing_pred_ds.dims[1].label = dim_1_name
-
- lon_ds = h5f_out.create_dataset('longitude', data=lons, dtype='f4')
- lon_ds.attrs.create('units', data='degrees_east')
- lon_ds.attrs.create('long_name', data='icing prediction longitude')
- lon_ds.dims[0].label = dim_0_name
- lon_ds.dims[1].label = dim_1_name
-
- lat_ds = h5f_out.create_dataset('latitude', data=lats, dtype='f4')
- lat_ds.attrs.create('units', data='degrees_north')
- lat_ds.attrs.create('long_name', data='icing prediction latitude')
- lat_ds.dims[0].label = dim_0_name
- lat_ds.dims[1].label = dim_1_name
-
- proj_ds = h5f_out.create_dataset('Projection', data=0, dtype='b')
- proj_ds.attrs.create('long_name', data='Himawari Imagery Projection')
- proj_ds.attrs.create('grid_mapping_name', data='geostationary')
- proj_ds.attrs.create('sweep_angle_axis', data='y')
- proj_ds.attrs.create('units', data='rad')
- proj_ds.attrs.create('semi_major_axis', data=6378.137)
- proj_ds.attrs.create('semi_minor_axis', data=6356.7523)
- proj_ds.attrs.create('inverse_flattening', data=298.257)
- proj_ds.attrs.create('perspective_point_height', data=35785.863)
- proj_ds.attrs.create('latitude_of_projection_origin', data=0.0)
- proj_ds.attrs.create('longitude_of_projection_origin', data=140.7)
- proj_ds.attrs.create('CFAC', data=20466275)
- proj_ds.attrs.create('LFAC', data=20466275)
- proj_ds.attrs.create('COFF', data=2750.5)
- proj_ds.attrs.create('LOFF', data=2750.5)
-
- if x is not None:
- x_ds = h5f_out.create_dataset('x', data=x, dtype='f8')
- x_ds.dims[0].label = dim_0_name
- x_ds.attrs.create('units', data='rad')
- x_ds.attrs.create('standard_name', data='projection_x_coordinate')
- x_ds.attrs.create('long_name', data='GOES PUG W-E fixed grid viewing angle')
- x_ds.attrs.create('scale_factor', data=5.58879902955962e-05)
- x_ds.attrs.create('add_offset', data=-0.153719917308037)
- x_ds.attrs.create('CFAC', data=20466275)
- x_ds.attrs.create('COFF', data=2750.5)
-
- y_ds = h5f_out.create_dataset('y', data=y, dtype='f8')
- y_ds.dims[0].label = dim_1_name
- y_ds.attrs.create('units', data='rad')
- y_ds.attrs.create('standard_name', data='projection_y_coordinate')
- y_ds.attrs.create('long_name', data='GOES PUG S-N fixed grid viewing angle')
- y_ds.attrs.create('scale_factor', data=-5.58879902955962e-05)
- y_ds.attrs.create('add_offset', data=0.153719917308037)
- y_ds.attrs.create('LFAC', data=20466275)
- y_ds.attrs.create('LOFF', data=2750.5)
-
- if elems is not None:
- elem_ds = h5f_out.create_dataset('elems', data=elems, dtype='i2')
- elem_ds.dims[0].label = dim_0_name
- line_ds = h5f_out.create_dataset('lines', data=lines, dtype='i2')
- line_ds.dims[0].label = dim_1_name
- pass
-
- h5f_out.close()
-
-
-def write_icing_file_nc4(clvrx_str_time, output_dir, preds_dct, probs_dct,
- x, y, lons, lats, elems, lines, satellite='GOES16', domain='CONUS',
- has_time=False, use_nan=False, prob_thresh=0.5, bt_10_4=None):
- outfile_name = output_dir + 'icing_prediction_'+clvrx_str_time+'.nc'
- rootgrp = Dataset(outfile_name, 'w', format='NETCDF4')
-
- rootgrp.setncattr('Conventions', 'CF-1.7')
-
- dim_0_name = 'x'
- dim_1_name = 'y'
- time_dim_name = 'time'
- geo_coords = 'time y x'
-
- dim_0 = rootgrp.createDimension(dim_0_name, size=x.shape[0])
- dim_1 = rootgrp.createDimension(dim_1_name, size=y.shape[0])
- dim_time = rootgrp.createDimension(time_dim_name, size=1)
-
- tvar = rootgrp.createVariable('time', 'f8', time_dim_name)
- tvar[0] = get_timestamp(clvrx_str_time)
- tvar.units = 'seconds since 1970-01-01 00:00:00'
-
- if not has_time:
- var_dim_list = [dim_1_name, dim_0_name]
- else:
- var_dim_list = [time_dim_name, dim_1_name, dim_0_name]
-
- prob_s = []
-
- flt_lvls = list(preds_dct.keys())
- for flvl in flt_lvls:
- preds = preds_dct[flvl]
-
- icing_pred_ds = rootgrp.createVariable('icing_prediction_level_'+flt_level_ranges_str[flvl], 'i2', var_dim_list)
- icing_pred_ds.setncattr('coordinates', geo_coords)
- icing_pred_ds.setncattr('grid_mapping', 'Projection')
- icing_pred_ds.setncattr('missing', -1)
- if has_time:
- preds = preds.reshape((1, y.shape[0], x.shape[0]))
- icing_pred_ds[:,] = preds
-
- for flvl in flt_lvls:
- probs = probs_dct[flvl]
- prob_s.append(probs)
-
- icing_prob_ds = rootgrp.createVariable('icing_probability_level_'+flt_level_ranges_str[flvl], 'f4', var_dim_list)
- icing_prob_ds.setncattr('coordinates', geo_coords)
- icing_prob_ds.setncattr('grid_mapping', 'Projection')
- if not use_nan:
- icing_prob_ds.setncattr('missing', -1.0)
- else:
- icing_prob_ds.setncattr('missing', np.nan)
- if has_time:
- probs = probs.reshape((1, y.shape[0], x.shape[0]))
- if use_nan:
- probs = np.where(probs < prob_thresh, np.nan, probs)
- icing_prob_ds[:,] = probs
-
- prob_s = np.stack(prob_s, axis=-1)
- max_prob = np.max(prob_s, axis=2)
- if use_nan:
- max_prob = np.where(max_prob < prob_thresh, np.nan, max_prob)
- if has_time:
- max_prob = max_prob.reshape(1, y.shape[0], x.shape[0])
-
- icing_prob_ds = rootgrp.createVariable('max_icing_probability_column', 'f4', var_dim_list)
- icing_prob_ds.setncattr('coordinates', geo_coords)
- icing_prob_ds.setncattr('grid_mapping', 'Projection')
- if not use_nan:
- icing_prob_ds.setncattr('missing', -1.0)
- else:
- icing_prob_ds.setncattr('missing', np.nan)
- icing_prob_ds[:,] = max_prob
-
- prob_s = np.where(prob_s < prob_thresh, -1.0, prob_s)
- max_lvl = np.where(np.all(prob_s == -1, axis=2), -1, np.argmax(prob_s, axis=2))
- if use_nan:
- max_lvl = np.where(max_lvl == -1, np.nan, max_lvl)
- if has_time:
- max_lvl = max_lvl.reshape((1, y.shape[0], x.shape[0]))
-
- icing_pred_ds = rootgrp.createVariable('max_icing_probability_level', 'i2', var_dim_list)
- icing_pred_ds.setncattr('coordinates', geo_coords)
- icing_pred_ds.setncattr('grid_mapping', 'Projection')
- icing_pred_ds.setncattr('missing', -1)
- icing_pred_ds[:,] = max_lvl
-
- if bt_10_4 is not None:
- bt_ds = rootgrp.createVariable('bt_10_4', 'f4', var_dim_list)
- bt_ds.setncattr('coordinates', geo_coords)
- bt_ds.setncattr('grid_mapping', 'Projection')
- bt_ds[:,] = bt_10_4
-
- lon_ds = rootgrp.createVariable('longitude', 'f4', [dim_1_name, dim_0_name])
- lon_ds.units = 'degrees_east'
- lon_ds[:,] = lons
-
- lat_ds = rootgrp.createVariable('latitude', 'f4', [dim_1_name, dim_0_name])
- lat_ds.units = 'degrees_north'
- lat_ds[:,] = lats
-
- cf_nav_dct = get_cf_nav_parameters(satellite, domain)
-
- if satellite == 'H08':
- long_name = 'Himawari Imagery Projection'
- elif satellite == 'H09':
- long_name = 'Himawari Imagery Projection'
- elif satellite == 'GOES16':
- long_name = 'GOES-16/17 Imagery Projection'
-
- proj_ds = rootgrp.createVariable('Projection', 'b')
- proj_ds.setncattr('long_name', long_name)
- proj_ds.setncattr('grid_mapping_name', 'geostationary')
- proj_ds.setncattr('sweep_angle_axis', cf_nav_dct['sweep_angle_axis'])
- proj_ds.setncattr('semi_major_axis', cf_nav_dct['semi_major_axis'])
- proj_ds.setncattr('semi_minor_axis', cf_nav_dct['semi_minor_axis'])
- proj_ds.setncattr('inverse_flattening', cf_nav_dct['inverse_flattening'])
- proj_ds.setncattr('perspective_point_height', cf_nav_dct['perspective_point_height'])
- proj_ds.setncattr('latitude_of_projection_origin', cf_nav_dct['latitude_of_projection_origin'])
- proj_ds.setncattr('longitude_of_projection_origin', cf_nav_dct['longitude_of_projection_origin'])
-
- if x is not None:
- x_ds = rootgrp.createVariable(dim_0_name, 'f8', [dim_0_name])
- x_ds.units = 'rad'
- x_ds.setncattr('axis', 'X')
- x_ds.setncattr('standard_name', 'projection_x_coordinate')
- x_ds.setncattr('long_name', 'fixed grid viewing angle')
- x_ds.setncattr('scale_factor', cf_nav_dct['x_scale_factor'])
- x_ds.setncattr('add_offset', cf_nav_dct['x_add_offset'])
- x_ds[:] = x
-
- y_ds = rootgrp.createVariable(dim_1_name, 'f8', [dim_1_name])
- y_ds.units = 'rad'
- y_ds.setncattr('axis', 'Y')
- y_ds.setncattr('standard_name', 'projection_y_coordinate')
- y_ds.setncattr('long_name', 'fixed grid viewing angle')
- y_ds.setncattr('scale_factor', cf_nav_dct['y_scale_factor'])
- y_ds.setncattr('add_offset', cf_nav_dct['y_add_offset'])
- y_ds[:] = y
-
- if elems is not None:
- elem_ds = rootgrp.createVariable('elems', 'i2', [dim_0_name])
- elem_ds[:] = elems
- line_ds = rootgrp.createVariable('lines', 'i2', [dim_1_name])
- line_ds[:] = lines
- pass
-
- rootgrp.close()
-
-
-def write_icing_file_nc4_viirs(clvrx_str_time, output_dir, preds_dct, probs_dct, lons, lats,
- has_time=False, use_nan=False, prob_thresh=0.5, bt_10_4=None):
- outfile_name = output_dir + 'icing_prediction_'+clvrx_str_time+'.nc'
- rootgrp = Dataset(outfile_name, 'w', format='NETCDF4')
-
- rootgrp.setncattr('Conventions', 'CF-1.7')
-
- dim_0_name = 'x'
- dim_1_name = 'y'
- time_dim_name = 'time'
- geo_coords = 'longitude latitude'
-
- dim_1_len, dim_0_len = lons.shape
-
- dim_0 = rootgrp.createDimension(dim_0_name, size=dim_0_len)
- dim_1 = rootgrp.createDimension(dim_1_name, size=dim_1_len)
- dim_time = rootgrp.createDimension(time_dim_name, size=1)
-
- tvar = rootgrp.createVariable('time', 'f8', time_dim_name)
- tvar[0] = get_timestamp(clvrx_str_time)
- tvar.units = 'seconds since 1970-01-01 00:00:00'
-
- if not has_time:
- var_dim_list = [dim_1_name, dim_0_name]
- else:
- var_dim_list = [time_dim_name, dim_1_name, dim_0_name]
-
- prob_s = []
-
- flt_lvls = list(preds_dct.keys())
- for flvl in flt_lvls:
- preds = preds_dct[flvl]
-
- icing_pred_ds = rootgrp.createVariable('icing_prediction_level_'+flt_level_ranges_str[flvl], 'i2', var_dim_list)
- icing_pred_ds.setncattr('coordinates', geo_coords)
- icing_pred_ds.setncattr('grid_mapping', 'Projection')
- icing_pred_ds.setncattr('missing', -1)
- if has_time:
- preds = preds.reshape((1, dim_1_len, dim_0_len))
- icing_pred_ds[:,] = preds
-
- for flvl in flt_lvls:
- probs = probs_dct[flvl]
- prob_s.append(probs)
-
- icing_prob_ds = rootgrp.createVariable('icing_probability_level_'+flt_level_ranges_str[flvl], 'f4', var_dim_list)
- icing_prob_ds.setncattr('coordinates', geo_coords)
- icing_prob_ds.setncattr('grid_mapping', 'Projection')
- if not use_nan:
- icing_prob_ds.setncattr('missing', -1.0)
- else:
- icing_prob_ds.setncattr('missing', np.nan)
- if has_time:
- probs = probs.reshape((1, dim_1_len, dim_0_len))
- if use_nan:
- probs = np.where(probs < prob_thresh, np.nan, probs)
- icing_prob_ds[:,] = probs
-
- prob_s = np.stack(prob_s, axis=-1)
- max_prob = np.max(prob_s, axis=2)
- if use_nan:
- max_prob = np.where(max_prob < prob_thresh, np.nan, max_prob)
- if has_time:
- max_prob = max_prob.reshape(1, dim_1_len, dim_0_len)
-
- icing_prob_ds = rootgrp.createVariable('max_icing_probability_column', 'f4', var_dim_list)
- icing_prob_ds.setncattr('coordinates', geo_coords)
- icing_prob_ds.setncattr('grid_mapping', 'Projection')
- if not use_nan:
- icing_prob_ds.setncattr('missing', -1.0)
- else:
- icing_prob_ds.setncattr('missing', np.nan)
- icing_prob_ds[:,] = max_prob
-
- prob_s = np.where(prob_s < prob_thresh, -1.0, prob_s)
- max_lvl = np.where(np.all(prob_s == -1, axis=2), -1, np.argmax(prob_s, axis=2))
- if use_nan:
- max_lvl = np.where(max_lvl == -1, np.nan, max_lvl)
- if has_time:
- max_lvl = max_lvl.reshape((1, dim_1_len, dim_0_len))
-
- icing_pred_ds = rootgrp.createVariable('max_icing_probability_level', 'i2', var_dim_list)
- icing_pred_ds.setncattr('coordinates', geo_coords)
- icing_pred_ds.setncattr('grid_mapping', 'Projection')
- icing_pred_ds.setncattr('missing', -1)
- icing_pred_ds[:,] = max_lvl
-
- if bt_10_4 is not None:
- bt_ds = rootgrp.createVariable('bt_10_4', 'f4', var_dim_list)
- bt_ds.setncattr('coordinates', geo_coords)
- bt_ds.setncattr('grid_mapping', 'Projection')
- bt_ds[:,] = bt_10_4
-
- lon_ds = rootgrp.createVariable('longitude', 'f4', [dim_1_name, dim_0_name])
- lon_ds.units = 'degrees_east'
- lon_ds[:,] = lons
-
- lat_ds = rootgrp.createVariable('latitude', 'f4', [dim_1_name, dim_0_name])
- lat_ds.units = 'degrees_north'
- lat_ds[:,] = lats
-
- proj_ds = rootgrp.createVariable('Projection', 'b')
- proj_ds.setncattr('grid_mapping_name', 'latitude_longitude')
-
- rootgrp.close()
-
-
def write_cld_prods_file_nc4(clvrx_str_time, outfile_name, cloud_fraction, cloud_frac_opd,
x, y, elems, lines, satellite='GOES16', domain='CONUS',
has_time=False):
@@ -1680,40 +1019,3 @@ def median_filter_2d(z, kernel_size=3):
def median_filter_2d_single(z, kernel_size=3):
return medfilt2d(z, kernel_size=kernel_size)
-
-
-def get_training_parameters(day_night='DAY', l1b_andor_l2='both', satellite='GOES16', use_dnb=False):
- if day_night == 'DAY':
- train_params_l2 = ['cld_height_acha', 'cld_geo_thick', 'cld_temp_acha', 'cld_press_acha', 'supercooled_cloud_fraction',
- 'cld_emiss_acha', 'conv_cloud_fraction', 'cld_reff_dcomp', 'cld_opd_dcomp', 'iwc_dcomp', 'lwc_dcomp']
-
- if satellite == 'GOES16':
- train_params_l1b = ['temp_10_4um_nom', 'temp_11_0um_nom', 'temp_12_0um_nom', 'temp_13_3um_nom', 'temp_3_75um_nom',
- 'temp_6_2um_nom', 'temp_6_7um_nom', 'temp_7_3um_nom', 'temp_8_5um_nom', 'temp_9_7um_nom',
- 'refl_0_47um_nom', 'refl_0_65um_nom', 'refl_0_86um_nom', 'refl_1_38um_nom', 'refl_1_60um_nom']
- # 'refl_2_10um_nom']
- elif satellite == 'H08':
- train_params_l1b = ['temp_10_4um_nom', 'temp_12_0um_nom', 'temp_8_5um_nom', 'temp_3_75um_nom', 'refl_2_10um_nom',
- 'refl_1_60um_nom', 'refl_0_86um_nom', 'refl_0_47um_nom']
- else:
- train_params_l2 = ['cld_height_acha', 'cld_geo_thick', 'cld_temp_acha', 'cld_press_acha', 'supercooled_cloud_fraction',
- 'cld_emiss_acha', 'conv_cloud_fraction', 'cld_reff_acha', 'cld_opd_acha']
-
- if use_dnb is True:
- train_params_l2 = ['cld_height_acha', 'cld_geo_thick', 'cld_temp_acha', 'cld_press_acha', 'supercooled_cloud_fraction',
- 'cld_emiss_acha', 'conv_cloud_fraction', 'cld_reff_dcomp', 'cld_opd_dcomp', 'iwc_dcomp', 'lwc_dcomp']
-
- if satellite == 'GOES16':
- train_params_l1b = ['temp_10_4um_nom', 'temp_11_0um_nom', 'temp_12_0um_nom', 'temp_13_3um_nom', 'temp_3_75um_nom',
- 'temp_6_2um_nom', 'temp_6_7um_nom', 'temp_7_3um_nom', 'temp_8_5um_nom', 'temp_9_7um_nom']
- elif satellite == 'H08':
- train_params_l1b = ['temp_10_4um_nom', 'temp_12_0um_nom', 'temp_8_5um_nom', 'temp_3_75um_nom']
-
- if l1b_andor_l2 == 'both':
- train_params = train_params_l1b + train_params_l2
- elif l1b_andor_l2 == 'l1b':
- train_params = train_params_l1b
- elif l1b_andor_l2 == 'l2':
- train_params = train_params_l2
-
- return train_params, train_params_l1b, train_params_l2
--
GitLab