Skip to content
Snippets Groups Projects
Forked from Ray Garcia / goesr
82 commits behind the upstream repository.
cmi_changer.py 47.54 KiB
#!/usr/bin/env python
# -*- coding: utf-8 -*-
"""
cmi_changer.py
==============

PURPOSE
Write S/CMI data from PUG L1b 

:author: R.K.Garcia <rayg@ssec.wisc.edu>
:copyright: 2017 by University of Wisconsin Regents, see AUTHORS for more details
:license: GPLv3, see LICENSE for more details
"""
import argparse
import logging
import multiprocessing as mp
import os
import re
import socket
import sys
import unittest
from uuid import uuid1, uuid5, UUID
from collections import OrderedDict, defaultdict
from datetime import datetime
from functools import lru_cache

import numpy as np
from pyproj import Proj

import goesr.l1b as l1b
import goesr.rockfall as rf

LOG = logging.getLogger(__name__)

try:
    from goesr._version import VERSION, VERSION_UUID
    VERSION_UUID = UUID(VERSION_UUID)
except ImportError:
    VERSION = "developer build"
    VERSION_UUID = uuid1()
    LOG.warning("could not find goesr._version VERSION string or VERSION_UUID")


def product_uuid(input_name=None):
    """generate a repeatable UUID using the software version UUID plus the input filename/glob,
    on the expectation that those are sufficiently unique
    """
    return uuid1() if input_name is None else uuid5(VERSION_UUID, os.path.split(input_name)[-1])

# zy_xxxx-rrr-Bnn-MnCnn-Tnnn_Gnn_sYYYYDDDhhmmss _cYYYYDDDhhmmss.nc
# OR_HFD-020-B14-M1C07-T055_GH8_s2015181030000_c2015181031543
# ref Table 3.3.4.1.2-1 Sectorized CMI File Naming Convention Fields on NOAA VLAB wiki
FMT_SCMI_NAME = "{environment:1s}{data_type:1s}_{region:s}-{resolution:3s}-B{bits:02d}-M{mode:1d}C{channel:02d}-T{tile:03d}_{satellite:s}_s{scene_time:13s}_c{creation_time:13s}.nc"

# OR_ABI-L2-CMIPC-M3C13_G16_s20171230007204_e20171230009588_c20171230010022.nc
FMT_CMIP_NAME = "{environment:1s}{data_type:1s}_{instrument:3s}-L2-CMIP{region:s}-M{mode:1d}C{channel:02d}_{satellite:s}_s{scene_start_time}_e{scene_end_time}_c{creation_time}.nc"

DEFAULT_SCMI_VAR_NAME = "Sectorized_CMI"
DEFAULT_CMIP_VAR_NAME = "CMI"

FMT_SCMI_TIMESTAMP = '%Y%j%H%M%S'
FMT_CMIP_TIMESTAMP = '%Y%j%H%M%S%f'  # previously '%Y%m%d%H%M%S'; note that format YYYYjjjHHMMSST requires truncation; ref PUG vol5 Table A-1
FMT_CMIP_LONG_TIMESTAMP = "%Y-%m-%dT%H:%M:%S.%f"

# don't pass any variables outside this set thru to SCMI files
SCMI_VARS_WHITELIST = {DEFAULT_SCMI_VAR_NAME, 'DQF', 'y', 'x', 'goes_imager_projection'}
CMIP_VARS_WHITELIST = {DEFAULT_CMIP_VAR_NAME,
                       'DQF',
                       # 'algorithm_dynamic_input_data_container',
                       # 'algorithm_product_version_container',
                       'band_id',
                       'band_wavelength',
                       'earth_sun_distance_anomaly_in_AU',
                       'esun',
                       'geospatial_lat_lon_extent',
                       'goes_imager_projection',
                       'kappa0',
                       'max_reflectance_factor',
                       'mean_reflectance_factor',
                       'min_reflectance_factor',
                       'nominal_satellite_height',
                       'nominal_satellite_subpoint_lat',
                       'nominal_satellite_subpoint_lon',
                       'outlier_pixel_count',
                       'percent_uncorrectable_GRB_errors',
                       'percent_uncorrectable_L0_errors',
                       'planck_bc1',
                       'planck_bc2',
                       'planck_fk1',
                       'planck_fk2',
                       'std_dev_reflectance_factor',
                       't',
                       'time_bounds',
                       'total_number_of_points',
                       'valid_pixel_count',
                       'x',
                       'x_image',
                       'x_image_bounds',
                       'y',
                       'y_image',
                       'y_image_bounds'
                       }

CMIP_ATTRS_WHITELIST = {'naming_authority',
                        # 'Conventions',
                        'Metadata_Conventions',
                        'standard_name_vocabulary',
                        'institution',
                        'project',
                        # 'production_site',
                        # 'production_environment',
                        'spatial_resolution',
                        'orbital_slot',
                        'platform_ID',
                        'instrument_type',
                        'scene_id',
                        'instrument_ID',
                        # 'dataset_name',
                        # 'iso_series_metadata_id',
                        # 'title',
                        # 'summary',
                        # 'keywords',
                        # 'keywords_vocabulary',
                        'license',
                        # 'processing_level',
                        # 'date_created',
                        'cdm_data_type',
                        'time_coverage_start',
                        'time_coverage_end',
                        'timeline_id',
                        'production_data_source',
                        # 'id'
                        }
# NONHERITABLE_L1B_ATTR_NAMES = {'production_site', 'summary', 'dataset', 'title', 'processing_level', 'keywords'}

# region = 'HFD',  # HFD / EFD / WFD / ECONUS / WCONUS / HIREGI / PRREG / AKREGI
REGION_FOR_SCENE = {
    'Full Disk': 'TFD',  # Test Full Disk?
    'CONUS': 'TCONUS',
}

SATELLITE_FOR_PLATFORM = {
    'G16': 'GOES-16',
    'G17': 'GOES-17',
    }

# ref PUG Vol5 L2+ product tables
CMIP_KEYWORDS = {'bt': "SPECTRAL/ENGINEERING > INFRARED WAVELENGTHS > BRIGHTNESS TEMPERATURE",
                 'refl': "ATMOSPHERE > ATMOSPHERIC RADIATION > REFLECTANCE, SPECTRAL/ENGINEERING > VISIBLE WAVELENGTHS > REFLECTANCE"}

CMIP_SUMMARY = {
    'bt': "Single emissive band Cloud and Moisture Imagery Products are digital maps of clouds, moisture and atmospheric windows at IR bands.",
    'refl': "Single reflective band Cloud and Moisture Imagery Products are digital maps of clouds, moisture, and atmospheric windows at visible and near-IR bands."}


def cmi_att_extract():
    from glob import glob
    import netCDF4 as nc4
    import re
    attrology = {}
    for fn in glob('OR*nc'):
        nc = nc4.Dataset(fn)
        ch, = map(int, re.findall('C(\d+)', fn))
        cmi = nc['CMI']
        attrology[ch] = [(a, getattr(cmi, a)) for a in cmi.ncattrs()]
    return attrology

# from libHimawari HCastFile.cxx
# this is included here because having a stable central-wavelength is useful for AWIPS2
# HSD format includes a CWL but HCast uses a nominal value table
AHI_NOMINAL_CENTRAL_WAVELENGTHS_UM = [ 0.46, 0.51, 0.64, 0.86, 1.61, 2.26, 3.89, 6.24, 6.94,
                7.35, 8.59, 9.64, 10.41, 11.24, 12.38, 13.28 ]

# override file-provided central_wavelength values to match AWIPS2:
AWIPS2_CENTRAL_WAVELENGTH = {
    ('G16', 8): 6.19,
}
AWIPS2_CENTRAL_WAVELENGTH.update(dict( (('H8', band), cwl) for (band,cwl) in zip(range(1, 17), AHI_NOMINAL_CENTRAL_WAVELENGTHS_UM) ))

# generated on changeo:~rayg/Workspace/CMIPF-samples/ using PDA OR_ABI-L2-CMIPF-M3C01_G16_s20180012145389_e20180012156156_c20180012156226.nc and its kin

ABI_CMI_ATTS = {
    1: [('_FillValue', -1),
        ('long_name', 'ABI L2+ Cloud and Moisture Imagery reflectance factor'),
        ('standard_name',
         'toa_lambertian_equivalent_albedo_multiplied_by_cosine_solar_zenith_angle'),
        ('_Unsigned', 'true'),
        ('sensor_band_bit_depth', 10),
        ('valid_range', np.array([0, 4095], dtype=np.int16)),
        ('scale_factor', 0.00031746001),
        ('add_offset', 0.0),
        ('units', '1'),
        ('resolution', 'y: 0.000028 rad x: 0.000028 rad'),
        ('coordinates', 'band_id band_wavelength t y x'),
        ('grid_mapping', 'goes_imager_projection'),
        ('cell_methods', 't: point area: point'),
        ('ancillary_variables', 'DQF')],
    2: [('_FillValue', -1),
        ('long_name', 'ABI L2+ Cloud and Moisture Imagery reflectance factor'),
        ('standard_name',
         'toa_lambertian_equivalent_albedo_multiplied_by_cosine_solar_zenith_angle'),
        ('_Unsigned', 'true'),
        ('sensor_band_bit_depth', 12),
        ('valid_range', np.array([0, 4095], dtype=np.int16)),
        ('scale_factor', 0.00031746001),
        ('add_offset', 0.0),
        ('units', '1'),
        ('resolution', 'y: 0.000014 rad x: 0.000014 rad'),
        ('coordinates', 'band_id band_wavelength t y x'),
        ('grid_mapping', 'goes_imager_projection'),
        ('cell_methods', 't: point area: point'),
        ('ancillary_variables', 'DQF')],
    3: [('_FillValue', -1),
        ('long_name', 'ABI L2+ Cloud and Moisture Imagery reflectance factor'),
        ('standard_name',
         'toa_lambertian_equivalent_albedo_multiplied_by_cosine_solar_zenith_angle'),
        ('_Unsigned', 'true'),
        ('sensor_band_bit_depth', 10),
        ('valid_range', np.array([0, 4095], dtype=np.int16)),
        ('scale_factor', 0.00031746001),
        ('add_offset', 0.0),
        ('units', '1'),
        ('resolution', 'y: 0.000028 rad x: 0.000028 rad'),
        ('coordinates', 'band_id band_wavelength t y x'),
        ('grid_mapping', 'goes_imager_projection'),
        ('cell_methods', 't: point area: point'),
        ('ancillary_variables', 'DQF')],
    4: [('_FillValue', -1),
        ('long_name', 'ABI L2+ Cloud and Moisture Imagery reflectance factor'),
        ('standard_name',
         'toa_lambertian_equivalent_albedo_multiplied_by_cosine_solar_zenith_angle'),
        ('_Unsigned', 'true'),
        ('sensor_band_bit_depth', 11),
        ('valid_range', np.array([0, 4095], dtype=np.int16)),
        ('scale_factor', 0.00031746001),
        ('add_offset', 0.0),
        ('units', '1'),
        ('resolution', 'y: 0.000056 rad x: 0.000056 rad'),
        ('coordinates', 'band_id band_wavelength t y x'),
        ('grid_mapping', 'goes_imager_projection'),
        ('cell_methods', 't: point area: point'),
        ('ancillary_variables', 'DQF')],
    5: [('_FillValue', -1),
        ('long_name', 'ABI L2+ Cloud and Moisture Imagery reflectance factor'),
        ('standard_name',
         'toa_lambertian_equivalent_albedo_multiplied_by_cosine_solar_zenith_angle'),
        ('_Unsigned', 'true'),
        ('sensor_band_bit_depth', 10),
        ('valid_range', np.array([0, 4095], dtype=np.int16)),
        ('scale_factor', 0.00031746001),
        ('add_offset', 0.0),
        ('units', '1'),
        ('resolution', 'y: 0.000028 rad x: 0.000028 rad'),
        ('coordinates', 'band_id band_wavelength t y x'),
        ('grid_mapping', 'goes_imager_projection'),
        ('cell_methods', 't: point area: point'),
        ('ancillary_variables', 'DQF')],
    6: [('_FillValue', -1),
        ('long_name', 'ABI L2+ Cloud and Moisture Imagery reflectance factor'),
        ('standard_name',
         'toa_lambertian_equivalent_albedo_multiplied_by_cosine_solar_zenith_angle'),
        ('_Unsigned', 'true'),
        ('sensor_band_bit_depth', 10),
        ('valid_range', np.array([0, 4095], dtype=np.int16)),
        ('scale_factor', 0.00031746001),
        ('add_offset', 0.0),
        ('units', '1'),
        ('resolution', 'y: 0.000056 rad x: 0.000056 rad'),
        ('coordinates', 'band_id band_wavelength t y x'),
        ('grid_mapping', 'goes_imager_projection'),
        ('cell_methods', 't: point area: point'),
        ('ancillary_variables', 'DQF')],
    7: [('_FillValue', -1),
        ('long_name', 'ABI L2+ Cloud and Moisture Imagery brightness temperature'),
        ('standard_name', 'toa_brightness_temperature'),
        ('_Unsigned', 'true'),
        ('sensor_band_bit_depth', 14),
        ('valid_range', np.array([0, 16383], dtype=np.int16)),
        ('scale_factor', 0.01309618),
        ('add_offset', 197.31),
        ('units', 'K'),
        ('resolution', 'y: 0.000056 rad x: 0.000056 rad'),
        ('coordinates', 'band_id band_wavelength t y x'),
        ('grid_mapping', 'goes_imager_projection'),
        ('cell_methods', 't: point area: point'),
        ('ancillary_variables', 'DQF')],
    8: [('_FillValue', -1),
        ('long_name', 'ABI L2+ Cloud and Moisture Imagery brightness temperature'),
        ('standard_name', 'toa_brightness_temperature'),
        ('_Unsigned', 'true'),
        ('sensor_band_bit_depth', 12),
        ('valid_range', np.array([0, 4095], dtype=np.int16)),
        ('scale_factor', 0.042249858),
        ('add_offset', 138.05),
        ('units', 'K'),
        ('resolution', 'y: 0.000056 rad x: 0.000056 rad'),
        ('coordinates', 'band_id band_wavelength t y x'),
        ('grid_mapping', 'goes_imager_projection'),
        ('cell_methods', 't: point area: point'),
        ('ancillary_variables', 'DQF')],
    9: [('_FillValue', -1),
        ('long_name', 'ABI L2+ Cloud and Moisture Imagery brightness temperature'),
        ('standard_name', 'toa_brightness_temperature'),
        ('_Unsigned', 'true'),
        ('sensor_band_bit_depth', 11),
        ('valid_range', np.array([0, 4095], dtype=np.int16)),
        ('scale_factor', 0.042339109),
        ('add_offset', 137.7),
        ('units', 'K'),
        ('resolution', 'y: 0.000056 rad x: 0.000056 rad'),
        ('coordinates', 'band_id band_wavelength t y x'),
        ('grid_mapping', 'goes_imager_projection'),
        ('cell_methods', 't: point area: point'),
        ('ancillary_variables', 'DQF')],
    10: [('_FillValue', -1),
         ('long_name',
          'ABI L2+ Cloud and Moisture Imagery brightness temperature'),
         ('standard_name', 'toa_brightness_temperature'),
         ('_Unsigned', 'true'),
         ('sensor_band_bit_depth', 12),
         ('valid_range', np.array([0, 4095], dtype=np.int16)),
         ('scale_factor', 0.049889188),
         ('add_offset', 126.91),
         ('units', 'K'),
         ('resolution', 'y: 0.000056 rad x: 0.000056 rad'),
         ('coordinates', 'band_id band_wavelength t y x'),
         ('grid_mapping', 'goes_imager_projection'),
         ('cell_methods', 't: point area: point'),
         ('ancillary_variables', 'DQF')],
    11: [('_FillValue', -1),
         ('long_name',
          'ABI L2+ Cloud and Moisture Imagery brightness temperature'),
         ('standard_name', 'toa_brightness_temperature'),
         ('_Unsigned', 'true'),
         ('sensor_band_bit_depth', 12),
         ('valid_range', np.array([0, 4095], dtype=np.int16)),
         ('scale_factor', 0.05216432),
         ('add_offset', 127.69),
         ('units', 'K'),
         ('resolution', 'y: 0.000056 rad x: 0.000056 rad'),
         ('coordinates', 'band_id band_wavelength t y x'),
         ('grid_mapping', 'goes_imager_projection'),
         ('cell_methods', 't: point area: point'),
         ('ancillary_variables', 'DQF')],
    12: [('_FillValue', -1),
         ('long_name',
          'ABI L2+ Cloud and Moisture Imagery brightness temperature'),
         ('standard_name', 'toa_brightness_temperature'),
         ('_Unsigned', 'true'),
         ('sensor_band_bit_depth', 11),
         ('valid_range', np.array([0, 4095], dtype=np.int16)),
         ('scale_factor', 0.047270339),
         ('add_offset', 117.49),
         ('units', 'K'),
         ('resolution', 'y: 0.000056 rad x: 0.000056 rad'),
         ('coordinates', 'band_id band_wavelength t y x'),
         ('grid_mapping', 'goes_imager_projection'),
         ('cell_methods', 't: point area: point'),
         ('ancillary_variables', 'DQF')],
    13: [('_FillValue', -1),
         ('long_name',
          'ABI L2+ Cloud and Moisture Imagery brightness temperature'),
         ('standard_name', 'toa_brightness_temperature'),
         ('_Unsigned', 'true'),
         ('sensor_band_bit_depth', 12),
         ('valid_range', np.array([0, 4095], dtype=np.int16)),
         ('scale_factor', 0.06145332),
         ('add_offset', 89.620003),
         ('units', 'K'),
         ('resolution', 'y: 0.000056 rad x: 0.000056 rad'),
         ('coordinates', 'band_id band_wavelength t y x'),
         ('grid_mapping', 'goes_imager_projection'),
         ('cell_methods', 't: point area: point'),
         ('ancillary_variables', 'DQF')],
    14: [('_FillValue', -1),
         ('long_name',
          'ABI L2+ Cloud and Moisture Imagery brightness temperature'),
         ('standard_name', 'toa_brightness_temperature'),
         ('_Unsigned', 'true'),
         ('sensor_band_bit_depth', 12),
         ('valid_range', np.array([0, 4095], dtype=np.int16)),
         ('scale_factor', 0.059850749),
         ('add_offset', 96.190002),
         ('units', 'K'),
         ('resolution', 'y: 0.000056 rad x: 0.000056 rad'),
         ('coordinates', 'band_id band_wavelength t y x'),
         ('grid_mapping', 'goes_imager_projection'),
         ('cell_methods', 't: point area: point'),
         ('ancillary_variables', 'DQF')],
    15: [('_FillValue', -1),
         ('long_name',
          'ABI L2+ Cloud and Moisture Imagery brightness temperature'),
         ('standard_name', 'toa_brightness_temperature'),
         ('_Unsigned', 'true'),
         ('sensor_band_bit_depth', 12),
         ('valid_range', np.array([0, 4095], dtype=np.int16)),
         ('scale_factor', 0.05956082),
         ('add_offset', 97.379997),
         ('units', 'K'),
         ('resolution', 'y: 0.000056 rad x: 0.000056 rad'),
         ('coordinates', 'band_id band_wavelength t y x'),
         ('grid_mapping', 'goes_imager_projection'),
         ('cell_methods', 't: point area: point'),
         ('ancillary_variables', 'DQF')],
    16: [('_FillValue', -1),
         ('long_name',
          'ABI L2+ Cloud and Moisture Imagery brightness temperature'),
         ('standard_name', 'toa_brightness_temperature'),
         ('_Unsigned', 'true'),
         ('sensor_band_bit_depth', 10),
         ('valid_range', np.array([0, 4095], dtype=np.int16)),
         ('scale_factor', 0.055081531),
         ('add_offset', 92.699997),
         ('units', 'K'),
         ('resolution', 'y: 0.000056 rad x: 0.000056 rad'),
         ('coordinates', 'band_id band_wavelength t y x'),
         ('grid_mapping', 'goes_imager_projection'),
         ('cell_methods', 't: point area: point'),
         ('ancillary_variables', 'DQF')]}



def scmi_product(
        region='XFD',  # HFD / EFD / WFD / ECONUS / WCONUS / HIREGI / PRREG / AKREGI
        resolution='040',  # technically may not be valid to have 4km?
        bits=0,  # 8..14
        mode=1,  # ABI mode
        channel=0,
        **kwargs):  # channel number, 1..16
    # bits = _ahi_bit_depth(channel) if bits==0 else bits
    assert (bits > 0)
    return "{region:s}-{resolution:3s}-B{bits:02d}-M{mode:1d}C{channel:02d}".format(**locals())


def scmi_filename(
        environment='D',  # Integrated Test, Development, Operational
        data_type='T',  # Real-time Playback Simulated Test
        region='XFD',  # HFD / EFD / WFD / ECONUS / WCONUS / HIREGI / PRREG / AKREGI
        resolution='040',  # technically may not be valid to have 4km?
        bits=0,  # 8..14
        mode=1,  # ABI mode
        channel=0,  # channel number, 1..16
        tile=None,  # 001..### upper left to lower right
        satellite='G16',  # G16, G17, ... H8
        scene_time=None,  # datetime object
        creation_time=None,    # now, datetime object
        **kwargs):
    scene_time = scene_time.strftime(FMT_SCMI_TIMESTAMP)
    if creation_time is None:
        creation_time = datetime.utcnow()
    creation_time = creation_time.strftime(FMT_SCMI_TIMESTAMP)
    assert (1 <= channel <= 16)
    assert (bits)
    # bits = _ahi_bit_depth(channel) if bits==0 else bits
    return FMT_SCMI_NAME.format(**locals())


def cmip_filename(
        environment='D',  # Integrated Test, Development, Operational
        data_type='T',  # Real-time Playback Simulated Test
        region='F',  # C / F / M1 / M2 / JP01 / JP02 / etc
        instrument="ABI",  # ABI / AHI / etc
        mode=1,  # ABI mode
        channel=0,  # channel number, 1..16
        satellite='G16',  # G16, G17, ... H8
        scene_start_time=None,  # datetime object, actual observation start
        scene_end_time=None,  # datetime, actual observation end
        creation_time=None,  # datetime
        **kwargs):
    scene_start_time = scene_start_time.strftime(FMT_CMIP_TIMESTAMP)[:14]  # YYYYjjjHHMMSST
    scene_end_time = scene_end_time.strftime(FMT_CMIP_TIMESTAMP)[:14]
    if creation_time is None:
        creation_time = datetime.utcnow()
    creation_time = creation_time.strftime(FMT_CMIP_TIMESTAMP)[:14]
    assert (1 <= channel <= 16)
    # bits = _ahi_bit_depth(channel) if bits==0 else bits
    return FMT_CMIP_NAME.format(**locals())


class axi_scmi_helper(object):
    """
    helper routines for miscellaneous PVDA SCMI filters when generating from PUG L1B files
    """
    as_cmi = None  # whether to use CMI or SCMI convention
    pug = None  # pug object and its attendant netcdf instance
    ny, nx = None, None
    n_tiles = None
    tile_count = None
    creation_time = None
    proj_height = None  # meters, used to translate y and x into nadir-meters used by proj4
    proj = None  # proj4 object
    region = None  # used for filename prefix
    scene = None  # used for .scene_id, since pug .scene_id doesn't work on all HSD/HCAST
    ydim, xdim = 'y', 'x'
    primary_var_name = None  # determined by as_cmi

    def __init__(self, pug: l1b.PugFile, n_tiles=(0, 0), max_pixels=0,
                 environment='D', data_type='T', region=None, scene=None, as_cmi=False):
        self.environment = environment
        self.data_type = data_type
        self.as_cmi = as_cmi
        self.pug = pug
        ny, nx = n_tiles if not as_cmi else (1, 1)
        py, px = self.pug.shape
        def calc_tiles(extent, tiles, mx=max_pixels):
            tiles = int((extent + mx - 1) / mx)
            for p in range(32):  # crudely find an even divisor
                if 0 == (extent % tiles):
                    return tiles
                tiles += 1
            raise ValueError('did not find an even divisor for extent={}'.format(extent))
        if ny <= 0 and max_pixels > 0:
            ny = calc_tiles(py, ny)
        if nx <= 0 and max_pixels > 0:
            nx = calc_tiles(px, nx)
        LOG.info('using ny,nx = {},{}'.format(ny, nx))
        if 0 != (py % ny):
            raise ValueError('cannot evenly divide height {} into {} tiles'.format(py, ny))
        if 0 != (px % nx):
            raise ValueError('cannot evenly divide width {} into {} tiles'.format(px, nx))
        self.ny, self.nx = self.n_tiles = ny, nx
        self.tile_count = ny * nx
        self.creation_time = datetime.utcnow()
        projparams = dict(self.pug.proj4_params)
        self.proj_height = projparams['h']
        self.proj = Proj(projparams=projparams)
        self.scene = scene
        self.region = region or REGION_FOR_SCENE.get(self.scene or self.pug.scene_id, '???')
        if self.region == '???':
            raise ValueError('could not automatically determine region')
        self.primary_var_name = DEFAULT_CMIP_VAR_NAME if as_cmi else DEFAULT_SCMI_VAR_NAME

    @property
    def product_info(self):
        return dict(environment=self.environment,
                    data_type=self.data_type,
                    region=self.region,
                    resolution='%03d' % int(self.pug.resolution / 100),
                    instrument=self.pug.instrument_name,
                    bits=self.pug.bit_depth,
                    mode=self.pug.mode,
                    channel=self.pug.band,
                    satellite=re.match('((?:G|H)\d+)', self.pug.platform_id).group(1),
                    scene_time=self.pug.sched_time,
                    scene_start_time=self.pug.time_span[0],
                    scene_end_time=self.pug.time_span[1],
                    creation_time=self.creation_time,
                    )

    @property
    def product_name(self):
        return scmi_product(**self.product_info) if not self.as_cmi else None # cmip_product(**self.product_info)

    @lru_cache(256)
    def working_filename(self, path, iy, ix, ny, nx):
        """
        used by rf.slice_yx to form the output file name based on tile index 
        note that this is not the final filename, since we may need to count tiles first! 
        :param path: original file path (discarded)
        :param iy: offset in y direction 0..n-1
        :param ix: offset in x direction 0..m-1
        :param ny: total number of tiles in y direction
        :param nx: total number of tiles in x direction
        :return: desired filename to send output to, ultimately used in rf.nc_write
        """
        return '.' + scmi_filename(tile=(self.nx * iy + ix), **self.product_info)

    # def _calculate_slice_numbers(self, mask):
    #     """
    #     use slice_mask to figure out canonical numbers for slices, starting from 1
    #     :return:
    #     """
    #     # convert from boolean to integer and flatten, inverting mask
    #     f = np.array(~mask, dtype=np.int32).ravel()
    #     n = np.cumsum(f)
    #     n = n.reshape(mask.shape)
    #     return n

    def final_filename(self, tile):
        """
        return final filename, now that we know how many slices were empty
        :param tile: tile number, 1..n
        :return:
        """
        return scmi_filename(tile=tile, **self.product_info) if not self.as_cmi else cmip_filename(**self.product_info)

    def annotate_atomic_paths(self, pvda):
        """
        assign final names to files, based on the surviving non-empty tiles
        requires that annotate_abort_on_empty_var is upstream,
        as well as aborted_slice_mask
        :param pvda: 
        :return: 
        """
        for P,V,D,A in pvda:
            sm = A.get('.slice_mask', None)
            sl = A.get('.slice_yx', None)
            sn = A.get('.slice_numbers', None)
            if (sm is not None) and (sl is not None) and (sn is not None):
                # LOG.debug('using slice numbering {}'.format(sn))
                LOG.debug('indexing is {}'.format(sl))
                # LOG.debug('slice mask is {}'.format(sm))
                A = OrderedDict(A)
                iy, ix, ny, nx = sl
                A['.atomic_path'] = p = self.final_filename(sn[iy,ix])
                # also add SCMI global attribute for total tile count
                if not self.as_cmi:
                    A['number_product_tiles'] = np.int16(np.max(sn.ravel()))
                # also set dataset_name to be the final filename (rayg/goesr#21)
                _, fn = os.path.split(p)
                A['dataset_name'] = fn
                LOG.debug('assigning final name {} to {}'.format(p, P[0]))
            yield P, V, D,A

    def sub_cmi_atts(self, atts):
        """
        replace bt/refl attributes with SCMI
        :param atts: incoming variable attributes, which we discard
        :return: substituted attributes for this band, lifted from CMI files
        """
        cmi_atts = ABI_CMI_ATTS[self.pug.band]
        zult = OrderedDict(cmi_atts)
        # atts.update(OrderedDict(atts))
        zult['.dtype'] = np.int16  # creation annotation saying we want scaled integers, see rf._nc_split_creation_atts
        return zult

    def rad_to_bt_or_refl(self, rad):
        LOG.info('transforming radiance for C{:02d}'.format(self.pug.band))
        kind, btrefl, _ = self.pug.convert_radiances(rad[:])
        LOG.info('{} generated from radiance for C{:02d}'.format(kind, self.pug.band))
        return btrefl

    def trim_to_dynamic_range(self, V, A):
        "return value array masked according to bit depth, using scale_factor and add_offset"
        bit_depth = self.pug.bit_depth
        vmin = A['add_offset']
        vmax = vmin + A['scale_factor'] * 2**bit_depth
        return np.ma.masked_array(V, (V < vmin) | (V >= vmax))

    def trim_below_add_offset(self, V, A):
        "mask off values less than add_offset due to _Unsigned"
        return np.ma.masked_array(V, (V < A['add_offset']))

    def pvda_trim_below_add_offset(self, s, var_name=None):
        var_name = var_name or self.primary_var_name
        for p, v, d, a in s:
            if (p[-1] == var_name) and (v is not None):
                v = self.trim_below_add_offset(v, a)
            yield p, v, d, a

    def misc_cmip_fixes(self, pvda):
        """Minor corrections to CMI y and x missing attributes,
        DQF long_name replacement,
        and remove num_star_looks dimension
        """
        for P, V, D, A in pvda:
            if 'num_star_looks' in D:
                D = dict(D)
                del D['num_star_looks']
            if P[-1] in {'y', 'x', 't'}:
                a = dict(A)
                a['axis'] = P[-1].upper()
                yield P, V, D, a
            elif P[-1] == 'DQF':
                a = dict(A)
                a['long_name'] = 'ABI L2+ Cloud and Moisture Imagery reflectance factor data quality flags'
                yield P, V, D, a
            else:
                yield P, V, D, A

    @property
    def transform_table(self):
        """
        variable transforms for l1b to go to SCMI, used by a rf.transform_vars
        this includes setting scale_factor, add_offset, and _dtype override as variable attributes
        :return: dictionary of transforms usable by rockfall.transform_vars 
        """

        if hasattr(self.pug, 'primary_var_name'):  # GOESR
            return {
                # rename and convert Rad to default SCMI variable name, substituting CMI atts for L1B
                self.pug.primary_var_name: (lambda p: p[:-1] + (self.primary_var_name,),  # replace Rad with Sectorized_CMI
                                        self.rad_to_bt_or_refl,  # convert Rad to BT or Refl
                                        lambda dims: dims,  # Dimensions unchanged
                                        self.sub_cmi_atts)  # Attributes replaced with annotated CMI atts for non-CMI input
            }
        else:
            return {}

    def yx_to_latlon(self, y, x):
        """
        convert y and x angles to lat and lon, using proj4
        :param y: radians in projection space
        :param x: 
        :return: lat,lon
        """
        h = self.proj_height
        lon, lat = self.proj(x * h, y * h, inverse=True)
        return lat, lon

    @property
    def product_center_latlon(self):
        y, x = self.pug.yx_center
        return self.yx_to_latlon(y, x)

    @property
    def central_wavelength(self):
        from_file = self.pug.central_wavelength
        key = (self.pug.platform_id, self.pug.band)
        actual = AWIPS2_CENTRAL_WAVELENGTH.get(key, None)
        if actual is not None:
            LOG.info('using central_wavelength override value {} = {} (changed from {})'.format(repr(key), actual, from_file))
            return actual
        return from_file

    def _global_atts_for_tile(self, iy, ix, ny, nx, yang, xang):
        """
        generate global attribute OrderedDict for tile (iy,ix) with its navigation data
        :param iy: tile number 0..ny-1
        :param ix: tile number 0..nx-1
        :param ny: total y slices
        :param nx: total x slices
        :param yang: y fixed grid angle array in radians for this tile 
        :param xang: x fixed grid angle array in radians for this tile
        :return: attributes to add to the tile file 
        """
        # tile dimensions
        th, = yang.shape
        tw, = xang.shape
        # product dimensions
        ph, pw = self.pug.shape
        # nadir dimensions in meters
        dym, dxm = self.pug.shape_m

        # find tile center in y/x angles
        if (th % 2) == 0:  # boundary between pixels, not center of pixel
            ycen = 0.5 * (yang[int(th / 2)] + yang[int(th / 2) + 1])
        else:
            ycen = yang[int(th / 2)]
        if (tw % 2) == 0:
            xcen = 0.5 * (xang[int(tw / 2)] + xang[int(tw / 2) + 1])
        else:
            xcen = xang[int(tw / 2)]

        # get center lat and lon of tile
        latcen, loncen = self.yx_to_latlon(ycen, xcen)
        latprodcen, lonprodcen = self.product_center_latlon
        rez_km = float(self.pug.resolution) / 1000.0

        zult = dict(
            product_name=self.product_name,  # "HFD-010-B11-M1C01", only used for SCMI
            ICD_version="SE-08_7034704_GS_AWIPS_Ext_ICD_RevB.3" if not self.as_cmi else None,
            Conventions="CF-1.6",
        )

        if not self.as_cmi:  # SCMI-specifics
            satellite_id=self.pug.platform_id 
            if satellite_id in SATELLITE_FOR_PLATFORM:
                satellite_id = SATELLITE_FOR_PLATFORM[satellite_id]
            else:
                LOG.warning('Platform/satellite_id not replaced %s' % satellite_id)
           
            s = dict(
                abi_mode=np.int16(self.pug.mode),
                bit_depth=np.int16(self.pug.bit_depth),
                central_wavelength=np.float32(self.central_wavelength),
                channel_id=np.int16(self.pug.band),  # 1,
                product_tile_height=np.int32(th) if not self.as_cmi else None,  # 1100,
                product_tile_width=np.int32(tw) if not self.as_cmi else None,  # 1375,
                production_location=socket.gethostname(),  # "MSC",
                # projection=None,
                periodicity=10,  # 10?
                tile_center_longitude=np.float32(loncen) if not self.as_cmi else None,  # 88.0022078322,
                tile_center_latitude=np.float32(latcen) if not self.as_cmi else None,  # 62.41709885,
                
                satellite_id=satellite_id,
                    
                tile_column_offset=np.int32(ix * tw),  # 2750,
                tile_row_offset=np.int32(iy * th),  # 0,
                start_date_time=self.pug.sched_time.strftime('%Y%j%H%M%S'),  # 2015181030000,  # %Y%j%H%M%S
                pixel_y_size=np.float32(np.abs(dym) / 1000.0),  # km
                pixel_x_size=np.float32(np.abs(dxm) / 1000.0),  # km
                product_rows=np.int32(ph),  # 11000,
                product_columns=np.int32(pw),  # 11000,
                product_center_latitude=np.float32(latprodcen),
                product_center_longitude=np.float32(lonprodcen),  # 140.7,
                # number_product_tiles=np.int16(self.tile_count),  # 76,  added in annotate_atomic_paths based on actual yield!
                satellite_altitude=self.pug.nav['perspective_point_height'] / 1000.0,  # 35785.863
                satellite_longitude=self.pug.nav['longitude_of_projection_origin'],
                # 140.7,  FIXME: ssp_longitude preferable for AHI
                satellite_latitude=self.pug.nav['latitude_of_projection_origin'],
                # FIXME: ssp_latitude preferred for AHI
                source_scene=self.scene or self.pug.scene_id,
                source_spatial_resolution=np.float32(rez_km),  # km
                request_spatial_resolution=np.float32(rez_km),
                title="Sectorized Cloud and Moisture {} Imagery".format(self.scene),
            )
            zult.update(s)

        else:  # CMIP-specifics
            inf = self.product_info
            creation_time = inf['creation_time'].strftime(FMT_CMIP_LONG_TIMESTAMP)[:21] + 'Z'
            product_id = product_uuid(self.pug.path)
            c = dict(
                title="{} L2 Cloud and Moisture Imagery".format(self.pug.instrument_name),
                scene_id=self.scene or self.pug.scene_id,
                production_environment=inf['environment'] + inf['data_type'],
                production_site=socket.gethostname(),
                Conventions='CF-1.7',
                summary=CMIP_SUMMARY[self.pug.bt_or_refl],
                keywords_vocabulary="NASA Global Change Master Directory (GCMD) Earth Science Keywords, Version 7.0.0.0.0",
                created_by="AXI-tools cmi_changer " + VERSION,
                date_created=creation_time,
                # dataset=os.path.split(self.pug.path)[-1],
                processing_level='National Aeronautics and Space Administration (NASA) L2',
                keywords=CMIP_KEYWORDS[self.pug.bt_or_refl],
                id=str(product_id),
                # dataset_name=self.final_filename(0),

                # id
                # institution
                iso_series_metadata_id="8c9e8150-3692-11e3-aa6e-0800200c9a66",  # PUG Vol5
                # keywords
                # keywords_vocabulary
                # license
                # production_data_source
                # production_environment
                # production_site
                # project
                # summary
            )
            zult.update(c)
        return zult

    def append_global_atts(self, pvda, additional_attributes={}):
        """
        given a series of sliced CMI consistent with our data, create global attributes
        assumes that a slice_yx stage is providing .slice_yx annotation in 'A' stream
        :param pvda: path-value-dimensions-attributes sequence
        :yield: pdva + global attributes for SCMI tiles 
        """
        S = defaultdict(dict)  # table of y and x by slice
        for P, V, D, A in pvda:
            # record y and x slices so we can make navigation summary per tile
            if V is not None and len(P) == 2 and P[1] in {self.ydim,
                                                          self.xdim}:  # then this frame is y or x slice for this tile
                key = (P[0],) + tuple(A.get('.slice_yx', (0, 0, 1, 1)))
                S[key][P[1]] = V
            yield P, V, D, A
        # now generate the metadata for each tile as root additions
        mt = OrderedDict()
        for (path, iy, ix, ny, nx), yx in S.items():
            a = self._global_atts_for_tile(iy, ix, ny, nx, yx[self.ydim], yx[self.xdim])
            a.update(additional_attributes)
            LOG.info('assigned additional attributes for ch{:02d} in {}'.format(self.pug.band, path))
            yield (path,), None, mt, a


def axi2scmi(path, n_tiles=(1, 1), max_pixels=0,
             environment='D', data_type='T', region=None, scene=None,
             collate=None, as_cmi=False, autoscale=False, attributes={}, debug=False):
    """
    Convert ABI/AHI to Sectorized CMI for AWIPS2
    Let's make CMIchangas!
    :param path: source file
    :param n_tiles: (ytiles, xtiles)
    :param region: region string to use in filename generation; will attempt autodetection if not provided 
    :param max_pixels: max pixels in any direction; used to determine tile_counts where they are 0 
    :return: 
    """
    # open the file and collect information on it
    if path.endswith('.nc'):
        is_l1b = 'L1b-Rad' in os.path.split(path)[-1]
        LOG.info('attaching GOES-R PUG {} file'.format('L1b' if is_l1b else 'CMI'))
        cmi = l1b.PugL1bTools(path) if is_l1b else l1b.PugCmiTools(path)
        is_cmi = not is_l1b
        is_ahi = False
    else:  # assume we have an AHI HSD or HCAST stem specifying one or more files
        LOG.info('attaching AHI HSD/HCAST')
        try:
            from himawari.ahi2cmi import HimawariAHIasCMIP
        except ImportError as not_available:
            LOG.error('no himwari module available, see https://gitlab.ssec.wisc.edu/rayg/himawari')
            raise
        cmi = HimawariAHIasCMIP(path, as_radiances=False, data_name=DEFAULT_SCMI_VAR_NAME if not as_cmi else DEFAULT_CMIP_VAR_NAME)
        is_cmi = True
        is_ahi = True
    if is_ahi and (scene is None):
        LOG.error("scene identifier is needed for AHI")
        raise ValueError("set scene type when processing AHI data")
    # make a helper which does SCMI-specific calculations and attributes for all the tiles
    helper = axi_scmi_helper(cmi, n_tiles, max_pixels, environment, data_type, region, scene, as_cmi=as_cmi)
    # if n_tiles dimensions are 0s and max_pixels is >0, scmi will calculate the number of tiles
    n_tiles = helper.n_tiles
    # read the input image into a sequence of path, value, dimension, attribute tuples matching PUG conventions
    seq = cmi.walk(as_cmi=as_cmi)
    # convert radiances to bt/refl and rename the variable; also substitute CMI attributes including scale/offset and proper int16 dtype
    if not is_cmi:  # ahi PUG adapter does the conversion on its own unless initialized with as_radiances=True
        seq = rf.transform_vars(seq, **helper.transform_table)
    # discard variables we don't want, we want to do this after the rad transformation since the name is Sectorized_CMI
    seq = rf.discard_unwanted_vars(seq, inclusion_set=CMIP_VARS_WHITELIST if as_cmi else SCMI_VARS_WHITELIST)
    if as_cmi:  # strip a variety of attributes from L1b that are inappropriate for CMIP
        seq = rf.discard_unwanted_attrs(seq, inclusion_set=CMIP_ATTRS_WHITELIST, global_attrs_only=True)
        seq = helper.misc_cmip_fixes(seq)
    # if autoscale is enabled, put autoscale on the CMI variable
    if autoscale or is_ahi:  # we always want to autoscale AHI data due to slight channel differences
        LOG.info('adding autoscale stage for CMI data')
        seq = rf.autoscale(seq, variable_names_matching={DEFAULT_SCMI_VAR_NAME, DEFAULT_CMIP_VAR_NAME}, dtype=np.uint16, dmin=0, dmax=32767)
    else:
        # trim values below guidebook-provided add_offset, to avoid sparkles in night side reflectances
        seq = helper.pvda_trim_below_add_offset(seq)
    # slice into pieces, name the slices to SCMI convention based on A['.slice_yx']; uses default slicing dimensions of y and x
    seq = rf.slice_yx(seq, n_tiles=n_tiles, path_format=helper.working_filename)
    if not as_cmi:
        # cancel any slices that result in empty radiance
        seq = rf.annotate_abort_on_empty_var(seq, helper.primary_var_name)
    # collect navigation information for each tile, then put in additional global attributes before closing the file
    seq = helper.append_global_atts(seq, attributes)
    # note which slices have been aborted, we're going to use that for renaming final tiles
    seq = rf.aborted_slice_mask(seq)
    # tell nc_write to rename the files when we're done, given what's been found about surviving non-empty tiles
    seq = helper.annotate_atomic_paths(seq)
    # add a debug stage at any point by the way to dump everything to screen as it goes
    if debug:
        seq = rf.debug(seq)
    if collate is None:
        # generate all the tiles simply
        for filename in rf.nc_write(seq):  # run the filter stack and generate output files
            print("> {}".format(filename))
    else:  # eat a ton of memory and return all the final frames for alternate uses, sorted by output file path (p[0])
        stack = defaultdict(list)
        for p,v,d,a in seq:
            stack[p[0]].append((p,v,d,a))
        return stack

l1b_makin_the_chimichangas = axi2scmi

RE_ATT = re.compile(r'^\s*(\w+)\s*=\s*(.*)$')

def attribution(args):
    """Given strings in the form name=expression, yield a series of attributename: value
    """
    od = OrderedDict()
    if args is None:
        return od

    for arg in args:
        m = RE_ATT.match(arg)
        if not m:
            raise ValueError("unable to separate {} into name and expression".format(arg))
        name, expr = m.groups()
        try:
            value = eval(expr)
        except:
            LOG.warning("could not evaluate '{}', assuming string attribute".format(expr))
            value = expr
        od[name] = value
    return od

class tests(unittest.TestCase):
    data_file = os.environ.get('TEST_DATA', os.path.expanduser("~/Data/test_files/thing.dat"))

    def setUp(self):
        pass

    def test_something(self):
        pass


def _debug(typ, value, tb):
    "enable with sys.excepthook = debug"
    if not sys.stdin.isatty():
        sys.__excepthook__(typ, value, tb)
    else:
        import traceback, pdb
        traceback.print_exception(typ, value, tb)
        # …then start the debugger in post-mortem mode.
        pdb.post_mortem(tb)  # more “modern”


class chimicherry(object):
    def __init__(self, **kwargs):
        self.kwargs=kwargs
    def __call__(self, *args):
        axi2scmi(*args, **self.kwargs)

def cherrychanga(seq):
    for filename in rf.nc_write(seq):  # run the filter stack and generate output files
        print("> {}".format(filename))

def main():
    parser = argparse.ArgumentParser(
        description="slice and convert ABI L1B, ABI CMIP, AHI HSD or AHI HCAST scenes into SCMI tiles for AWIPS2",
        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")
    # http://docs.python.org/2.7/library/argparse.html#nargs
    # parser.add_argument('--stuff', nargs='5', dest='my_stuff',
    #                    help="one or more random things")
    parser.add_argument('-Y', '--ytiles', type=int, default=0,
                        help="number of tiles in Y dimension")
    parser.add_argument('-X', '--xtiles', type=int, default=0,
                        help="number of tiles in X dimension")
    parser.add_argument('-p', '--pixels', type=int, default=2000,
                        help="max number of pixels in any direction, for tile counts not specified")
    parser.add_argument('-P', '--parallel', type=int, default=1,
                        help="max number of processors to use, default 1; 0 for all-available")
    parser.add_argument('-R', '--region', help='region tag to use in filename, e.g. ECONUS')
    parser.add_argument('-E', '--environment', default='DT',
                        help='environment and data type to use in filename, e.g. DT or OR')
    parser.add_argument('-S', '--scene', help='scene_id tag to use in attributes, e.g. "Full Disk". Default is to use file contents to guess. Needed for AHI.')
    parser.add_argument('-C', '--cmi', dest='as_cmi', action='store_true',
                        help="follow CMI product conventions over SCMI")
    parser.add_argument('-A', '--autoscale', dest='autoscale', action='store_true',
                        help="create scale_factor and add_offset from data content instead of lookup table")
    parser.add_argument('-a', '--attribute', dest='attributes', action='append',
                        help="assign global attribute to all output files given name=expression, e.g. --attribute software_version=\"v0.0.0\"")
    # parser.add_argument('-F', '--fast', action='store_true', default=False,
    #                     help="sacrifice memory for speed")
    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

    additional_attributes = attribution(args.attributes)
    if additional_attributes:
        LOG.info("additional attributes: {!r}".format(additional_attributes))

    chimi = chimicherry(n_tiles=(args.ytiles, args.xtiles), max_pixels=args.pixels,
                        environment=args.environment[0], data_type=args.environment[1],
                        region=args.region, scene=args.scene, autoscale=args.autoscale,
                        as_cmi=args.as_cmi, attributes=additional_attributes,
                        debug=args.debug)

    if args.parallel==1 or len(args.inputs)<=1:
        run_all = None
    elif args.parallel==0:
        multi = mp.Pool()
        run_all = multi.map
        LOG.info('using {} cores'.format(os.cpu_count()))
    else:
        multi = mp.Pool(args.parallel)
        run_all = multi.map
        LOG.info('using {} cores'.format(args.parallel))

    # collection = list(run_all(chimi, args.inputs))
    #
    # if args.fast:
    #     # we haven't actually written any files, let's re-use the pool and write those out
    #     # assumption: we can pass entire pvda sequences back through pipes to master
    #     # desired outcome: run input-parallel, THEN run output-parallel
    #     D = {}
    #     for d in collection:
    #         D.update(d)
    #     run_all(D.values(), cherrychanga)

    if run_all:
        run_all(chimi, args.inputs)
    else:
        for inp in args.inputs:
            chimi(inp)

    return 0


if __name__ == '__main__':
    sys.exit(main())