Почему прямая проекция pyproj.Proj, по-видимому, не объясняет происхождение лат? - PullRequest
1 голос
/ 13 июня 2019

Меня смущает, как проекции, использующие 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).

Может кто-нибудь объяснить, почему это так?Мои знания о системах проекции / координат ограничены, но я приложил все усилия, чтобы понять Картографическая справка PROJ и связанную страницу викибук .

Большое спасибо заранее любомукто может помочь мне разобраться и направить меня в правильном направлении: -)

РЕДАКТИРОВАТЬ: Обновление 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), как и раньше, поэтому вопрос в том, какая другая информация необходима для строки проекции?

Ответы [ 2 ]

1 голос
/ 16 июня 2019

Я могу независимо подтвердить, что результат должен быть (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, в чем я сомневаюсь.

0 голосов
/ 20 июня 2019

Оказывается, я был глупым!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
Добро пожаловать на сайт PullRequest, где вы можете задавать вопросы и получать ответы от других членов сообщества.
...