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

unfscking py/ vs himawari/ in progress

parent c2efc0fd
......@@ -359,12 +359,12 @@ class HimawariScene(object):
except AssertionError:
self._repalb = alb = None
print "counts[%d,%d] = %.4f" % (line, column, counts[line,column])
print "rad[%d,%d] = %.4f" % (line, column, rad[line,column])
print("counts[%d,%d] = %.4f" % (line, column, counts[line, column]))
print("rad[%d,%d] = %.4f" % (line, column, rad[line, column]))
if temp is not None:
print "btemp[%d,%d] = %.4f" % (line, column, temp[line,column])
print("btemp[%d,%d] = %.4f" % (line, column, temp[line, column]))
elif alb is not None:
print "alb[%d,%d] = %.4f" % (line, column, alb[line,column])
print("alb[%d,%d] = %.4f" % (line, column, alb[line, column]))
def downsample(data, factor, sentinel=np.NAN, first_row=0, first_col=0):
......@@ -406,13 +406,13 @@ def test_partial_scene(stem, remove_count=2):
while remove_count>0:
remove_count -= 1
gbye = randint(0, len(tolink)-1)
print "dropping %s" % tolink[gbye]
print("dropping %s" % tolink[gbye])
del tolink[gbye]
for apath in sorted(tolink):
_,fn = os.path.split(apath)
lpath = os.path.join(dirname, fn)
os.symlink(apath, lpath)
print "%s -> %s" % (apath,lpath)
print("%s -> %s" % (apath, lpath))
toload = os.path.join(dirname, os.path.split(stem)[-1])
fob = HimawariScene(toload)
rmtree(dirname)
......@@ -623,28 +623,31 @@ TEST_COORDS = [ (7332,3667),
def test_jma_nav(stem):
hime = HimawariScene(stem) if isinstance(stem,str) else stem
print "extracting coords"
print("extracting coords")
sy,sx,scale = hime.coords(unscaled=True)
print "computing geo"
print("computing geo")
geo = hime.geo()
print "extracting counts"
print("extracting counts")
counts = hime.counts()
print "extracting radiances"
print("extracting radiances")
rad = hime.radiances()
lines,columns = hime.shape
coords = [(c,l) for (c,l) in TEST_COORDS if (c<columns) and (l<lines)]
for c,l in coords:
pc,pl = (c-1),(l-1) # convert to offsets
print "(Pix,Lin)(%8.1f,%8.1f) --> (Lon,Lat)(%9.3f,%9.3f)" % (sx[pl,pc],sy[pl,pc], geo.longitude[pl,pc], geo.latitude[pl,pc])
print("(Pix,Lin)(%8.1f,%8.1f) --> (Lon,Lat)(%9.3f,%9.3f)" % (
sx[pl, pc], sy[pl, pc], geo.longitude[pl, pc], geo.latitude[pl, pc]))
for c,l in coords:
pc,pl = (c-1),(l-1) # convert to offsets properly
# pc,pl = (c-0),(l-1) # FIXME convert to offsets with extra WTF
print "hPix=%10.5f hLin=%10.5f pl=%d pc=%d" % (sx[pl,pc],sy[pl,pc],pl,pc)
print "(Lon,Lat)(%9.3f,%9.3f) --> count value = %d" % (geo.longitude[pl,pc], geo.latitude[pl,pc], counts[pl,pc])
print("hPix=%10.5f hLin=%10.5f pl=%d pc=%d" % (sx[pl, pc], sy[pl, pc], pl, pc))
print("(Lon,Lat)(%9.3f,%9.3f) --> count value = %d" % (
geo.longitude[pl, pc], geo.latitude[pl, pc], counts[pl, pc]))
for c,l in coords:
pc,pl = (c-1),(l-1) # convert to offsets properly
# pc,pl = (c-0),(l-1) # FIXME convert to offsets with extra WTF
print "(Lon,Lat)(%9.3f,%9.3f) --> radiance = %7.4f" % (geo.longitude[pl,pc], geo.latitude[pl,pc], rad[pl,pc])
print(
"(Lon,Lat)(%9.3f,%9.3f) --> radiance = %7.4f" % (geo.longitude[pl, pc], geo.latitude[pl, pc], rad[pl, pc]))
return True
......@@ -653,7 +656,7 @@ def test_line_times(stem):
hime = HimawariScene(stem) if isinstance(stem,str) else stem
from functools import reduce
base,lt = hime.line_times
print lt
print(lt)
return(np.all(lt[1:]>=lt[:-1]))
# for a,b in zip(lt[1:], lt[:-1]):
# assert(a>=b)
......@@ -661,11 +664,11 @@ def test_line_times(stem):
def test_all(args): # FUTURE: make this a real test pattern
for fn in (args or ['sample/HS_H08_20130710_0300_B01_FLDK_R10']):
hs = HimawariScene(fn)
print fn, hs.shape
print(fn, hs.shape)
del hs
# check decompression cache
hs = HimawariScene(fn)
print "Satellite: %r" % hs.satellite_name
print("Satellite: %r" % hs.satellite_name)
assert(test_line_times(hs))
assert(test_rad_refl(hs))
assert(test_jma_nav(hs))
......
......@@ -14,8 +14,7 @@ REFERENCES
# FIXME: qualityflag
# FIXME: central wavelength to distinguish green band. consider renaming green band to its own ABI_Scaled_b## or AHI_Scaled?
# FIXME: export GOES y and x scale_factor=1.0 and add_offset=0.0 as double or float fully scaled, using conversion from CGMS y/x to GOES y/x
# FIXME: long_name
# FIXME: long_name
REQUIRES
libHimawari.so compiled for your platform
......@@ -134,8 +133,7 @@ NATIVE_REZ = [0, 10,10, 5, 10, 20, 20, 20, 20, 20, 20, 20, 20, 20, 20, 20, 20]
def FACOFF_matches(d, nav):
if d is None: return False
globals().update(d)
if CFAC != nav.CFAC or LFAC != nav.LFAC or COFF != nav.COFF or LOFF != nav.LOFF: return False
if d['CFAC'] != nav.CFAC or d['LFAC'] != nav.LFAC or d['COFF'] != nav.COFF or d['LOFF'] != nav.LOFF: return False
return True
def native_FACOFF(nav):
......@@ -172,13 +170,14 @@ class NetCDF4_writer(object):
_include_fgf = True
_include_rad = False
row_dim_name, col_dim_name = 'y', 'x'
y_var_name, x_var_name = 'cgms_y', 'cgms_x'
y_var_name, x_var_name = 'y', 'x'
bt_var_name = 'brightness_temp'
rad_var_name = 'radiance'
alb_var_name = 'albedo' # reflectance?
lat_var_name = 'latitude'
lon_var_name = 'longitude'
line_time_var_name = 'line_time_offset'
projection_var_name = 'Projection'
bt = None
rad = None
alb = None
......@@ -208,38 +207,37 @@ class NetCDF4_writer(object):
# FUTURE: skip the uint16 -> float -> uint16 conversion, making sure to preserve missing values
self.rad = self._nc.createVariable(self.rad_var_name, 'i2' if not self.rad_as_float else 'f4', dimensions=(self.row_dim_name, self.col_dim_name), fill_value=self.imissing if not self.rad_as_float else self.missing)
self.rad.coordinates = geo_coords if self._include_geo else fgf_coords
self.rad.grid_mapping = self.projection_var_name
self.rad.units = 'W m-2 sr-1 um-1'
else:
if self._kind == 'brightness_temp':
self.bt = self._nc.createVariable(self.bt_var_name, 'f4', dimensions=(self.row_dim_name, self.col_dim_name), fill_value=self.missing)
self.bt.coordinates = geo_coords if self._include_geo else fgf_coords
self.bt.grid_mapping = self.projection_var_name
self.bt.units = 'K'
elif self._kind == 'albedo':
self.alb = self._nc.createVariable(self.alb_var_name, 'f4', dimensions=(self.row_dim_name, self.col_dim_name), fill_value=self.missing)
self.alb.coordinates = geo_coords if self._include_geo else fgf_coords
self.alb.grid_mapping = self.projection_var_name
self.alb.units = '1'
if self._include_fgf:
self.fgf_y = self._nc.createVariable(self.y_var_name, 'f8', dimensions=(self.row_dim_name,))
self.fgf_y = self._nc.createVariable(self.y_var_name, 'i4', dimensions=(self.row_dim_name,))
self.fgf_y.units='radians'
self.fgf_y.standard_name = "projection_y_coordinate"
self.fgf_y.long_name = "CGMS N/S fixed grid viewing angle (not interchangeable with GOES y)"
self.fgf_y.long_name = "GOES PUG S-N fixed grid viewing angle"
self.fgf_x = self._nc.createVariable(self.x_var_name, 'f8', dimensions=(self.col_dim_name,))
self.fgf_x = self._nc.createVariable(self.x_var_name, 'i4', dimensions=(self.col_dim_name,))
self.fgf_x.units='radians'
self.fgf_x.standard_name = "projection_x_coordinate"
self.fgf_x.long_name = "CGMS E/W fixed grid viewing angle (not interchangeable with GOES x)"
# FUTURE: include compatibility 'y' and 'x', though there's a nonlinear transformation from CGMS to GOES y/x angles.
# This requires that the scale_factor and add_offset are 1.0 and 0.0 respectively,
# which violates some uses that use the line/column unscaled form expected by some applications.
self.fgf_x.long_name = "GOES PUG W-E fixed grid viewing angle"
if self._include_geo:
self.lat = self._nc.createVariable(self.lat_var_name, 'f4', dimensions=(self.row_dim_name, self.col_dim_name), fill_value=self.missing)
self.lat.units = 'degrees_north'
self.lon = self._nc.createVariable(self.lon_var_name, 'f4', dimensions=(self.row_dim_name, self.col_dim_name), fill_value=self.missing)
self.lon.units = 'degrees_east'
self.projection = self._nc.createVariable('Projection', 'c', fill_value='\0')
self.projection = self._nc.createVariable(self.projection_var_name, 'c', fill_value='\0')
self.line_time = self._nc.createVariable(self.line_time_var_name, 'f8', dimensions=(self.row_dim_name,))
self.line_time.units = 'seconds POSIX'
self.line_time.long_name = "POSIX epoch seconds elapsed since base_time for image line"
......@@ -264,9 +262,10 @@ class NetCDF4_writer(object):
def set_fgf(self, fgf, downsample_factor=1):
if self.fgf_y is not None:
self.fgf_y[:] = fgf.y
self.fgf_y.scale_factor = fgf.my * float(downsample_factor)
self.fgf_y.add_offset = fgf.by
self.fgf_y[:] = fgf.y
LOG.warning('negating y scale/offset to conform to northward-positive PUG expectation')
self.fgf_y.scale_factor = -fgf.my * float(downsample_factor)
self.fgf_y.add_offset = -fgf.by
if self.fgf_x is not None:
self.fgf_x[:] = fgf.x
self.fgf_x.scale_factor = fgf.mx * float(downsample_factor)
......@@ -320,30 +319,29 @@ class NetCDF4_writer(object):
def set_projection_attrs(self, nav, facoff = None):
"""
ref: http://www.unidata.ucar.edu/software/thredds/current/netcdf-java/reference/StandardCoordinateTransforms.html
assign projection attributes per GRB standard
"""
Projection = self.projection
Projection.grid_mapping_name = "geostationary"
Projection.long_name = "Himawari Imagery Projection" # "GOES ABI Fixed Grid"
Projection._CoordinateTransformType = "Projection"
cax = "%s %s" % (self.lat_var_name, self.lon_var_name) if self._include_geo else "%s %s" % (self.y_var_name, self.x_var_name)
Projection._CoordinateAxes = cax
Projection.sweep_angle_axis = self.col_dim_name
Projection.long_name = "Himawari Imagery Projection" # "GOES ABI Fixed Grid" for GOES-R
Projection.sweep_angle_axis = self.row_dim_name # "y" for Himawari, "x" for GOES-R
Projection.units = "radians"
# calculate invflat 'f' such that rpol = req - req/invflat
a = nav.earth_equatorial_radius
b = nav.earth_polar_radius
a = nav.earth_equatorial_radius # in km from HSD guide
b = nav.earth_polar_radius # km
try:
f = 1.0 / (1.0 - b/a)
except ZeroDivisionError as hcast_probably_did_this:
f = 0.0
if np.isnan(f) or f==0.0:
LOG.warning('invalid projection parameters, hello HimawariCast - using hardcoded values from Harris sample')
Projection.semi_major_axis = a = 6378169.0 / 1000.0 # m ==> km
Projection.semi_minor_axis = b = 6356583.8 / 1000.0
LOG.warning('invalid projection parameters, hello HimawariCast - using hardcoded values from HSD guide p12')
Projection.semi_major_axis = a = 6378.1370 # m ==> km
Projection.semi_minor_axis = b = 6356.7523
f = 1.0 / (1.0 - b/a)
Projection.inverse_flattening = np.float32(f) # 298.2572f ;
Projection.perspective_point_height = 35785831.0 / 1000.0
Projection.perspective_point_height = 42164.0 - a # 35785831.0 / 1000.0
Projection.description = "HimawariCast nominal projection values"
else:
Projection.semi_major_axis = nav.earth_equatorial_radius # 6378.137f ;
......@@ -378,6 +376,9 @@ class NetCDF4_writer(object):
self._nc.distance_from_earth_center_to_satellite = nav.distance_from_earth_center_to_satellite
self._nc.creator = "hsd2nc.py"
self._nc.creation_time = datetime.utcnow().strftime('%Y-%m-%dT%H:%M:%S')
self._nc.Conventions = "CF-1.7"
self._nc.cdm_data_type = "Image"
self._nc.standard_name_vocabulary = "CF Standard Name Table (v25, 05 July 2013)"
def set_time_offsets(self, start_time, time_offsets):
self.line_time.base_time = start_time
......@@ -410,7 +411,8 @@ def hsd2nc(scene, filename='test.nc', include_rad=True,
scene_kind = scene.kind
start = meta.start_time
print "start_time is %04d-%02d-%02dT%02d:%02d:%02d = %d" % (start.year, start.month, start.day, start.hour, start.minute, start.second, start.unix_time)
print("start_time is %04d-%02d-%02dT%02d:%02d:%02d = %d" % (
start.year, start.month, start.day, start.hour, start.minute, start.second, start.unix_time))
downsample_factor = 1 if downsample_factor<=0 else int(downsample_factor)
stride = 1 if stride<=0 else stride
......@@ -446,7 +448,7 @@ def hsd2nc(scene, filename='test.nc', include_rad=True,
ts = time.time()
zult = himawari.downsample(np.ma.fix_invalid(q, copy=False, fill_value=np.NAN), factor=downsample_factor)
te = time.time()
print "%.3f seconds elapsed during downsampling" % (te-ts)
print("%.3f seconds elapsed during downsampling" % (te-ts))
return np.ma.masked_array(zult, np.isnan(zult))
qkw = {
......@@ -533,7 +535,7 @@ def test_quicklook_radiance(filename='test.nc', var='radiance'):
try:
qlc.map_image('test.png', swath=data, lat=lat, lon=lon)
except:
print "quicklook error"
print("quicklook error")
return {'data': data, 'lat': lat, 'lon':lon}
......
......@@ -406,7 +406,7 @@ def test_partial_scene(stem, remove_count=2):
while remove_count>0:
remove_count -= 1
gbye = randint(0, len(tolink)-1)
print "dropping %s" % tolink[gbye]
print("dropping %s" % tolink[gbye])
del tolink[gbye]
for apath in sorted(tolink):
_,fn = os.path.split(apath)
......
Supports Markdown
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