PVDA adapter to produce himawari CMIP-equivalent

Path, Value, Dimensions, Attributes approach
see goesr.rockfall.nc_write et al
requires much testing, but once it works we can use goesr SCMI generator with minor modifications

himawari/HimawariScene.py

@@ -87,6 +87,16 @@ class HimawariScene(object):
         name = ffi.string(meta.satellite_name).decode('UTF-8').strip(' \0')
         return name
 
+    @property
+    def processing_center_name(self):
+        meta = self.metadata
+        return ffi.string(meta.processing_center_name).decode('UTF-8').strip(' \0')
+
+    @property
+    def observation_area(self):
+        meta = self.metadata
+        return ffi.string(meta.observation_area).decode('UTF-8').strip(' \0')
+
     def __init__(self, path, sentinel_invalid=None, sentinel_outside_scan=None,
                  downsample_to_band=0, nav_options=DEFAULT_NAVOPT):
         sentinels = ffi.new('struct HimawariSentinels *')

himawari/pug_cmip_pvda.py

+#!/usr/bin/env python
+# -*- coding: utf-8 -*-
+"""
+PURPOSE
+Generate path-value-dimensions-attributes frames (PVDA) for use with goesr.rockfall
+This can be used with nc_write consumer to generate CMIP PUG-compatible Himawari files
+
+REFERENCES
+
+REQUIRES
+
+:author: R.K.Garcia <rkgarcia@wisc.edu>
+:copyright: 2017 by University of Wisconsin Regents, see AUTHORS for more details
+:license: GPLv3, see LICENSE for more details
+"""
+import os, sys
+import logging, unittest, re
+import numpy as np
+from datetime import datetime, timedelta
+from collections import OrderedDict
+from himawari.HimawariScene import HimawariScene
+
+LOG = logging.getLogger(__name__)
+
+DEFAULT_CMIP_NAME = 'CMI'
+DEFAULT_CMIP_PROJECTION_NAME = 'goes_imager_projection'
+DEFAULT_BAND_DIM_NAME = 'band'
+DEFAULT_YX_BOUNDS_NAME = 'number_of_image_bounds'
+DEFAULT_TIME_BOUNDS_NAME = 'number_of_time_bounds'
+
+class HimawariAHIasCMIP(object):
+    """
+    yield path, value, dimensions, attributes frames used by goesr.cmi_changer
+    """
+    path = None
+    _hs = None   # HimawariScene object
+    _as_radiances = False  # whether ot use bt/refl or rad
+    _data_name = False
+    _projection_name = DEFAULT_CMIP_PROJECTION_NAME
+    _when = None
+
+    def p(self, *strings):
+        return (self.path,) + tuple(strings)
+
+    def d(self, *dimnames):
+        """
+        return dimension dictionary subset
+        :param dimnames: names of dimensions in order
+        :return:
+        """
+        D = self._pug_global_dims
+        return OrderedDict((k,D[k]) for k in dimnames)
+
+    def __init__(self, path, as_radiances=False, data_name=None):
+        self.path = path
+        self._hs = HimawariScene(path)
+        self._as_radiances = as_radiances
+        self._data_name = data_name or (self.bt_or_refl if not as_radiances else 'rad')
+        self._when = datetime.utcnow()
+
+    @property
+    def bt_or_refl(self):
+        return {'albedo': 'refl', 'brightness_temp': 'bt'}[self._hs.kind]
+
+    @property
+    def long_name(self):
+        if not self._as_radiances:
+            data_kind = {'albedo': 'Reflectance', 'brightness_temp': 'Brightness Temperature'}[self._hs.kind]
+        else:
+            data_kind = 'Radiances'
+        return 'Himawari ' + data_kind
+
+    @property
+    def standard_name(self):
+        if not self._as_radiances:
+            data_kind = {'albedo': 'toa_lambertian_equivalent_albedo_multiplied_by_cosine_solar_zenith_angle',
+                         'brightness_temp': 'toa_brightness_temperature'}[self._hs.kind]
+        else:
+            data_kind = 'toa_outgoing_radiance_per_unit_wavenumber'
+        return data_kind
+
+    @property
+    def units(self):
+        if not self._as_radiances:
+            units = {'albedo': '1',
+                     'brightness_temp': 'K'}[self._hs.kind]
+        else:
+            units = 'mW m-2 sr-1 (cm-1)-1'
+        return units
+
+    @property
+    def data_attrs(self):
+        """
+		Rad:_FillValue = 4095s ;
+		Rad:long_name = "ABI L1b Radiances" ;
+		Rad:standard_name = "toa_outgoing_radiance_per_unit_wavenumber" ;
+		Rad:_Unsigned = "true" ;
+		Rad:sensor_band_bit_depth = 12b ;
+		Rad:valid_range = 0s, 4094s ;
+		Rad:scale_factor = 0.04572892f ;
+		Rad:add_offset = -1.6443f ;
+		Rad:units = "mW m-2 sr-1 (cm-1)-1" ;
+		Rad:resolution = "y: 0.000056 rad x: 0.000056 rad" ;
+		Rad:coordinates = "band_id band_wavelength t y x" ;
+		Rad:grid_mapping = "goes_imager_projection" ;
+		Rad:cell_methods = "t: point area: point" ;
+		Rad:ancillary_variables = "DQF" ;
+
+		CMI:_FillValue = -1s ;
+		CMI:long_name = "ABI L2+ Cloud and Moisture Imagery brightness temperature" ;
+		CMI:standard_name = "toa_brightness_temperature" ;
+		CMI:_Unsigned = "true" ;
+		CMI:sensor_band_bit_depth = 12b ;
+		CMI:valid_range = 0s, 4095s ;
+		CMI:scale_factor = 0.03931624f ;
+		CMI:add_offset = 173.15f ;
+		CMI:units = "K" ;
+		CMI:resolution = "y: 0.000056 rad x: 0.000056 rad" ;
+		CMI:coordinates = "band_id band_wavelength t y x" ;
+		CMI:grid_mapping = "goes_imager_projection" ;
+		CMI:cell_methods = "t: point area: point" ;
+		CMI:ancillary_variables = "DQF" ;
+
+		CMI:_FillValue = -1s ;
+		CMI:long_name = "ABI L2+ Cloud and Moisture Imagery reflectance factor" ;
+		CMI:standard_name = "toa_lambertian_equivalent_albedo_multiplied_by_cosine_solar_zenith_angle" ;
+		CMI:_Unsigned = "true" ;
+		CMI:sensor_band_bit_depth = 10b ;
+		CMI:valid_range = 0s, 4095s ;
+		CMI:scale_factor = 0.0002442f ;
+		CMI:add_offset = 0.f ;
+		CMI:units = "1" ;
+		CMI:resolution = "y: 0.000028 rad x: 0.000028 rad" ;
+		CMI:coordinates = "band_id band_wavelength t y x" ;
+		CMI:grid_mapping = "goes_imager_projection" ;
+		CMI:cell_methods = "t: point area: point" ;
+		CMI:ancillary_variables = "DQF" ;
+        :return: dictionary of GOES-R PUG-complaint CMIP attributes for data
+        """
+        return {
+            'long_name': self.long_name,
+            'standard_name': self.standard_name,
+            'units': self.units,
+            'grid_mapping': self._projection_name,
+            'coordinates': 'band_id band_wavelength t y x',
+            # 'sensor_band_bit_depth':
+            # 'resolution':
+        }
+
+    @property
+    def data_dims(self):
+        return self.d('y', 'x')
+
+    @staticmethod
+    def radii(nav):
+        # 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
+        try:
+            f = 1.0 / (1.0 - b/a)
+        except ZeroDivisionError as hcast_probably_did_this:
+            f = np.NAN
+        if np.isnan(f) or f==0.0:
+            LOG.warning('invalid projection parameters, hello HimawariCast - using hardcoded values from HSD guide p12')
+            a = 6378137.0
+            b = 6356752.3
+            f = 1.0 / (1.0 - b/a)
+            inverse_flattening = f # 298.2572f ;
+            perspective_point_height = 42164000.0 - a  # 35785831.0 / 1000.0
+            # 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 ;
+        return {
+            'semi_major_axis': a,
+            'semi_minor_axis': b,
+            'perspective_point_height': perspective_point_height,
+            'inverse_flattening': inverse_flattening,
+        }
+
+    @property
+    def projection_attrs(self):
+        """
+        goes_imager_projection:long_name = "GOES-R ABI fixed grid projection" ;
+		goes_imager_projection:grid_mapping_name = "geostationary" ;
+		goes_imager_projection:perspective_point_height = 35786023. ;
+		goes_imager_projection:semi_major_axis = 6378137. ;
+		goes_imager_projection:semi_minor_axis = 6356752.31414 ;
+		goes_imager_projection:inverse_flattening = 298.2572221 ;
+		goes_imager_projection:latitude_of_projection_origin = 0. ;
+		goes_imager_projection:longitude_of_projection_origin = -89.5 ;
+		goes_imager_projection:sweep_angle_axis = "x" ;
+
+        real(c_double) :: sub_lon = 0.0
+        integer(c_int32_t) :: LFAC = 0
+        integer(c_int32_t) :: CFAC = 0
+        real(c_float) :: LOFF = 0.0
+        real(c_float) :: COFF = 0.0
+        real(c_double) :: distance_from_earth_center_to_virtual_satellite = 0.0
+        real(c_double) :: earth_equatorial_radius = 0.0
+        real(c_double) :: earth_polar_radius = 0.0
+        real(c_float) :: ssp_latitude = 0.0
+        real(c_float) :: ssp_longitude = 0.0
+        real(c_double) :: nadir_latitude = 0.0
+        real(c_double) :: nadir_longitude = 0.0
+        real(c_double) :: distance_from_earth_center_to_satellite = 0.0
+
+        :return: dictionary of PUG-compliant CMIP attributes
+        """
+        nav = self._hs.navigation
+        zult = {
+            'long_name': 'AHI fixed grid projection',
+            'grid_mapping_name': 'geostationary',
+            'latitude_of_projection_origin': 0.0,
+            'longitude_of_projection_origin': nav.sub_lon,
+            'sweep_angle_axis': 'y',
+        }
+        zult.update(self.radii(nav))
+        return zult
+
+
+    @property
+    def data(self):
+        kind = self.bt_or_refl
+        if self._as_radiances:
+            data = self._hs.radiances()
+        elif kind == 'bt':
+            data = self._hs.brightnessTemps()
+        elif kind == 'refl':
+            data = self._hs.albedos()
+        return data
+
+
+    @property
+    def resolution_nadir_meters(self):
+        meta = self.hs.metadata
+        lres, cres = meta.lines_res_meters, meta.columns_res_meters
+        res = min(lres, cres)
+        if res <= 0:
+            raise ValueError('resolution cannot be zero - HimawariCast data?')
+        return res
+
+    @property
+    def satellite_name(self):
+        return self._hs.satellite_name
+
+    @property
+    def orbital_slot(self):
+        return self._hs.satellite_name
+
+    @property
+    def platform_id(self):
+        sat = self._hs.satellite_name
+        # convert to 'H8', 'H9'
+        return ''.join(re.match(r'([A-Z]).*?(\d+)', sat).groups())
+
+    @property
+    def scene_id(self):
+        return self._hs.observation_area
+        # FUTURE: translate from FLDK et al to Full Disk et al?
+
+    @property
+    def datetime_range(self):
+        meta = self._hs.metadata
+        when = datetime(meta.year, meta.month, meta.day, meta.hour, meta.minute, meta.second, meta.microsecond)
+        _, line_times = self._hs.line_times
+        min_sec_offset = np.nanmin(line_times)
+        max_sec_offset = np.nanmax(line_times)
+        return when + timedelta(seconds=min_sec_offset), when + timedelta(seconds=max_sec_offset)
+
+    @property
+    def _pug_global_attrs(self):
+        UTCFMT = '%Y-%m-%dT%H:%M:%SZ'
+        start,stop = self.datetime_range
+        zult = dict(
+            naming_authority="gov.nesdis.noaa",
+            Conventions = "CF-1.7",
+            Metadata_Conventions = "Unidata Dataset Discovery v1.0",
+            standard_name_vocabulary = "CF Standard Name Table (v25, 05 July 2013)",
+            institution = None ,
+            project = None,
+            production_site = None,
+            production_environment = None,
+            spatial_resolution = "%.1fkm at nadir" % (np.round(self.resolution_nadir_meters / 100.0) / 10.0),
+            orbital_slot = self.orbital_slot,
+            platform_ID = self.platform_id,
+            instrument_type = "%s Advanced Himawari Imager" % self.satellite_name,
+            scene_id = self.scene_id,
+            instrument_ID = "AHI",
+            dataset_name = os.path.split(self.path)[-1],
+            iso_series_metadata_id = None,
+            title = "AHI L2 Cloud and Moisture Imagery",
+            summary = None ,
+            keywords = None ,
+            keywords_vocabulary = None ,
+            license = None, # "Unclassified data.  Access is restricted to approved users only.",
+            processing_level = "National Aeronautics and Space Administration (NASA) L2",
+            date_created = self._when.strftime(UTCFMT),
+            cdm_data_type = "Image",
+            time_coverage_start = start.strftime(UTCFMT),
+            time_coverage_end = stop.strftime(UTCFMT),
+            timeline_id = None,
+            production_data_source = None,
+            id = None,
+        )
+        return dict((k,v) for (k,v) in zult.items() if v is not None)
+
+    @property
+    def _pug_global_dims(self):
+        z = OrderedDict()
+        ey, ex = self._hs.shape
+        z['y'] = ey
+        z['x'] = ex
+        z[DEFAULT_BAND_DIM_NAME] = 1
+        z[DEFAULT_YX_BOUNDS_NAME] = 2
+        return z
+
+    def pvda_yx(self):
+
+        # get the m & b values, plus the actual values
+        _,_,mb = self._hs.coords(unscaled=True)
+        y,x = self._hs.coords()
+
+        def yx_pvda(name, v, m, b):
+            assert(name in {'y', 'x'})
+            a = {
+                '.dtype': np.int16,  # annotation telling goesr.rockfall.nc_write what kind of var to make
+                'units': 'rad',
+                'long_name': "AHI fixed grid projection %s-coordinate" % name,
+                'standard_name': "projection_%s_coordinate" % name,
+                'scale_factor': m,
+                'add_offset': b
+            }
+            # y/x
+            yield self.p(name), v, self.d(name), a
+
+            # y_image/x_image single value
+            v_image = np.nanmean(v)
+            a['long_name'] = "AHI fixed grid projection %s-coordinate center of image" % name,
+            yield self.p('%s_image' % name), v_image, {}, a
+
+            # y_image_bounds
+            vrange = np.array([np.nanmin(v), np.nanmax(v)])
+            nsew = {
+                'y': 'north/south',
+                'x': 'west/east'
+            }[name]
+            a = {}
+            a['long_name'] = "AHI fixed grid projection %s-coordinate %s extent of image" % (name, nsew)
+            yield self.p('%s_image_bounds' % name), vrange, self.d(DEFAULT_YX_BOUNDS_NAME), a
+
+        # y and x variables and their bounds and centerpoints
+        yx_pvda('y', y, mb.my, mb.by)
+        yx_pvda('x', x, mb.mx, mb.bx)
+
+    def pvda_nav(self):
+        nav = self._hs.navigation
+
+        a = dict(
+            long_name = "nominal satellite subpoint latitude (platform latitude)",
+            standard_name = "latitude",
+            units = "degrees_north"
+        )
+        yield self.p('nominal_satellite_subpoint_lat'), np.float32(nav.ssp_latitude), {}, a
+
+        a = dict(
+            long_name="nominal satellite subpoint longitude (platform longitude)",
+            standard_name = "longitude",
+            units = "degrees_east"
+        )
+        yield self.p('nominal_satellite_subpoint_lon'), np.float32(nav.ssp_longitude), {}, a
+
+        a = dict(
+            long_name="nominal satellite height above GRS 80 ellipsoid (platform altitude)",
+            standard_name = "height_above_reference_ellipsoid",
+            units = "km"
+        )
+        nsh = nav.distance_from_earth_center_to_satellite - nav.earth_equatorial_radius
+        assert(nsh >= 1000.0)
+        yield self.p('nominal_satellite_height', np.float32(nsh), {}, a)
+
+        # FIXME geospatial
+        """
+        float geospatial_lat_lon_extent ;
+		geospatial_lat_lon_extent:long_name = "geospatial latitude and longitude references" ;
+		geospatial_lat_lon_extent:geospatial_westbound_longitude = -140.6163f ;
+		geospatial_lat_lon_extent:geospatial_northbound_latitude = 52.76771f ;
+		geospatial_lat_lon_extent:geospatial_eastbound_longitude = -49.17929f ;
+		geospatial_lat_lon_extent:geospatial_southbound_latitude = 14.00016f ;
+		geospatial_lat_lon_extent:geospatial_lat_center = 29.294f ;
+		geospatial_lat_lon_extent:geospatial_lon_center = -91.406f ;
+		geospatial_lat_lon_extent:geospatial_lat_nadir = 0.f ;
+		geospatial_lat_lon_extent:geospatial_lon_nadir = -89.5f ;
+		geospatial_lat_lon_extent:geospatial_lat_units = "degrees_north" ;
+		geospatial_lat_lon_extent:geospatial_lon_units = "degrees_east" ;
+		"""
+
+    def pvda_misc(self):
+        band = np.array([self._hs.band], dtype=np.uint8)
+
+        a = dict(
+            long_name="AHI band number",
+            standard_name = "sensor_band_identifier",
+            units = "1"
+        )
+        yield self.p('band_id'), band, self.d(DEFAULT_BAND_DIM_NAME), a
+
+        a = dict(
+            long_name="AHI band central wavelength",
+            standard_name = "sensor_band_central_radiation_wavelength",
+            units = "um"
+        )
+        bwl_um = self._hs.metadata.central_wavelength / 1e6
+        bwl = np.array([bwl_um], dtype=np.float32)
+        yield self.p('band_id'), band, self.d(DEFAULT_BAND_DIM_NAME), a
+
+
+    def __call__(self):
+        """
+        iterate PVDA frames equivalent to a PUG nc_walk
+        """
+        # root
+        yield self.p(), None, self._pug_global_dims, self._pug_global_attrs
+
+        # nav structure
+        yield self.p(DEFAULT_CMIP_PROJECTION_NAME), np.int32(0), {}, self.projection_attrs
+
+        # y and x variables and friends
+        self.pvda_yx()
+        self.pvda_nav()
+
+        # misc
+        self.pvda_misc()
+
+        # we almost forgot! the data itself
+        v = self.data
+        d = self.data_dims
+        a = self.data_attrs
+        yield self.p(self._data_name), v, d, a
+
+
+
+
+PATH_TEST_DATA = os.environ.get('TEST_DATA', os.path.expanduser("~/Data/test_files/thing.dat"))
+
+
+class tests(unittest.TestCase):
+    def setUp(self):
+        pass
+
+    def test_something(self):
+        pass
+
+
+def _debug(type, value, tb):
+    "enable with sys.excepthook = debug"
+    if not sys.stdin.isatty():
+        sys.__excepthook__(type, value, tb)
+    else:
+        import traceback, pdb
+        traceback.print_exception(type, value, tb)
+        # …then start the debugger in post-mortem mode.
+        pdb.post_mortem(tb)  # more “modern”
+
+
+def main():
+    import argparse
+    parser = argparse.ArgumentParser(
+        description="PURPOSE",
+        epilog="",
+        fromfile_prefix_chars='@')
+    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',
+                        help="enable interactive PDB debugger on exception")
+    parser.add_argument('inputs', nargs='*',
+                        help="input files to process")
+    args = parser.parse_args()
+
+    levels = [logging.ERROR, logging.WARN, logging.INFO, logging.DEBUG]
+    logging.basicConfig(level=levels[min(3, args.verbosity)])
+
+    if args.debug:
+        sys.excepthook = _debug
+
+    if not args.inputs:
+        unittest.main()
+        return 0
+
+    for pn in args.inputs:
+        pass
+
+    return 0
+
+
+if __name__ == '__main__':
+    sys.exit(main())