NC4 version of write_icing_file

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