Sola/proj_info
1
0
Fork 0
mirror of https://github.com/SolaProject/proj_info.git synced 2025-10-20 15:27:25 +08:00
Code Issues Packages Projects Releases Wiki Activity

Create proj_info.py

Browse Source
This commit is contained in:
SolaProject
2023-03-19 14:55:37 +08:00
committed by GitHub
parent f7ac26625a
commit 43f90ca2ae
1 changed files with 275 additions and 0 deletions
Download Patch File Download Diff File Expand all files Collapse all files

275
proj_info.py Normal file
View File

@ -0,0 +1,275 @@
import numpy as np
from math import radians, cos, tan, log10, sin, sqrt, atan2, atan, degrees
"""
更新记录:
2022-09-23 11:59:01 Sola v1 编写源代码, 修正set_lc代码错误的问题
"""
EARTH_RADIUS_M = 6370000.
class proj_info(object):
"""参考WPS的geogrid源代码写的"""
def __init__(self, code, lat1=None, lon1=None, lat0=None, lon0=None,
dx=None, dy=None, latinc=None, loninc=None, stdlon=None,
truelat1=None, truelat2=None, phi=None, lambda1=None,
ixdim=None, jydim=None, stagger='HH', nlat=None, nlon=None, nxmin=None,
nxmax=None, hemi=None, cone=None, polei=None, polej=None,
rsw=None, knowni=None, knownj=None, re_m=EARTH_RADIUS_M,
init=False, wrap=False, rho0=None, nc=None, bigc=None, comp_ll=False,
gauss_lat=None) -> None:
self.code = code
self.lat1 = lat1 # SW latitude (1,1) in degrees (-90->90N) 格点(1, 1)纬度, 西南角, 度
self.lon1 = lon1 # SW longitude (1,1) in degrees (-180->180E) 格点(1, 1)经度, 西南角, 度
self.lat0 = lat0
self.lon0 = lon0
self.dx = dx # Grid spacing in meters at truelats, used x方向网格距, m
self.dy = dy # Grid spacing in meters at truelats, used y方向网格距, m
self.latinc = latinc
self.loninc = loninc
self.stdlon = stdlon # Longitude parallel to y-axis (-180->180E) 中央经线, 与网格上所有y方向的线平行
self.truelat1 = truelat1 # First true latitude (all projections) 标准纬线1, 所有投影都需要
self.truelat2 = truelat2 # Second true lat (LC only) 标准纬线2, 仅兰伯特投影需要
self.phi = phi
self.lambda1 = lambda1
self.ixdim = ixdim
self.jydim = jydim
self.stagger = stagger
self.nlat = nlat
self.nlon = nlon
self.nxmin = nxmin
self.nxmax = nxmax
self.hemi = hemi # 1 for NH, -1 for SH 判断位于哪个半球, 北半球则为1, 南半球则为-1
self.cone = cone # Cone factor for LC projections 兰伯特投影的角度因素
self.polei = polei # Computed i-location of pole point 极点的i坐标
self.polej = polej # Computed j-location of pole point 极点的j坐标
self.rsw = rsw # Computed radius to SW corner 西南角(左下角)半径
self.knowni = knowni # X-location of known lat/lon 已知的点位的i坐标
self.knownj = knownj # Y-location of known lat/lon 已知的点位的j坐标（一般就是ctl文件里面的那个点，也就是中心点）
self.re_m = re_m # Radius of spherical earth, meters 地球半径, 6370000, m
self.init = init
self.wrap = wrap
self.rho0 = rho0
self.nc = nc
self.bigc = bigc
self.comp_ll = comp_ll
self.gauss_lat = gauss_lat
class proj_LC(proj_info):
"""
参考WPS源码中的proj_LC改写, 因为WRF计算得到的网格与cartopy的不同
更新记录:
2022-09-22 22:07:51 Sola 编写源代码
"""
def __init__(self, code='PROJ_LC', truelat1=None, truelat2=None, lat1=None,
lon1=None, knowni=None, knownj=None, stdlon=None, dx=None,
dy=None, nx=None, ny=None, re_m=EARTH_RADIUS_M) -> None:
"""
初始化
必要参数:
code 投影编码
truelat1 标准纬线1
truelat2 标准纬线2
lat1 参考点纬度
lon1 参考点经度
stdlon 中央经线
dx x方向网格距(m)
nx x方向格点数
ny y方向格点数
可选参数:
knowni 参考点x方向坐标, 默认为网格中心
knownj 参考点y方向坐标, 默认为网格中心
dy y方向网格距(m), 默认与dx一致
re_m 地球半径, 默认为6370000
"""
if truelat1 is None or truelat2 is None or lat1 is None or lon1 is None\
or nx is None or ny is None or stdlon is None or dx is None:
print('[ERROR] cannot generate proj!')
if abs(lat1) > 90 or dx <= 0 or truelat1 > 90:
pass
dummy_lon1 = (lon1 + 180) % 360 - 180 # 限制经度范围
dummy_stdlon = (stdlon + 180) % 360 - 180 # 限制中央经线范围
if knowni is None and knownj is None:
knowni = (nx + 1) / 2
knownj = (ny + 1) / 2
if dy is None: # 设置dy, 如果dy不存在, 则利用dx给定
dy = dx
if truelat1 < 0: # 所在半球, 1为北半球, -1为南半球
hemi = -1
else:
hemi = 1
if abs(truelat2) > 90: # 如果标准纬线2超过范围, 则用标准纬线1赋值
truelat2 = truelat1
super().__init__(code=code, lat1=lat1, lon1=dummy_lon1, dx=dx, dy=dy,
stdlon=dummy_stdlon, truelat1=truelat1, truelat2=truelat2, hemi=hemi,
knowni=knowni, knownj=knownj, re_m=re_m) # 初始化各变量
self.rebydx = re_m / dx # 地球半径除以网格距
self.set_lc() # 计算其他变量
self.check_init() # 确认是否所有变量都计算完毕
def set_lc(self):
"""初始化兰伯特投影"""
# 这些参数都是从proj结构体中获得的, 目的是获得cone, 一个与圆锥角度相关的量, 应该是正弦值
self.lc_cone()
# 左下角网格的经度与中央经线的差
deltalon1 = self.lon1 - self.stdlon
# 将这个差值限制在 [-180, 180] 之间
deltalon1 = (deltalon1 + 180) % 360 - 180
# Convert truelat1 to radian and compute COS for later use, 将第一条标准纬线计算为余弦值, 接下来要用
ctl1r = cos(radians(self.truelat1))
# Compute the radius to our known lower-left (SW) corner 计算距离左下角的半径?
# 变量说明: rebydx: 地球半径/x方向网格距; hemi: 半球;
self.rsw = self.rebydx * ctl1r / self.cone * (tan(radians(90*self.hemi \
- self.lat1)/2)/tan(radians(90*self.hemi - self.truelat1)/2))**self.cone
# Find pole point 找到极点
arg = self.cone * radians(deltalon1)
self.polei = self.hemi * self.knowni - self.hemi * self.rsw * sin(arg)
self.polej = self.hemi * self.knownj + self.rsw * cos(arg)
def lc_cone(self):
"""计算切面圆锥的角度"""
# 首先, 看是割线投影还是切线投影
if abs(self.truelat1 - self.truelat2) > 0.1:
self.cone = log10(cos(radians(self.truelat1))) - log10(cos(radians(self.truelat2)))
self.cone = self.cone/(log10(tan(radians(90-abs(self.truelat1))/2))\
- log10(tan(radians(90-abs(self.truelat2))/2)))
else:
self.cone = sin(radians(abs(self.truelat1)))
def llij_lc(self, lon, lat):
"""通过经纬度计算坐标"""
if not self.init:
print('[ERROR] proj cannot use!')
deltalon = lon - self.stdlon # 计算经度与中央经线差
deltalon = (deltalon + 180) % 360 - 180 # 限定范围 -180~180
ctl1r = cos(radians(self.truelat1)) # 剩下就看不懂了
rm = self.rebydx * ctl1r / self.cone * (tan(radians(90*self.hemi \
- lat)/2)/tan(radians(90*self.hemi - self.truelat1)/2))**self.cone
arg = self.cone * radians(deltalon)
i = self.polei + self.hemi * rm * sin(arg)
j = self.polej - rm * cos(arg)
i = self.hemi * i
j = self.hemi * j
return i, j
def llijs_lc(self, lon, lat):
"""通过经纬度序列计算坐标"""
if not self.init:
print('[ERROR] proj cannot use!')
deltalon = lon - self.stdlon # 计算经度与中央经线差
deltalon = (deltalon + 180) % 360 - 180 # 限定范围 -180~180
ctl1r = np.cos(np.radians(self.truelat1)) # 剩下就看不懂了
rm = self.rebydx * ctl1r / self.cone * (np.tan(np.radians(90*self.hemi \
- lat)/2)/np.tan(np.radians(90*self.hemi - self.truelat1)/2))**self.cone
arg = self.cone * np.radians(deltalon)
i = self.polei + self.hemi * rm * np.sin(arg)
j = self.polej - rm * np.cos(arg)
i = self.hemi * i
j = self.hemi * j
return i, j
def ijll_lc(self, i, j):
"""通过坐标计算经纬度"""
if not self.init:
print('[ERROR] proj cannot use!')
chi1 = radians(90 - self.hemi * self.truelat1)
chi2 = radians(90 - self.hemi * self.truelat2)
inew = self.hemi * i
jnew = self.hemi * j
xx = inew - self.polei
yy = self.polej - jnew
r2 = xx**2 + yy**2
r = sqrt(r2) / self.rebydx
if r2 == 0.:
lat = self.hemi * 90
lon = self.stdlon
else:
lon = self.stdlon + degrees(atan2(self.hemi*xx, yy))/self.cone
lon = (lon + 180) % 360 - 180
if chi1 == chi2:
chi = 2 * atan((r/tan(chi1))**(1/self.cone) * tan(chi1*0.5))
else:
chi = 2 * atan((r*self.cone/sin(chi1))**(1/self.cone) * tan(chi1*0.5))
lat = (90 - degrees(chi)) * self.hemi
return lon, lat
def ijlls_lc(self, i, j):
"""
通过坐标计算经纬度
2022-09-28 18:22:26 Sola 修正计算chi时, xx, yy未进行筛选的问题
"""
if not self.init:
print('[ERROR] proj cannot use!')
chi1 = np.radians(90 - self.hemi * self.truelat1)
chi2 = np.radians(90 - self.hemi * self.truelat2)
inew = self.hemi * i
jnew = self.hemi * j
xx = inew - self.polei
yy = self.polej - jnew
r2 = xx**2 + yy**2
r = np.sqrt(r2) / self.rebydx
lon = np.zeros(r2.shape)
lat = np.zeros(r2.shape)
select_array = r2==0
lat[select_array] = self.hemi * 90
lon[select_array] = self.stdlon
lon[~select_array] = self.stdlon + np.degrees(np.arctan2(
self.hemi*xx[~select_array], yy[~select_array]))/self.cone
lon[~select_array] = (lon[~select_array] + 180) % 360 - 180
chi = np.zeros(r2.shape)
chi[chi1==chi2] = 2 * np.arctan((r/np.tan(chi1))**(1/self.cone) * np.tan(chi1*0.5))
chi[chi1!=chi2] = 2 * np.arctan((r*self.cone/np.sin(chi1))**(1/self.cone) * np.tan(chi1*0.5))
lat[~select_array] = (90 - np.degrees(chi[~select_array])) * self.hemi
return lon, lat
def check_init(self):
"""确认是否所有变量都已经准备好了"""
if self.cone is None or self.rsw is None or self.polei is None or \
self.polej is None:
print('[ERROR] cannot set proj_lc!')
else:
self.init = True
def grid_id_float(self, lon, lat):
"""返回以0为开始的下标"""
i,j = self.llij_lc(lon, lat)
return i - 1, j - 1
def grid_ids_float(self, lon_array, lat_array):
"""返回以0为开始的下标数组"""
i_array, j_array = self.llijs_lc(lon_array, lat_array)
return i_array - 1, j_array - 1
def grid_lonlat(self, ix, iy):
"""返回对应网格(以0为下标开始)的经纬度"""
lon, lat = self.ijll_lc(ix + 1, iy + 1)
return lon, lat
def grid_lonlats(self, ix_array, iy_array):
"""返回对应网格数组(以0为下标开始)的经纬度数组"""
lon_array, lat_array = self.ijlls_lc(ix_array + 1, iy_array + 1)
return lon_array, lat_array
def transform_point(self, lon, lat, proj_useless=None):
"""返回对应经纬度坐标的网格坐标(m)"""
ix, iy = self.llij_lc(lon, lat)
cx, cy1 = self.llij_lc(self.stdlon, self.truelat1)
cx, cy2 = self.llij_lc(self.stdlon, self.truelat2)
cy = (cy1 + cy2) / 2
return (ix - cx) * self.dx, (iy - cy) * self.dy
if __name__ == '__main__':
proj = proj_LC(truelat1=45, truelat2=15, lat1=30, lon1=108, stdlon=108, dx=3000, dy=3000, nx=2025, ny=2025)
proj.llij_lc(108, 30)
proj.ijll_lc(1013, 1013)
x = np.arange(1, 2000, 10)
y = np.arange(1, 2000, 10)
print(proj.ijlls_lc(x, y))
print(proj.ijll_lc(x[0], y[0]), proj.ijll_lc(x[-1], y[-1]))
x, y = proj.llij_lc(0, 90)
lon, lat = proj.ijll_lc(x, y)
x, y = proj.llij_lc(lon, lat)
lon1, lat1 = proj.ijll_lc(x, y)
print(lon, lat, x, y, lon1, lat1)