From 7aac9e0662a76d0b13ff436204641f25fda6c47d Mon Sep 17 00:00:00 2001 From: tomrink <rink@ssec.wisc.edu> Date: Sun, 23 Jan 2022 15:55:20 -0600 Subject: [PATCH] NC4 version of write_icing_file --- modules/util/util.py | 144 +++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 144 insertions(+) diff --git a/modules/util/util.py b/modules/util/util.py index f946e656..09ed0be8 100644 --- a/modules/util/util.py +++ b/modules/util/util.py @@ -9,6 +9,7 @@ from collections import namedtuple import os import h5py import pickle +from netCDF4 import Dataset from util.setup import ancillary_path LatLonTuple = namedtuple('LatLonTuple', ['lat', 'lon']) @@ -741,3 +742,146 @@ def write_icing_file(clvrx_str_time, output_dir, preds_dct, probs_dct, x, y, lon pass h5f_out.close() + + +def write_icing_file_nc4(clvrx_str_time, output_dir, preds_dct, probs_dct, x, y, lons, lats, elems, lines): + outfile_name = output_dir + 'icing_prediction_'+clvrx_str_time+'.nc' + rootgrp = Dataset(outfile_name, 'w', format='NETCDF4') + + dim_0_name = 'x' + dim_1_name = 'y' + + dim_0 = rootgrp.createDimension(dim_0_name, size=x.shape[0]) + dim_1 = rootgrp.createDimension(dim_1_name, size=y.shape[0]) + + 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 + + icing_pred_ds = rootgrp.createVariable('icing_prediction_level_'+flt_level_ranges_str[flvl], 'i2', [dim_1_name, dim_0_name]) + icing_pred_ds.setncattr('coordinates', [dim_1_name, dim_0_name]) + icing_pred_ds.setncattr('grid_mapping', 'Projection') + icing_pred_ds.setncattr('missing', -1) + icing_pred_ds[:,] = preds + + 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 + + icing_prob_ds = rootgrp.createVariable('icing_probability_level_'+flt_level_ranges_str[flvl], 'f4', [dim_1_name, dim_0_name]) + icing_prob_ds.setncattr('coordinates', [dim_1_name, dim_0_name]) + icing_prob_ds.setncattr('grid_mapping', 'Projection') + icing_prob_ds.setncattr('missing', -1.0) + icing_prob_ds[:,] = probs + + 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 + + icing_prob_ds = rootgrp.createVariable('max_icing_probability_column', 'f4', [dim_1_name, dim_0_name]) + icing_prob_ds.setncattr('coordinates', [dim_1_name, dim_0_name]) + icing_prob_ds.setncattr('grid_mapping', 'Projection') + icing_prob_ds.setncattr('missing', -1.0) + icing_prob_ds[:,] = max_prob + + 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 + + icing_pred_ds = rootgrp.createVariable('max_icing_probability_level', 'i2', [dim_1_name, dim_0_name]) + icing_pred_ds.setncattr('coordinates', [dim_1_name, dim_0_name]) + icing_pred_ds.setncattr('grid_mapping', 'Projection') + icing_pred_ds.setncattr('missing', -1) + icing_pred_ds[:,] = max_lvl + + # 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) + + proj_ds = rootgrp.createVariable('Projection', 'b') + proj_ds.setncattr('long_name', 'Himawari Imagery Projection') + proj_ds.setncattr('grid_mapping_name', 'geostationary') + proj_ds.setncattr('sweep_angle_axis', 'y') + proj_ds.setncattr('semi_major_axis', 6378.137) + proj_ds.setncattr('semi_minor_axis', 6356.7523) + proj_ds.setncattr('inverse_flattening', 298.257) + proj_ds.setncattr('perspective_point_height', 35785.863) + proj_ds.setncattr('latitude_of_projection_origin', 0.0) + proj_ds.setncattr('longitude_of_projection_origin', 140.7) + + if x is not None: + x_ds = rootgrp.createVariable(dim_0_name, 'f8', [dim_0_name]) + x_ds.units = 'rad' + x_ds.setncattr('standard_name', 'projection_x_coordinate') + x_ds.setncattr('long_name', 'GOES PUG W-E fixed grid viewing angle') + x_ds.setncattr('scale_factor', 5.58879902955962e-05) + x_ds.setncattr('add_offset', -0.153719917308037) + x_ds[:] = x + + y_ds = rootgrp.createVariable(dim_1_name, 'f8', [dim_1_name]) + y_ds.units = 'rad' + y_ds.setncattr('standard_name', 'projection_y_coordinate') + y_ds.setncattr('long_name', 'GOES PUG W-E fixed grid viewing angle') + y_ds.setncattr('scale_factor', -5.58879902955962e-05) + y_ds.setncattr('add_offset', 0.153719917308037) + 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() \ No newline at end of file -- GitLab