From 60b9e7b80b31cfa4d67b86e5edbf6eda90dd9905 Mon Sep 17 00:00:00 2001 From: tomrink <rink@ssec.wisc.edu> Date: Fri, 29 Apr 2022 11:32:37 -0500 Subject: [PATCH] snapshot... --- modules/util/util.py | 114 +++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 114 insertions(+) diff --git a/modules/util/util.py b/modules/util/util.py index f781d510..000f6adc 100644 --- a/modules/util/util.py +++ b/modules/util/util.py @@ -1047,3 +1047,117 @@ def write_icing_file_nc4(clvrx_str_time, output_dir, preds_dct, probs_dct, pass rootgrp.close() + + +def write_icing_file_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 = '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 + + # 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']) + + rootgrp.close() -- GitLab