def project(self, lon, lat):
"""Project lon, lat into grid coordinates."""
if self.proj_name == 'LAMBERT':
ellipsoids = {
'WGS-84': 'WGS84',
'GRS-80': 'GRS80',
'WGS-72': 'WGS72',
'Australian': 'aust_SA',
'Krasovsky': 'krass',
'International': 'new_intl',
'Hayford-1909': 'intl',
'Clarke-1880': 'clrk80',
'Clarke-1866': 'clrk66',
'Airy': 'airy',
'Bessel': 'bessel',
# 'Hayford-1830':
'Sphere': 'sphere'
}
try:
ellps = ellipsoids[self.proj_ellipsoid]
except KeyError:
raise ValueError(
'Ellipsoid not supported: {}'.format(self.proj_ellipsoid))
p = Proj(proj='lcc', lat_0=self.orig_lat, lon_0=self.orig_lon,
lat_1=self.first_std_paral, lat_2=self.second_std_paral,
ellps=ellps)
elif self.proj_name == 'SIMPLE':
p = Proj(proj='eqc', lat_0=self.orig_lat, lon_0=self.orig_lon)
else:
raise ValueError('Projection not supported: {}'.format(
self.proj_name))
x, y = p(lon, lat)
x = np.array(x)
y = np.array(y)
x[np.isnan(lon)] = np.nan
y[np.isnan(lat)] = np.nan
x /= 1000.
y /= 1000.
# inverse rotation
theta = np.radians(self.map_rot)
x1 = x*np.cos(theta) + y*np.sin(theta)
y1 = -x*np.sin(theta) + y*np.cos(theta)
x = x1
y = y1
# Try to return the same type of lon, lat
if not isinstance(lon, np.ndarray):
if isinstance(lon, Iterable):
x = type(lon)(x)
else:
x = float(x)
if not isinstance(lat, np.ndarray):
if isinstance(lat, Iterable):
y = type(lat)(y)
else:
y = float(y)
return x, y
评论列表
文章目录