Something went wrong on our end
Forked from
Ray Garcia / goesr
82 commits behind the upstream repository.
-
Ray Garcia authoredRay Garcia authored
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())