Skip to content

GitLab

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

reformat pasted tabs out of here

parent 6f5bc9ae
Hide whitespace changes
Inline Side-by-side
himawari/pug_cmip_pvda.py
View file @ bfc6916c
......@@ -28,12 +28,13 @@ DEFAULT_BAND_DIM_NAME = 'band'
DEFAULT_YX_BOUNDS_NAME = 'number_of_image_bounds'
DEFAULT_TIME_BOUNDS_NAME = 'number_of_time_bounds'
class HimawariAHIasCMIP(object):
"""
yield path, value, dimensions, attributes frames used by goesr.cmi_changer
"""
path = None
_hs = None # HimawariScene object
_hs = None # HimawariScene object
_as_radiances = False # whether ot use bt/refl or rad
_data_name = False
_projection_name = DEFAULT_CMIP_PROJECTION_NAME
......@@ -49,7 +50,7 @@ class HimawariAHIasCMIP(object):
:return:
"""
D = self._pug_global_dims
return OrderedDict((k,D[k]) for k in dimnames)
return OrderedDict((k, D[k]) for k in dimnames)
def __init__(self, path, as_radiances=False, data_name=None):
self.path = path
......@@ -91,50 +92,50 @@ class HimawariAHIasCMIP(object):
@property
def data_attrs(self):
"""
Rad:_FillValue = 4095s ;
Rad:long_name = "ABI L1b Radiances" ;
Rad:standard_name = "toa_outgoing_radiance_per_unit_wavenumber" ;
Rad:_Unsigned = "true" ;
Rad:sensor_band_bit_depth = 12b ;
Rad:valid_range = 0s, 4094s ;
Rad:scale_factor = 0.04572892f ;
Rad:add_offset = -1.6443f ;
Rad:units = "mW m-2 sr-1 (cm-1)-1" ;
Rad:resolution = "y: 0.000056 rad x: 0.000056 rad" ;
Rad:coordinates = "band_id band_wavelength t y x" ;
Rad:grid_mapping = "goes_imager_projection" ;
Rad:cell_methods = "t: point area: point" ;
Rad:ancillary_variables = "DQF" ;
CMI:_FillValue = -1s ;
CMI:long_name = "ABI L2+ Cloud and Moisture Imagery brightness temperature" ;
CMI:standard_name = "toa_brightness_temperature" ;
CMI:_Unsigned = "true" ;
CMI:sensor_band_bit_depth = 12b ;
CMI:valid_range = 0s, 4095s ;
CMI:scale_factor = 0.03931624f ;
CMI:add_offset = 173.15f ;
CMI:units = "K" ;
CMI:resolution = "y: 0.000056 rad x: 0.000056 rad" ;
CMI:coordinates = "band_id band_wavelength t y x" ;
CMI:grid_mapping = "goes_imager_projection" ;
CMI:cell_methods = "t: point area: point" ;
CMI:ancillary_variables = "DQF" ;
CMI:_FillValue = -1s ;
CMI:long_name = "ABI L2+ Cloud and Moisture Imagery reflectance factor" ;
CMI:standard_name = "toa_lambertian_equivalent_albedo_multiplied_by_cosine_solar_zenith_angle" ;
CMI:_Unsigned = "true" ;
CMI:sensor_band_bit_depth = 10b ;
CMI:valid_range = 0s, 4095s ;
CMI:scale_factor = 0.0002442f ;
CMI:add_offset = 0.f ;
CMI:units = "1" ;
CMI:resolution = "y: 0.000028 rad x: 0.000028 rad" ;
CMI:coordinates = "band_id band_wavelength t y x" ;
CMI:grid_mapping = "goes_imager_projection" ;
CMI:cell_methods = "t: point area: point" ;
CMI:ancillary_variables = "DQF" ;
Rad:_FillValue = 4095s ;
Rad:long_name = "ABI L1b Radiances" ;
Rad:standard_name = "toa_outgoing_radiance_per_unit_wavenumber" ;
Rad:_Unsigned = "true" ;
Rad:sensor_band_bit_depth = 12b ;
Rad:valid_range = 0s, 4094s ;
Rad:scale_factor = 0.04572892f ;
Rad:add_offset = -1.6443f ;
Rad:units = "mW m-2 sr-1 (cm-1)-1" ;
Rad:resolution = "y: 0.000056 rad x: 0.000056 rad" ;
Rad:coordinates = "band_id band_wavelength t y x" ;
Rad:grid_mapping = "goes_imager_projection" ;
Rad:cell_methods = "t: point area: point" ;
Rad:ancillary_variables = "DQF" ;
CMI:_FillValue = -1s ;
CMI:long_name = "ABI L2+ Cloud and Moisture Imagery brightness temperature" ;
CMI:standard_name = "toa_brightness_temperature" ;
CMI:_Unsigned = "true" ;
CMI:sensor_band_bit_depth = 12b ;
CMI:valid_range = 0s, 4095s ;
CMI:scale_factor = 0.03931624f ;
CMI:add_offset = 173.15f ;
CMI:units = "K" ;
CMI:resolution = "y: 0.000056 rad x: 0.000056 rad" ;
CMI:coordinates = "band_id band_wavelength t y x" ;
CMI:grid_mapping = "goes_imager_projection" ;
CMI:cell_methods = "t: point area: point" ;
CMI:ancillary_variables = "DQF" ;
CMI:_FillValue = -1s ;
CMI:long_name = "ABI L2+ Cloud and Moisture Imagery reflectance factor" ;
CMI:standard_name = "toa_lambertian_equivalent_albedo_multiplied_by_cosine_solar_zenith_angle" ;
CMI:_Unsigned = "true" ;
CMI:sensor_band_bit_depth = 10b ;
CMI:valid_range = 0s, 4095s ;
CMI:scale_factor = 0.0002442f ;
CMI:add_offset = 0.f ;
CMI:units = "1" ;
CMI:resolution = "y: 0.000028 rad x: 0.000028 rad" ;
CMI:coordinates = "band_id band_wavelength t y x" ;
CMI:grid_mapping = "goes_imager_projection" ;
CMI:cell_methods = "t: point area: point" ;
CMI:ancillary_variables = "DQF" ;
:return: dictionary of GOES-R PUG-complaint CMIP attributes for data
"""
return {
......@@ -154,23 +155,23 @@ class HimawariAHIasCMIP(object):
@staticmethod
def radii(nav):
# calculate invflat 'f' such that rpol = req - req/invflat
a = nav.earth_equatorial_radius * 1000.0 # km from HSD guide -> m
b = nav.earth_polar_radius * 1000.0 # km -> m
a = nav.earth_equatorial_radius * 1000.0 # km from HSD guide -> m
b = nav.earth_polar_radius * 1000.0 # km -> m
try:
f = 1.0 / (1.0 - b/a)
f = 1.0 / (1.0 - b / a)
except ZeroDivisionError as hcast_probably_did_this:
f = np.NAN
if np.isnan(f) or f==0.0:
if np.isnan(f) or f == 0.0:
LOG.warning('invalid projection parameters, hello HimawariCast - using hardcoded values from HSD guide p12')
a = 6378137.0
b = 6356752.3
f = 1.0 / (1.0 - b/a)
inverse_flattening = f # 298.2572f ;
f = 1.0 / (1.0 - b / a)
inverse_flattening = f # 298.2572f ;
perspective_point_height = 42164000.0 - a # 35785831.0 / 1000.0
# description = "HimawariCast nominal projection values"
else:
inverse_flattening = f # 298.2572f ;
perspective_point_height = nav.distance_from_earth_center_to_virtual_satellite - nav.earth_equatorial_radius # 35786.03f ;
inverse_flattening = f # 298.2572f ;
perspective_point_height = nav.distance_from_earth_center_to_virtual_satellite - nav.earth_equatorial_radius # 35786.03f ;
return {
'semi_major_axis': a,
'semi_minor_axis': b,
......@@ -182,14 +183,14 @@ class HimawariAHIasCMIP(object):
def projection_attrs(self):
"""
goes_imager_projection:long_name = "GOES-R ABI fixed grid projection" ;
goes_imager_projection:grid_mapping_name = "geostationary" ;
goes_imager_projection:perspective_point_height = 35786023. ;
goes_imager_projection:semi_major_axis = 6378137. ;
goes_imager_projection:semi_minor_axis = 6356752.31414 ;
goes_imager_projection:inverse_flattening = 298.2572221 ;
goes_imager_projection:latitude_of_projection_origin = 0. ;
goes_imager_projection:longitude_of_projection_origin = -89.5 ;
goes_imager_projection:sweep_angle_axis = "x" ;
goes_imager_projection:grid_mapping_name = "geostationary" ;
goes_imager_projection:perspective_point_height = 35786023. ;
goes_imager_projection:semi_major_axis = 6378137. ;
goes_imager_projection:semi_minor_axis = 6356752.31414 ;
goes_imager_projection:inverse_flattening = 298.2572221 ;
goes_imager_projection:latitude_of_projection_origin = 0. ;
goes_imager_projection:longitude_of_projection_origin = -89.5 ;
goes_imager_projection:sweep_angle_axis = "x" ;
real(c_double) :: sub_lon = 0.0
integer(c_int32_t) :: LFAC = 0
......@@ -218,7 +219,6 @@ class HimawariAHIasCMIP(object):
zult.update(self.radii(nav))
return zult
@property
def data(self):
kind = self.bt_or_refl
......@@ -230,7 +230,6 @@ class HimawariAHIasCMIP(object):
data = self._hs.albedos()
return data
@property
def resolution_nadir_meters(self):
meta = self._hs.metadata
......@@ -262,45 +261,46 @@ class HimawariAHIasCMIP(object):
@property
def datetime_range(self):
meta = self._hs.metadata
when = lambda davros: datetime(davros.year, davros.month, davros.day, davros.hour, davros.minute, davros.second, davros.microsecond)
when = lambda davros: datetime(davros.year, davros.month, davros.day, davros.hour, davros.minute, davros.second,
davros.microsecond)
return when(meta.start_time), when(meta.end_time)
@property
def _pug_global_attrs(self):
UTCFMT = '%Y-%m-%dT%H:%M:%SZ'
start,stop = self.datetime_range
start, stop = self.datetime_range
zult = dict(
naming_authority="gov.nesdis.noaa",
Conventions = "CF-1.7",
Metadata_Conventions = "Unidata Dataset Discovery v1.0",
standard_name_vocabulary = "CF Standard Name Table (v25, 05 July 2013)",
institution = None ,
project = None,
production_site = None,
production_environment = None,
spatial_resolution = "%.1fkm at nadir" % (np.round(self.resolution_nadir_meters / 100.0) / 10.0),
orbital_slot = self.orbital_slot,
platform_ID = self.platform_id,
instrument_type = "%s Advanced Himawari Imager" % self.satellite_name,
scene_id = self.scene_id,
instrument_ID = "AHI",
dataset_name = os.path.split(self.path)[-1],
iso_series_metadata_id = None,
title = "AHI L2 Cloud and Moisture Imagery",
summary = None ,
keywords = None ,
keywords_vocabulary = None ,
license = None, # "Unclassified data. Access is restricted to approved users only.",
processing_level = "National Aeronautics and Space Administration (NASA) L2",
date_created = self._when.strftime(UTCFMT),
cdm_data_type = "Image",
time_coverage_start = start.strftime(UTCFMT),
time_coverage_end = stop.strftime(UTCFMT),
timeline_id = None,
production_data_source = None,
id = None,
Conventions="CF-1.7",
Metadata_Conventions="Unidata Dataset Discovery v1.0",
standard_name_vocabulary="CF Standard Name Table (v25, 05 July 2013)",
institution=None,
project=None,
production_site=None,
production_environment=None,
spatial_resolution="%.1fkm at nadir" % (np.round(self.resolution_nadir_meters / 100.0) / 10.0),
orbital_slot=self.orbital_slot,
platform_ID=self.platform_id,
instrument_type="%s Advanced Himawari Imager" % self.satellite_name,
scene_id=self.scene_id,
instrument_ID="AHI",
dataset_name=os.path.split(self.path)[-1],
iso_series_metadata_id=None,
title="AHI L2 Cloud and Moisture Imagery",
summary=None,
keywords=None,
keywords_vocabulary=None,
license=None, # "Unclassified data. Access is restricted to approved users only.",
processing_level="National Aeronautics and Space Administration (NASA) L2",
date_created=self._when.strftime(UTCFMT),
cdm_data_type="Image",
time_coverage_start=start.strftime(UTCFMT),
time_coverage_end=stop.strftime(UTCFMT),
timeline_id=None,
production_data_source=None,
id=None,
)
return dict((k,v) for (k,v) in zult.items() if v is not None)
return dict((k, v) for (k, v) in zult.items() if v is not None)
@property
def _pug_global_dims(self):
......@@ -315,11 +315,11 @@ class HimawariAHIasCMIP(object):
def pvda_yx(self):
# get the m & b values, plus the actual values
_,_,mb = self._hs.coords(unscaled=True)
y,x = self._hs.coords()
_, _, mb = self._hs.coords(unscaled=True)
y, x = self._hs.coords()
def yx_pvda(name, v, m, b):
assert(name in {'y', 'x'})
assert (name in {'y', 'x'})
a = {
'.dtype': np.int16, # annotation telling goesr.rockfall.nc_write what kind of var to make
'units': 'rad',
......@@ -357,64 +357,63 @@ class HimawariAHIasCMIP(object):
nav = self._hs.navigation
a = dict(
long_name = "nominal satellite subpoint latitude (platform latitude)",
standard_name = "latitude",
units = "degrees_north"
long_name="nominal satellite subpoint latitude (platform latitude)",
standard_name="latitude",
units="degrees_north"
)
yield self.p('nominal_satellite_subpoint_lat'), np.float32(nav.ssp_latitude), {}, a
a = dict(
long_name="nominal satellite subpoint longitude (platform longitude)",
standard_name = "longitude",
units = "degrees_east"
standard_name="longitude",
units="degrees_east"
)
yield self.p('nominal_satellite_subpoint_lon'), np.float32(nav.ssp_longitude), {}, a
a = dict(
long_name="nominal satellite height above GRS 80 ellipsoid (platform altitude)",
standard_name = "height_above_reference_ellipsoid",
units = "km"
standard_name="height_above_reference_ellipsoid",
units="km"
)
nsh = nav.distance_from_earth_center_to_satellite - nav.earth_equatorial_radius
assert(nsh >= 1000.0)
assert (nsh >= 1000.0)
yield self.p('nominal_satellite_height'), np.float32(nsh), {}, a
# FIXME geospatial
"""
float geospatial_lat_lon_extent ;
geospatial_lat_lon_extent:long_name = "geospatial latitude and longitude references" ;
geospatial_lat_lon_extent:geospatial_westbound_longitude = -140.6163f ;
geospatial_lat_lon_extent:geospatial_northbound_latitude = 52.76771f ;
geospatial_lat_lon_extent:geospatial_eastbound_longitude = -49.17929f ;
geospatial_lat_lon_extent:geospatial_southbound_latitude = 14.00016f ;
geospatial_lat_lon_extent:geospatial_lat_center = 29.294f ;
geospatial_lat_lon_extent:geospatial_lon_center = -91.406f ;
geospatial_lat_lon_extent:geospatial_lat_nadir = 0.f ;
geospatial_lat_lon_extent:geospatial_lon_nadir = -89.5f ;
geospatial_lat_lon_extent:geospatial_lat_units = "degrees_north" ;
geospatial_lat_lon_extent:geospatial_lon_units = "degrees_east" ;
"""
geospatial_lat_lon_extent:long_name = "geospatial latitude and longitude references" ;
geospatial_lat_lon_extent:geospatial_westbound_longitude = -140.6163f ;
geospatial_lat_lon_extent:geospatial_northbound_latitude = 52.76771f ;
geospatial_lat_lon_extent:geospatial_eastbound_longitude = -49.17929f ;
geospatial_lat_lon_extent:geospatial_southbound_latitude = 14.00016f ;
geospatial_lat_lon_extent:geospatial_lat_center = 29.294f ;
geospatial_lat_lon_extent:geospatial_lon_center = -91.406f ;
geospatial_lat_lon_extent:geospatial_lat_nadir = 0.f ;
geospatial_lat_lon_extent:geospatial_lon_nadir = -89.5f ;
geospatial_lat_lon_extent:geospatial_lat_units = "degrees_north" ;
geospatial_lat_lon_extent:geospatial_lon_units = "degrees_east" ;
"""
def pvda_misc(self):
band = np.array([self._hs.band], dtype=np.uint8)
a = dict(
long_name="AHI band number",
standard_name = "sensor_band_identifier",
units = "1"
standard_name="sensor_band_identifier",
units="1"
)
yield self.p('band_id'), band, self.d(DEFAULT_BAND_DIM_NAME), a
a = dict(
long_name="AHI band central wavelength",
standard_name = "sensor_band_central_radiation_wavelength",
units = "um"
standard_name="sensor_band_central_radiation_wavelength",
units="um"
)
bwl_um = self._hs.metadata.central_wavelength / 1e6
bwl = np.array([bwl_um], dtype=np.float32)
yield self.p('band_id'), band, self.d(DEFAULT_BAND_DIM_NAME), a
def __iter__(self):
"""
iterate PVDA frames equivalent to a PUG nc_walk
......@@ -433,7 +432,7 @@ class HimawariAHIasCMIP(object):
# misc
for each in self.pvda_misc():
assert(len(each)==4)
assert (len(each) == 4)
yield each
# we almost forgot! the data itself
......@@ -443,8 +442,6 @@ class HimawariAHIasCMIP(object):
yield self.p(self._data_name), v, d, a
PATH_TEST_DATA = os.environ.get('TEST_DATA', os.path.expanduser("~/Data/himawari/sample_bz2/HS_H08_20130710_0300_B16"))
......@@ -461,7 +458,6 @@ class tests(unittest.TestCase):
print("%d frames" % n_frames)
def _debug(type, value, tb):
"enable with sys.excepthook = debug"
if not sys.stdin.isatty():
......
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