为什么 pyproj.Proj 前向投影似乎不能解释经纬度原点?
Why does pyproj.Proj forward projection not seem to account for lat lon origin?
我对如何根据切点/经纬度原点定义使用 pyproj.Proj 的投影感到困惑。
考虑以下代码:
import pyproj
p = pyproj.Proj('+proj=tmerc +lat_0=55 +lon_0=-1 +a=6378137 +b=6356752.3 +units=m +no_defs')
x, y = p(55, -1)
既然我指定了纬度和经度的原点,我希望在指定这些坐标时我能够 assert x == 0 and y == 0
,但实际上我得到 (7571700.820174289, -6296411.725576388)
.
谁能解释一下为什么会这样?我对 projection/coordinate 系统的了解有限,但我已尽力了解 PROJ Cartographic help and a related wikibook page。
非常感谢任何能帮助我走上正确方向的人:-)
编辑:更新 1
感谢@lusitanica 和他们的有用回答,我现在尝试将比例因子设置为 1 并重新运行:
x, y = pyproj.Proj('+proj=tmerc +lat_0=55 +lon_0=-1 +k_0=1 +a=6378137 +b=6356752.3 +units=m +no_defs', preserve_units=True)(55, -1)
不幸的是,这和以前一样给出了 (7571700.820174289, -6296411.725576388)
,所以问题是投影字符串还需要哪些其他信息?
我可以独立确认结果应该是(0,0)
from math import *
# GRS-80
a = 6378137.
equad =0.00669437999
# Natural Origin
lat0=55.
lon0=-1.
############################################################################
# Meridian Arc
############################################################################
def arcmer(a,equad,lat1,lat2):
b=a*sqrt(1-equad)
n=(a-b)/(a+b)
a0=1.+((n**2)/4.)+((n**4)/64.)
a2=(3./2.)*(n-((n**3)/8.))
a4=(15./16.)*((n**2)-((n**4)/4.))
a6=(35./48.)*(n**3)
s1=a/(1+n)*(a0*lat1-a2*sin(2.*lat1)+a4*sin(4.*lat1)-a6*sin(6.*lat1))
s2=a/(1+n)*(a0*lat2-a2*sin(2.*lat2)+a4*sin(4.*lat2)-a6*sin(6.*lat2))
return s2-s1
#############################################################################
# Direct projection Gauss-Kruger
#############################################################################
def geogauss(lat,lon,a,equad,lat0,lon0):
lat0=radians(lat0)
lon0=radians(lon0)
lat=radians(lat)
lon=radians(lon)
lon=lon-lon0
N=a/sqrt(1-equad*(sin(lat))**2)
RO=a*(1-equad)/((1-equad*(sin(lat)**2))**(3./2.))
k1=(N/RO)+(4.*(N**2)/(RO**2))-((tan(lat))**2)
k2=(N/RO)-((tan(lat))**2)
k3=N/RO*(14.-58.*((tan(lat)))**2)+40.*((tan(lat))**2)+((tan(lat))**4)-9.
x=lon*N*cos(lat)+(lon**3)/6.*N*((cos(lat))**3)*k2+(lon**5)/120.*N*((cos(lat))**5)*k3
y=arcmer(a,equad,lat0,lat)+(lon**2)/2.*N*sin(lat)*cos(lat)+((lon**4)/24.)*N*sin(lat)*((cos(lat))**3)*k1
return x,y
lat=55.
lon=-1.
coordinates = geogauss(lat,lon,a,equad,lat0,lon0)
print lat,lon
print "x= %.3f" %coordinates[0]
print "y= %.3f" %coordinates[1]
结果:
55.0 -1.0
x= 0.000
y= 0.000
在我的代码中,我正在考虑 1
的自然起源的比例因子。
在您的 PROJ 字符串中我看不到任何内容。尝试将比例因子设置为 +k_0=1
如果这没有帮助,那么您的 PROJ 字符串中还缺少其他内容,因为结果确实 (0,0)
.
除非 PyProj 存在错误,我对此表示怀疑。
原来是我傻了! pyproj.Proj
期望输入 lon,lat 顺序不是 lat,lon :-(
因此以下工作:
import pyproj
lat0, lon0 = 55, -1
p = pyproj.Proj(f'+proj=tmerc +lat_0={lat0} +lon_0={lon0} +a=6378137 +b=6356752.3 +units=m +no_defs')
x, y = p(lon0, lat0)
print(x, y) # 0.0, 0.0
我对如何根据切点/经纬度原点定义使用 pyproj.Proj 的投影感到困惑。
考虑以下代码:
import pyproj
p = pyproj.Proj('+proj=tmerc +lat_0=55 +lon_0=-1 +a=6378137 +b=6356752.3 +units=m +no_defs')
x, y = p(55, -1)
既然我指定了纬度和经度的原点,我希望在指定这些坐标时我能够 assert x == 0 and y == 0
,但实际上我得到 (7571700.820174289, -6296411.725576388)
.
谁能解释一下为什么会这样?我对 projection/coordinate 系统的了解有限,但我已尽力了解 PROJ Cartographic help and a related wikibook page。
非常感谢任何能帮助我走上正确方向的人:-)
编辑:更新 1
感谢@lusitanica 和他们的有用回答,我现在尝试将比例因子设置为 1 并重新运行:
x, y = pyproj.Proj('+proj=tmerc +lat_0=55 +lon_0=-1 +k_0=1 +a=6378137 +b=6356752.3 +units=m +no_defs', preserve_units=True)(55, -1)
不幸的是,这和以前一样给出了 (7571700.820174289, -6296411.725576388)
,所以问题是投影字符串还需要哪些其他信息?
我可以独立确认结果应该是(0,0)
from math import *
# GRS-80
a = 6378137.
equad =0.00669437999
# Natural Origin
lat0=55.
lon0=-1.
############################################################################
# Meridian Arc
############################################################################
def arcmer(a,equad,lat1,lat2):
b=a*sqrt(1-equad)
n=(a-b)/(a+b)
a0=1.+((n**2)/4.)+((n**4)/64.)
a2=(3./2.)*(n-((n**3)/8.))
a4=(15./16.)*((n**2)-((n**4)/4.))
a6=(35./48.)*(n**3)
s1=a/(1+n)*(a0*lat1-a2*sin(2.*lat1)+a4*sin(4.*lat1)-a6*sin(6.*lat1))
s2=a/(1+n)*(a0*lat2-a2*sin(2.*lat2)+a4*sin(4.*lat2)-a6*sin(6.*lat2))
return s2-s1
#############################################################################
# Direct projection Gauss-Kruger
#############################################################################
def geogauss(lat,lon,a,equad,lat0,lon0):
lat0=radians(lat0)
lon0=radians(lon0)
lat=radians(lat)
lon=radians(lon)
lon=lon-lon0
N=a/sqrt(1-equad*(sin(lat))**2)
RO=a*(1-equad)/((1-equad*(sin(lat)**2))**(3./2.))
k1=(N/RO)+(4.*(N**2)/(RO**2))-((tan(lat))**2)
k2=(N/RO)-((tan(lat))**2)
k3=N/RO*(14.-58.*((tan(lat)))**2)+40.*((tan(lat))**2)+((tan(lat))**4)-9.
x=lon*N*cos(lat)+(lon**3)/6.*N*((cos(lat))**3)*k2+(lon**5)/120.*N*((cos(lat))**5)*k3
y=arcmer(a,equad,lat0,lat)+(lon**2)/2.*N*sin(lat)*cos(lat)+((lon**4)/24.)*N*sin(lat)*((cos(lat))**3)*k1
return x,y
lat=55.
lon=-1.
coordinates = geogauss(lat,lon,a,equad,lat0,lon0)
print lat,lon
print "x= %.3f" %coordinates[0]
print "y= %.3f" %coordinates[1]
结果:
55.0 -1.0
x= 0.000
y= 0.000
在我的代码中,我正在考虑 1
的自然起源的比例因子。
在您的 PROJ 字符串中我看不到任何内容。尝试将比例因子设置为 +k_0=1
如果这没有帮助,那么您的 PROJ 字符串中还缺少其他内容,因为结果确实 (0,0)
.
除非 PyProj 存在错误,我对此表示怀疑。
原来是我傻了! pyproj.Proj
期望输入 lon,lat 顺序不是 lat,lon :-(
因此以下工作:
import pyproj
lat0, lon0 = 55, -1
p = pyproj.Proj(f'+proj=tmerc +lat_0={lat0} +lon_0={lon0} +a=6378137 +b=6356752.3 +units=m +no_defs')
x, y = p(lon0, lat0)
print(x, y) # 0.0, 0.0