snapshot...

modules/util/util.py

@@ -11,6 +11,7 @@ import h5py
 import pickle
 from netCDF4 import Dataset
 from util.setup import ancillary_path
+from util.geos_nav import get_cf_nav_parameters
 
 LatLonTuple = namedtuple('LatLonTuple', ['lat', 'lon'])
 
@@ -744,7 +745,8 @@ def write_icing_file(clvrx_str_time, output_dir, preds_dct, probs_dct, x, y, lon
     h5f_out.close()
 
 
-def write_icing_file_nc4(clvrx_str_time, output_dir, preds_dct, probs_dct, x, y, lons, lats, elems, lines):
+def write_icing_file_nc4(clvrx_str_time, output_dir, preds_dct, probs_dct,
+                         x, y, lons, lats, elems, lines, satellite='GOES16', domain='CONUS'):
     outfile_name = output_dir + 'icing_prediction_'+clvrx_str_time+'.nc'
     rootgrp = Dataset(outfile_name, 'w', format='NETCDF4')
 
@@ -810,50 +812,41 @@ def write_icing_file_nc4(clvrx_str_time, output_dir, preds_dct, probs_dct, x, y,
     lat_ds.units = 'degrees_north'
     lat_ds[:,] = lats
 
-    # 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)
+    cf_nav_dct = get_cf_nav_parameters(satellite, domain)
+
+    if satellite == 'H08':
+        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', 'icing_prediction_2021-10-07_18:01.nc')
+    proj_ds.setncattr('long_name', long_name)
     proj_ds.setncattr('grid_mapping_name', 'geostationary')
-    proj_ds.setncattr('sweep_angle_axis', 'x')
-    proj_ds.setncattr('semi_major_axis', 6378137.0)
-    proj_ds.setncattr('semi_minor_axis', 6356752.3)
-    proj_ds.setncattr('inverse_flattening', 298.257)
-    proj_ds.setncattr('perspective_point_height', 35786023.0)
-    proj_ds.setncattr('latitude_of_projection_origin', 0.0)
-    proj_ds.setncattr('longitude_of_projection_origin', -75.0)
+    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', '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.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', '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.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: