Commit 40bd2bc3 authored by R.K.Garcia's avatar R.K.Garcia
Browse files

adapters to look likes goesr.l1b.PugL1bTools

parent e1f90751
......@@ -17,7 +17,7 @@ import os, sys
import logging, unittest, re
import numpy as np
from datetime import datetime, timedelta
from collections import OrderedDict
from collections import OrderedDict, namedtuple
from himawari.HimawariScene import HimawariScene
LOG = logging.getLogger(__name__)
......@@ -28,6 +28,14 @@ DEFAULT_BAND_DIM_NAME = 'band'
DEFAULT_YX_BOUNDS_NAME = 'number_of_image_bounds'
DEFAULT_TIME_BOUNDS_NAME = 'number_of_time_bounds'
geolocation = namedtuple('geolocation', ['lat', 'lon'])
# http://www.data.jma.go.jp/mscweb/en/himawari89/space_segment/doc/AHI8_performance_test_en.pdf
AHI_BIT_DEPTH = dict(zip(range(1,17), [11, 11, 11, 11,
11, 11, 14, 11,
11, 12, 12, 12,
12, 12, 12, 11]))
class HimawariAHIasCMIP(object):
"""
......@@ -40,18 +48,6 @@ class HimawariAHIasCMIP(object):
_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):
"""
Create a CMIP-compatibility file structure generator around HSD or HCAST data
......@@ -69,6 +65,19 @@ class HimawariAHIasCMIP(object):
self._data_name = data_name
self._when = datetime.utcnow()
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)
@property
def bt_or_refl(self):
return {'albedo': 'refl', 'brightness_temp': 'bt'}[self._hs.kind]
......@@ -158,7 +167,7 @@ class HimawariAHIasCMIP(object):
'units': self.units,
'grid_mapping': self._projection_name,
'coordinates': 'band_id band_wavelength t y x',
# 'sensor_band_bit_depth':
'sensor_band_bit_depth': self.bit_depth
# 'resolution':
}
......@@ -253,6 +262,18 @@ class HimawariAHIasCMIP(object):
raise ValueError('resolution cannot be zero - HimawariCast data?')
return res
@property
def resolution(self):
"""nominal resolution to 100m
"""
rm = self.resolution_nadir_meters
approx = np.round(float(rm) / 100.0)
return int(approx) * 100
@property
def bit_depth(self):
return AHI_BIT_DEPTH[self.band]
@property
def satellite_name(self):
return self._hs.satellite_name
......@@ -261,12 +282,6 @@ class HimawariAHIasCMIP(object):
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
......@@ -428,7 +443,7 @@ class HimawariAHIasCMIP(object):
bwl = np.array([bwl_um], dtype=np.float32)
yield self.p('band_id'), band, self.d(DEFAULT_BAND_DIM_NAME), a
def __iter__(self):
def walk(self):
"""
iterate PVDA frames equivalent to a PUG nc_walk
"""
......@@ -455,6 +470,157 @@ class HimawariAHIasCMIP(object):
a = self.data_attrs
yield self.p(self._data_name), v, d, a
#
# compatibility layer with goesr.l1b.PugL1bTools
#
@property
def mode(self):
return 0
@property
def bit_depth(self):
@property
def band(self):
return self._hs.band
@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 shape(self):
return self._hs.shape
@staticmethod
def _proj4_params(self, longitude_of_projection_origin=None, perspective_point_height=None, semi_major_axis=None,
semi_minor_axis=None, x_0=0, y_0=0, sweep_angle_axis='x', **nav):
return (('proj', 'geos'),
('lon_0', longitude_of_projection_origin),
('h', perspective_point_height),
('a', semi_major_axis),
('b', semi_minor_axis),
('x_0', x_0),
('y_0', y_0),
('sweep', sweep_angle_axis),
('units', 'm'),
# ('ellps', 'GRS80'), # redundant given +a and +b, PUG states GRS80 and not WGS84
)
@property
def proj4_params(self):
return self._proj4_params(**self.projection_attrs)
@property
def proj4_string(self):
return ' '.join(('+%s=%s' % q) for q in self.proj4_params)
@property
def sched_time(self):
"""
nominal observation schedule time for this image
:return:
"""
start, end = self.datetime_range
mid = start + ((end - start) / 2)
meta = self._hs.metadata
hhmm = meta.observation_timeline
hh, mm = int(hhmm / 100), hhmm % 100
return datetime(mid.year, mid.month, mid.day, hh, mm)
@property
def bt(self):
if 'bt' is self.bt_or_refl:
return self._hs.brightnessTemps()
LOG.warning('cannot request bt from non-emissive band')
return None
@property
def refl(self):
if 'refl' is self.bt_or_refl:
return self._hs.albedos()
LOG.warning('cannot request refl from non-reflectance band')
return None
@property
def rad(self):
return self._hs.radiances()
@property
def y(self):
y, x = self._hs.coords()
return y
@property
def x(self):
y, x = self._hs.coords()
return x
@property
def yx_center(self):
y, x = self._hs.coords()
return np.nanmean(y), np.nanmean(x)
@property
def perspective_point_height(self):
nav = self.projection_attrs
return float(nav['perspective_point_height'])
@property
def proj_y(self):
"""
projection coordinate as cartesian nadir-meters as used by PROJ.4
ref http://proj4.org/projections/geos.html
scanning_angle (radians) = projection_coordinate / h
"""
return np.require(self.y, dtype=np.float64) * self.perspective_point_height
@property
def proj_x(self):
"""
projection coordinate as cartesian nadir-meters as used by PROJ.4
ref http://proj4.org/projections/geos.html
scanning_angle (radians) = projection_coordinate / h
"""
return np.require(self.x, dtype=np.float64) * self.perspective_point_height
@property
def shape_m(self):
"""
:return: (y, x) nadir-meters nominal size of a cell
"""
_, _, mb = self._hs.coords(unscaled=True)
h = self.perspective_point_height
return mb.my * h, mb.mx * h
cell_size = shape_m
@property
def origin(self):
"""
:return: (y, x) nadir-meters nominal size of a cell
"""
h = self.perspective_point_height
yo, xo = h * self.y[0], h * self.x[0]
return yo, xo
@property
def geo(self):
from pyproj import Proj
proj = Proj(projparams=dict(self.proj4_params))
y_m = self.proj_y # nadir-meter
x_m = self.proj_x
w,h = len(x_m), len(y_m)
y = np.broadcast_to(y_m, (w, h)).transpose()
x = np.broadcast_to(x_m, (h, w))
lon, lat = proj(x,y, inverse=True)
return geolocation(lon=lon, lat=lat)
PATH_TEST_DATA = os.environ.get('TEST_DATA', os.path.expanduser("~/Data/himawari/sample_bz2/HS_H08_20130710_0300_B16"))
......@@ -467,7 +633,7 @@ class tests(unittest.TestCase):
from goesr.rockfall import debug
chimp = HimawariAHIasCMIP(PATH_TEST_DATA)
n_frames = 0
for stuff in debug(chimp):
for stuff in debug(chimp.walk()):
n_frames += 1
print("%d frames" % n_frames)
......
......@@ -65,6 +65,9 @@ setup(name='himawari',
# package_data={'himawari': ['data/TD/leapsec.dat']},
zip_safe=False,
install_requires=install_requires,
extras_require = {
'scmi': ['goesr'], # https://gitlab.ssec.wisc.edu/rayg/goesr
},
ext_modules = [_himawari],
# entry_points={
# 'console_scripts':
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment