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

Add .proj and .latlon_center to match goesr.l1b updates for SCMI

Technically optional since we only require this for GOESR MESOs
parent 361aeb9c
......@@ -498,34 +498,6 @@ class HimawariAHIasCMIP(object):
bwl = np.array([bwl_um], dtype=np.float32)
yield self.p('band_wavelength'), bwl, self.d(DEFAULT_BAND_DIM_NAME), a
def lat_lon_to_l_c(self, lat, lon, round=True):
"""
convert latitude-longitude all the way back to line-column by reversing projection and y/x scaling
:param lat: latitude value or array in degrees
:param lon: longitude value or array in degrees
:return: line, column array indexed 0..height-1, 0..width-1
"""
from pyproj import Proj
proj = Proj(projparams=dict(self.proj4_params))
# convert to nadir-meters using proj4
xm, ym = proj(lon, lat, inverse=False)
# divide out height to get back to projection view angles
h = self.perspective_point_height
y, x = ym / h, xm / h
_, _, scale = self._hs.coords(lines=1, columns=1, unscaled=True)
# extract scale factor and add offset to back out to line-from-0, column-from-0
# also note Y is inverse between AHI and GOES/PROJ4
ysf, yao = -scale.my, -scale.by
xsf, xao = scale.mx, scale.bx
# reverse scale factor and add offset to get line and column
# note that AHI coordinate system counts from 1! so we subtract 1.0 to get consistency with GOES
line, column = (y - yao) / ysf - 1.0, (x - xao) / xsf - 1.0
# round to nearest integer
if round:
return np.round(line).astype(np.int32), np.round(column).astype(np.int32)
else: # return floating point values
return line, column
def walk(self, extra_data_attrs={}, as_cmi=False):
"""
iterate PVDA frames equivalent to a PUG nc_walk
......@@ -714,9 +686,20 @@ class HimawariAHIasCMIP(object):
return yo, xo
@property
def geo(self):
def latlon_center(self):
proj = self.proj
y, x = self.yx_center
lon, lat = proj(x, y, inverse=True)
return geolocation(lon=lon, lat=lat)
@property
def proj(self):
from pyproj import Proj
proj = Proj(projparams=dict(self.proj4_params))
return Proj(projparams=dict(self.proj4_params))
@property
def geo(self):
proj = self.proj
y_m = self.proj_y # nadir-meter
x_m = self.proj_x
w,h = len(x_m), len(y_m)
......@@ -725,6 +708,33 @@ class HimawariAHIasCMIP(object):
lon, lat = proj(x,y, inverse=True)
return geolocation(lon=lon, lat=lat)
def lat_lon_to_l_c(self, lat, lon, round=True):
"""
convert latitude-longitude all the way back to line-column by reversing projection and y/x scaling
:param lat: latitude value or array in degrees
:param lon: longitude value or array in degrees
:return: line, column array indexed 0..height-1, 0..width-1
"""
proj = self.proj
# convert to nadir-meters using proj4
xm, ym = proj(lon, lat, inverse=False)
# divide out height to get back to projection view angles
h = self.perspective_point_height
y, x = ym / h, xm / h
_, _, scale = self._hs.coords(lines=1, columns=1, unscaled=True)
# extract scale factor and add offset to back out to line-from-0, column-from-0
# also note Y is inverse between AHI and GOES/PROJ4
ysf, yao = -scale.my, -scale.by
xsf, xao = scale.mx, scale.bx
# reverse scale factor and add offset to get line and column
# note that AHI coordinate system counts from 1! so we subtract 1.0 to get consistency with GOES
line, column = (y - yao) / ysf - 1.0, (x - xao) / xsf - 1.0
# round to nearest integer
if round:
return np.round(line).astype(np.int32), np.round(column).astype(np.int32)
else: # return floating point values
return line, column
PATH_TEST_DATA = os.environ.get('TEST_DATA', os.path.expanduser("~/Data/himawari/sample_bz2/HS_H08_20130710_0300_B16"))
......
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