as-radiances -R option to CMIP writer

himawari/pug_cmip_pvda.py

@@ -23,6 +23,7 @@ from himawari.HimawariScene import HimawariScene
 LOG = logging.getLogger(__name__)
 
 DEFAULT_CMIP_NAME = 'CMI'
+DEFAULT_CMIP_RAD_NAME = 'Rad'
 DEFAULT_CMIP_PROJECTION_NAME = 'goes_imager_projection'
 DEFAULT_BAND_DIM_NAME = 'band'
 DEFAULT_YX_BOUNDS_NAME = 'number_of_image_bounds'
@@ -61,11 +62,10 @@ class HimawariAHIasCMIP(object):
         if data_name is True:
             data_name = 'rad' if as_radiances else self.bt_or_refl
         if not data_name:
-            data_name = DEFAULT_CMIP_NAME
+            data_name = DEFAULT_CMIP_NAME if not as_radiances else DEFAULT_CMIP_RAD_NAME
         self._data_name = data_name
         self._when = datetime.utcnow()
 
-
     def p(self, *strings):
         return (self.path,) + tuple(strings)
 
@@ -177,6 +177,10 @@ class HimawariAHIasCMIP(object):
 
     @staticmethod
     def radii(nav):
+        """
+        provide miscellaneous navigation radii in meters, substituting unspecified values for HimawariCast
+        :return: dictionary of nav parameters
+        """
         # calculate invflat 'f' such that rpol = req - req/invflat
         a = nav.earth_equatorial_radius * 1000.0  # km from HSD guide -> m
         b = nav.earth_polar_radius * 1000.0  # km -> m
@@ -194,7 +198,7 @@ class HimawariAHIasCMIP(object):
             # description = "HimawariCast nominal projection values"
         else:
             inverse_flattening = f  # 298.2572f ;
-            perspective_point_height = nav.distance_from_earth_center_to_virtual_satellite - nav.earth_equatorial_radius  # 35786.03f ;
+            perspective_point_height = 1000.0 * (nav.distance_from_earth_center_to_virtual_satellite - nav.earth_equatorial_radius)  # 35786.03f ;
         return {
             'semi_major_axis': a,
             'semi_minor_axis': b,
@@ -272,7 +276,7 @@ class HimawariAHIasCMIP(object):
 
     @property
     def bit_depth(self):
-        return AHI_BIT_DEPTH[self.band]
+        return np.uint8(AHI_BIT_DEPTH[self.band])
 
     @property
     def satellite_name(self):
@@ -288,7 +292,7 @@ class HimawariAHIasCMIP(object):
         # FUTURE: translate from FLDK et al to Full Disk et al?
 
     @property
-    def datetime_range(self):
+    def time_span(self):
         meta = self._hs.metadata
         when = lambda davros: datetime(davros.year, davros.month, davros.day, davros.hour, davros.minute, davros.second,
                                        davros.microsecond)
@@ -297,7 +301,7 @@ class HimawariAHIasCMIP(object):
     @property
     def _pug_global_attrs(self):
         UTCFMT = '%Y-%m-%dT%H:%M:%SZ'
-        start, stop = self.datetime_range
+        start, stop = self.time_span
         zult = dict(
             naming_authority="gov.nesdis.noaa",
             Conventions="CF-1.7",
@@ -344,8 +348,8 @@ class HimawariAHIasCMIP(object):
     def pvda_yx(self):
 
         # get the m & b values, plus the actual values
-        _, _, mb = self._hs.coords(unscaled=True)
-        y, x = self._hs.coords()
+        _, _, mb = self._hs.coords(lines=1, columns=1, unscaled=True)
+        y, x = self.y, self.x
 
         def yx_pvda(name, v, m, b):
             assert (name in {'y', 'x'})
@@ -478,10 +482,6 @@ class HimawariAHIasCMIP(object):
     def mode(self):
         return 0
 
-    @property
-    def bit_depth(self):
-
-
     @property
     def band(self):
         return self._hs.band
@@ -525,7 +525,7 @@ class HimawariAHIasCMIP(object):
         nominal observation schedule time for this image
         :return:
         """
-        start, end = self.datetime_range
+        start, end = self.time_span
         mid = start + ((end - start) / 2)
         meta = self._hs.metadata
         hhmm = meta.observation_timeline
@@ -552,17 +552,17 @@ class HimawariAHIasCMIP(object):
 
     @property
     def y(self):
-        y, x = self._hs.coords()
-        return y
+        y, _ = self._hs.coords(columns=1)
+        return y.squeeze()
 
     @property
     def x(self):
-        y, x = self._hs.coords()
-        return x
+        _, x = self._hs.coords(lines=1)
+        return x.squeeze()
 
     @property
     def yx_center(self):
-        y, x = self._hs.coords()
+        y, x = self.y, self.x
         return np.nanmean(y), np.nanmean(x)
 
     @property
@@ -637,6 +637,18 @@ class tests(unittest.TestCase):
             n_frames += 1
         print("%d frames" % n_frames)
 
+def write_netcdf(path, as_radiances=False):
+    from goesr.rockfall import rename, nc_write
+    fnout = os.path.abspath(os.path.split(path)[-1] + '.nc')
+    lething = HimawariAHIasCMIP(path, as_radiances=as_radiances)
+    for fnwrit in nc_write(
+                            rename(
+                                lething.walk(),
+                                fnout
+                            )
+                        ):
+        print("> " + fnwrit)
+
 
 def _debug(type, value, tb):
     "enable with sys.excepthook = debug"
@@ -655,6 +667,8 @@ def main():
         description="PURPOSE",
         epilog="",
         fromfile_prefix_chars='@')
+    parser.add_argument('-R', '--radiances', dest='rad', action='store_true',
+                        help="write radiances (Rad aka L1b) instead of Brightness Temp or Reflectance (CMI)")
     parser.add_argument('-v', '--verbose', dest='verbosity', action="count", default=0,
                         help='each occurrence increases verbosity 1 level through ERROR-WARNING-INFO-DEBUG')
     parser.add_argument('-d', '--debug', dest='debug', action='store_true',
@@ -674,7 +688,7 @@ def main():
         return 0
 
     for pn in args.inputs:
-        pass
+        write_netcdf(pn, args.rad)
 
     return 0