r - 如何将经纬度转换为 NASA 的缓动网格?
问题描述
我正在处理从 nasa eosdis 下载的 he5 文件。
我使用 rhdf5 包成功读取了文件。
hdf 文件的子数据集由矩阵 dim 721 x 721 组成。
据我了解,该文件是作为easegrid 构建的,它没有关于坐标的任何信息,所以我无法找到我的工作场所在哪个网格(矩阵中的元素)中。
有什么方法可以转换纬度经度值以缓解网格?
太感谢了
解决方案
我发现了一个python软件,它能够在杨百翰大学的 FTP 服务器上的不同制图系统之间进行多次转换。
您可能希望ease2grid
在此文件中使用函数。
我还修复了文件上的几个 linting 问题。这是文件;
"""
Created on Sep 21, 2011
Revised on Apr 7, 2017 + include EASE2 support
@author: Bradley, DGL
"""
from __future__ import division
from numpy import cos, sin, tan, mod, sqrt, all
import numpy as np
from .ease2helper import ease2_map_info, easeconv_normalize_degrees
def latlon2pix(alon, alat, head):
"""Latitude/longitude to pixels
(x, y) = latlon2pix(lon,lat,head)
Convert a lat,lon coordinate (lon,lat) to an image pixel location
(x,y) (in floating point, matlab convention).
To compute integer pixel indices (ix,iy): check to insure
1 <= x < nsx+1 and 1 <= x < nsx+1 then ix=floor(x) iy=floor(y)
INPUTS:
lon,lat - longitude, latitude
head - header array from load sir
OUTPUTS:
x,y - pixel location (matlab coordinates y_matlab=nxy-y_sir+1)
"""
nsx = head[0] # noqa: F841
nsy = head[1]
iopt = head[16]
xdeg = head[2]
ydeg = head[3]
ascale = head[5]
bscale = head[6]
a0 = head[7]
b0 = head[8]
if iopt == -1: # image only (can't transform!)
x = ascale * (alon - a0)
y = bscale * (alat - b0)
elif iopt == 0: # rectalinear lat/lon
thelon = alon
thelat = alat
x = ascale * (thelon - a0)
y = bscale * (thelat - b0)
elif (iopt == 1) or (iopt == 2): # lambert
thelon, thelat = lambert1(alat, alon, ydeg, xdeg, iopt)
x = ascale * (thelon - a0)
y = bscale * (thelat - b0)
elif iopt == 5: # polar stereographic
thelon, thelat = polster(alon, alat, xdeg, ydeg)
x = (thelon - a0) / ascale
y = (thelat - b0) / bscale
elif (iopt == 8) or (iopt == 9) or (iopt == 10): # EASE2
thelon, thelat = ease2grid(iopt, alat, alon, ascale, bscale)
x = thelon + 1.0 - a0
y = thelat + 1.0 + b0
elif (iopt == 11) or (iopt == 12) or (iopt == 13): # EASE
thelon, thelat = easegrid(iopt, alat, alon, ascale)
thelon = thelon + xdeg
thelat = thelat + ydeg
x = thelon - (xdeg + a0)
y = thelat - (ydeg + b0)
else:
print("*** Unknown SIR transformation: %d" % iopt)
y = nsy - y - 1.0 # convert from matlab coordinates to SIR coordinates
return x, y
def lambert1(lat, lon, orglat, orglon, iopt):
"""Lambert azimuthal equal-area projection
(x,y)=lambert1(lat,lon,orglat,orglon,iopt)
Computes the transformation from lat/lon to x/y for the
lambert azimuthal equal-area projection
inputs:
lat (r): latitude +90 to -90 deg with north positive
lon (r): longitude 0 to +360 deg with east positive
or -180 to +180 with east more positive
orglat (r): origin parallel +90 to -90 deg with north positive
orglon (r): central meridian (longitude) 0 to +360 deg
or -180 to +180 with east more positive
iopt (i): earth radius option
for iopt=1 a fixed, nominal earth radius is used.
for iopt=2 the local radius of the earth is used.
outputs:
x,y (r): rectangular coordinates in km
see "map projections used by the u.s. geological survey"
geological survey bulletin 1532, pgs 157-173
for this routine, a spherical earth is assumed for the projection
the 1972 wgs ellipsoid model (bulletin pg 15).
the error will be small for small-scale maps.
"""
radearth = 6378.135 # equitorial earth radius
f = 298.26 # 1/f wgs 72 model values
dtr = 3.141592654 / 180.0
lon1 = mod(lon + 720.0, 360.0)
orglon1 = mod(orglon + 720.0, 360.0)
#
# compute local radius of the earth at center of image
#
eradearth = 6378.0 # use fixed nominal value
if iopt == 2: # local radius
era = 1.0 - 1.0 / f
eradearth = (
radearth
* era
/ sqrt(era * era * cos(orglat * dtr) ** 2 + sin(orglat * dtr) ** 2)
)
denom = (
1.0
+ sin(orglat * dtr) * sin(lat * dtr)
+ cos(orglat * dtr) * cos(lat * dtr) * cos(dtr * (lon1 - orglon1))
)
if all(denom > 0.0):
ak = sqrt(2.0 / denom)
else:
print("*** division error in lambert1 routine ***")
ak = 1.0
x = ak * cos(lat * dtr) * sin(dtr * (lon1 - orglon1))
y = ak * (
cos(dtr * orglat) * sin(dtr * lat)
- sin(dtr * orglat) * cos(dtr * lat) * cos(dtr * (lon1 - orglon1))
)
x = x * eradearth
y = y * eradearth
return x, y
def polster(alon, alat, xlam, slat):
"""Polar stereographic trasnformation
(x,y)=polster(lon,lat,xlam,slat)
computes the polar sterographic transformation for a lon,lat
input of (alon,alat) with reference origin lon,lat=(xlam,slat).
output is (x,y) in km
algorithm is the same as used for processing ers-1 sar images
as received from m. drinkwater (1994)
"""
# ported from polster.m by JPB 21 Sept 2011
e2 = 0.006693883
re = 6378.273
dtr = 3.141592654 / 180.0
e = sqrt(e2)
if slat < 0:
sn = -1.0
rlat = -alat
else:
sn = 1.0
rlat = alat
t = ((1.0 - e * sin(rlat * dtr)) / (1.0 + e * sin(rlat * dtr))) ** (e * 0.5)
ty = tan(dtr * (45.0 - 0.5 * rlat)) / t
if slat < 0:
rlat = -slat
else:
rlat = slat
t = ((1.0 - e * sin(dtr * rlat)) / (1.0 + e * sin(dtr * rlat))) ** (e * 0.5)
tx = tan(dtr * (45.0 - 0.5 * rlat)) / t
cm = cos(dtr * rlat) / sqrt(1.0 - e2 * sin(dtr * rlat) ** 2)
rho = re * cm * ty / tx
x = (sn * sin(dtr * (sn * alon - xlam))) * rho
y = -(sn * cos(dtr * (sn * alon - xlam))) * rho
return x, y
def easegrid(iopt, alat, alon, ascale):
"""EASE grid transformation
(thelon thelat)=easegrid(iopt,lat,lon,ascale)
computes the forward "ease" grid transform
given a lat,lon (alat,alon) and the scale (ascale) the image
transformation coordinates (thelon,thelat) are comuted
using the "ease grid" (version 1.0) transformation given in fortran
source code supplied by nsidc.
the radius of the earth used in this projection is imbedded into
ascale while the pixel dimension in km is imbedded in bscale
the base values are: radius earth= 6371.228 km
pixel dimen =25.067525 km
then, bscale = base_pixel_dimen
ascale = radius_earth/base_pixel_dimen
iopt is ease type: iopt=11=north, iopt=12=south, iopt=13=cylindrical
"""
# ported from easegrid.m by JPB 21 Sept 2011
pi2 = np.pi / 2.0
dtr = pi2 / 90.0
if iopt == 11: # ease grid north
thelon = ascale * sin(alon * dtr) * sin(dtr * (45.0 - 0.5 * alat))
thelat = ascale * cos(alon * dtr) * sin(dtr * (45.0 - 0.5 * alat))
elif iopt == 12: # ease grid south
thelon = ascale * sin(alon * dtr) * cos(dtr * (45.0 - 0.5 * alat))
thelat = ascale * cos(alon * dtr) * cos(dtr * (45.0 - 0.5 * alat))
elif iopt == 13: # ease cylindrical
thelon = ascale * pi2 * alon * cos(30.0 * dtr) / 90.0
thelat = ascale * sin(alat * dtr) / cos(30.0 * dtr)
return thelon, thelat
def ease2grid(iopt, alat, alon, ascale, bscale):
"""EASE2 grid transformation
(thelon thelat)=ease2grid(iopt,lat,lon,ascale,bscale)
given a lat,lon (alat,alon) and the scale (ascale) the image
transformation coordinates (thelon,thelat) are comuted
using the "ease2 grid" (version 2.0) transformation given in IDL
source code supplied by MJ Brodzik
RADIUS EARTH=6378.137 KM (WGS 84)
MAP ECCENTRICITY=0.081819190843 (WGS84)
inputs:
iopt: projection type 8=EASE2 N, 9-EASE2 S, 10=EASE2 T/M
alon, alat: lon, lat (deg) to convert (can be outside of image)
ascale and bscale should be integer valued)
ascale: grid scale factor (0..5) pixel size is (bscale/2^ascale)
bscale: base grid scale index (ind=int(bscale))
see ease2helper.py for definitions of isc and ind
outputs:
thelon: X coordinate in pixels (can be outside of image)
thelat: Y coordinate in pixels (can be outside of image)
"""
DTR = 0.01745329241994
ind = round(bscale)
isc = round(ascale)
dlon = alon
phi = DTR * alat
lam = dlon
# get base EASE2 map projection parameters
(
map_equatorial_radius_m,
map_eccentricity,
e2,
map_reference_latitude,
map_reference_longitude,
map_second_reference_latitude,
sin_phi1,
cos_phi1,
kz,
map_scale,
bcols,
brows,
r0,
s0,
epsilon,
) = ease2_map_info(iopt, isc, ind)
dlon = dlon - map_reference_longitude
dlon = easeconv_normalize_degrees(dlon)
lam = DTR * dlon
sin_phi = np.sin(phi)
q = (1.0 - e2) * (
(sin_phi / (1.0 - e2 * sin_phi * sin_phi))
- (1.0 / (2.0 * map_eccentricity))
* np.log(
(1.0 - map_eccentricity * sin_phi) / (1.0 + map_eccentricity * sin_phi)
)
)
if iopt == 8: # EASE2 grid north
qp = 1.0 - (
(1.0 - e2)
/ (2.0 * map_eccentricity)
* np.log((1.0 - map_eccentricity) / (1.0 + map_eccentricity))
)
rho = map_equatorial_radius_m * np.sqrt(qp - q)
if np.size(rho) > 1:
rho[np.abs(qp - q) < epsilon] = 0.0
else:
if np.abs(qp - q) < epsilon:
rho = 0
x = rho * np.sin(lam)
y = -rho * np.cos(lam)
elif iopt == 9: # EASE2 grid south
qp = 1.0 - (
(1.0 - e2)
/ (2.0 * map_eccentricity)
* np.log((1.0 - map_eccentricity) / (1.0 + map_eccentricity))
)
rho = map_equatorial_radius_m * np.sqrt(qp + q)
if np.size(rho) > 1:
rho[np.abs(qp - q) < epsilon] = 0.0
else:
if np.abs(qp - q) < epsilon:
rho = 0
x = rho * np.sin(lam)
y = rho * np.cos(lam)
elif iopt == 10: # EASE2 cylindrical
x = map_equatorial_radius_m * kz * lam
y = (map_equatorial_radius_m * q) / (2.0 * kz)
else:
print("*** invalid EASE2 projection specificaion in ease2grid")
thelon = r0 + (x / map_scale) + 0.5
thelat = s0 + (y / map_scale) + 0.5
return (thelon, thelat)
这段 Python 代码依赖于另一个名为 ease2helper 的 Python 代码。我在同一个 FTP 服务器上发现了这一点。
我还修复了这个文件的几个 linting 问题。这是ease2helper代码;
#!/usr/bin/env python
""" EASE2 grid helper utility functions for sirpy"""
######
# Imports
######
from __future__ import division
import numpy as np
##############################
# The EASE2 helper functions
##############################
def easeconv_normalize_degrees(dlon):
#
# Return dlon to within the range -180 <= dlon <= 180
# can handle array inputs
#
out = dlon
if np.size(out) > 1:
while (out < -180.0).sum() > 0:
out[out < -180.0] = out[out > 180.0] + 360.0
while (out > 180.0).sum() > 0:
out[out > 180.0] = out[out > 180.0] - 360.0
else:
while out < -180.0:
out = out + 360.0
while out > 180.0:
out = out - 360.0
return out
def ease2_map_info(iopt, isc, ind):
""" internally used routine
(map_equatorial_radius_m, map_eccentricity, \
e2, map_reference_latitude, map_reference_longitude, \
map_second_reference_latitude, sin_phi1, cos_phi1, kz, \
map_scale, bcols, brows, r0, s0, epsilon) = ease2_map_info(iopt, isc, nd)
defines EASE2 grid information
inputs:
iopt: projection type 8=EASE2 N, 9=EASE2 S, 10=EASE2 T/M
isc: scale factor 0..5 grid size is (basesize(ind))/2^isc
ind: base grid size index (map units per cell in m
NSIDC .grd file for isc=0
project type ind=0 ind=1 ind=2 ind=3
N EASE2_N25km EASE2_N30km EASE2_N36km EASE2_N24km
S EASE2_S25km EASE2_S30km EASE2_S36km EASE2_S24km
T/M EASE2_T25km EASE2_M25km EASE2_M36km EASE2_M24km
cell size (m) for isc=0 (scale is reduced by 2^isc)
project type ind=0 ind=1 ind=2 ind=3
N 25000.0 30000.0 36000.0 24000.0
S 25000.0 30000.0 36000.0 24000.0
T/M T25025.26 M25025.26000 M36032.220840584 M24021.480560389347
for a given base cell size (e.g., ind=0) isc is related to
NSIDC CETB .grd file names according to
isc N .grd name S .grd name T .grd name
0 EASE2_N25km EASE2_S25km EASE2_T25km
1 EASE2_N12.5km EASE2_S12.5km EASE2_T12.5km
2 EASE2_N6.25km EASE2_S6.25km EASE2_T6.25km
3 EASE2_N3.125km EASE2_S3.125km EASE2_T3.125km
4 EASE2_N1.5625km EASE2_S1.5625km EASE2_T1.5625km
outputs
map_equatorial_radius_m EASE2 Earth equitorial radius (km) [WGS84]
map_eccentricity EASE2 Earth eccentricity [WGS84]
map_reference_latitude Reference latitude (deg)
map_reference_longitude Reference longitude (deg)
map_second_reference_latitude Secondary reference longitude* (deg)
sin_phi1, cos_phi1 kz EASE2 Cylin parameters*
map_scale EASE2 map projection pixel size (km)
bcols, brows, EASE2 grid size in pixels
r0, s0 EASE2 base projection size in pixels
epsilon EASE2 near-polar test factor
"""
DTR = 0.01745329241994
m = 2 ** np.floor(isc) # compute power-law scale factor
map_equatorial_radius_m = 6378137.0 # WGS84
map_eccentricity = 0.081819190843 # WGS84
e2 = map_eccentricity * map_eccentricity
map_reference_longitude = 0.0
epsilon = 1.0e-6
# map-specific parameters
if iopt == 8: # EASE2 grid north
map_reference_latitude = 90.0
if ind == 1: # EASE2_N30km.gpd
base = 30000.0
nx = 600
ny = 600
elif ind == 2: # EASE2_N36km.gpd
base = 36000.0
nx = 500
ny = 500
elif ind == 3: # EASE2_N24km.gpd
base = 24000.0
nx = 750
ny = 750
else: # EASE2_N25km.gpd
base = 25000.0
nx = 720
ny = 720
map_second_reference_latitude = 0.0
sin_phi1 = 0.0
cos_phi1 = 1.0
kz = cos_phi1
elif iopt == 9: # EASE2 grid south
map_reference_latitude = -90.0
if ind == 1: # EASE2_S30km.gpd
base = 30000.0
nx = 600
ny = 600
elif ind == 2: # EASE2_S36km.gpd
base = 36000.0
nx = 500
ny = 500
elif ind == 3: # EASE2_S24km.gpd
base = 24000.0
nx = 750
ny = 750
else: # EASE2_S25km.gpd
base = 25000.0
nx = 720
ny = 720
map_second_reference_latitude = 0.0
sin_phi1 = 0.0
cos_phi1 = 1.0
kz = cos_phi1
elif iopt == 10: # EASE2 cylindrical
map_reference_latitude = 0.0
map_second_reference_latitude = 30.0
sin_phi1 = np.sin(DTR * map_second_reference_latitude)
cos_phi1 = np.cos(DTR * map_second_reference_latitude)
kz = cos_phi1 / np.sqrt(1.0 - e2 * sin_phi1 * sin_phi1)
if ind == 1: # EASE2_M25km.gpd
base = 25025.26000
nx = 1388
ny = 584
elif ind == 2: # EASE2_M36km.gpd
base = 36032.220840584
nx = 964
ny = 406
elif ind == 3: # EASE2_M24km.gpd
base = 24021.480560389347
nx = 1446
ny = 609
else: # EASE2_T25km.gpd
base = 25025.26000
nx = 1388
ny = 540
else:
print("*** invalid EASE2 projection code ***")
# grid info
if isc >= 0:
map_scale = base / m
bcols = np.ceil(nx * m)
brows = np.ceil(ny * m)
r0 = (nx * m - 1) / 2
s0 = (ny * m - 1) / 2
else:
map_scale = base * m
bcols = np.ceil(nx / np.float(m))
brows = np.ceil(ny / np.float(m))
r0 = (nx / np.float(m) - 1) / 2
s0 = (ny / np.float(m) - 1) / 2
return (
map_equatorial_radius_m,
map_eccentricity,
e2,
map_reference_latitude,
map_reference_longitude,
map_second_reference_latitude,
sin_phi1,
cos_phi1,
kz,
map_scale,
bcols,
brows,
r0,
s0,
epsilon,
)
推荐阅读
- r - 我可以获取 install.packages 将使用的 URL 吗?
- go - 代理处理函数中 []byte 的字符串表示
- swift - 如何将数据作为时间戳数据类型发布到 Firestore - swift
- python - Django V2.1 reverse 不是有效的视图函数或模式名称
- unit-testing - 无法在厨师食谱中运行 rspec 单元测试“touch_file”
- ios - Swift iOS - 选择新单元格时更改旧单元格的背景颜色
- android - React Native 未在发行版中编译
- android - NoClassDefFoundError:解析失败:Lcom/google/android/gms/common/internal/zzx
- php - WooCommerce 迷你购物车事件触发器
- python - 学习 Python 和 Scrapy