snapshot...

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()