Why does pyproj.Proj forward projection not seem to account for lat lon origin?

Question

I am getting confused about how projections using pyproj.Proj are defined with respect to a point of tangency / lat lon origin.

Consider the following code:

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)

Now given that I specified the origin for latitude and longitude, I would expect that when specifying those coordinates I would be able to assert x == 0 and y == 0, however I actually get (7571700.820174289, -6296411.725576388).

Can anyone explain why this is the case? My knowledge of projection/coordinate systems is limited, but I did my best to understand PROJ Cartographic help and a related wikibook page.

Many thanks in advance to any who can help set me straight and put me in the right direction :-)

EDIT: Update 1

Thanks to @lusitanica and their helpful answer, I have now tried setting the scale factor to 1 and rerunning:

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)

Unfortunately this gives (7571700.820174289, -6296411.725576388) as before, so the question is what other information is needed for the projection string?

score 1 · Answer 1 · edited Jun 20 '20 at 09:12

I can confirm independently that the result should be (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]

Result:

55.0 -1.0

x= 0.000

y= 0.000

In my code I'm considering a Scale Factor at Natural Origin of 1.

In your PROJ string I can't see any. Try to set the Scale factor to +k_0=1

If that doesn't help there's something else missing from your PROJ string because the result is realy (0,0).

Unless there is a bug with PyProj, which I doubt.

Thanks @lusitanica for the explanation and detailed alternative method! Sadly the scaling factor doesn't seem to do anything, so I'm thinking that there is some missing information needed for the projection string. I too doubt that there is a bug with PyProj itself (I'm using version 2.2.0), however, I won't rule this out (maybe with the underlying PROJ.4 library as `cs2cs` is a circa 2016 vintage and so should probably be updated). — Liam Deacon, Jun 17 '19 at 04:30
Try to isolate the problem. For example setting `+ellps=GRS80` instead of `a` and `b`. Also try with an empty string `''` to see if the parameters are being used. Tweak with all the possibilities in https://pyproj4.github.io/pyproj/stable/api/proj.html#pyproj-proj — LuisTavares, Jun 17 '19 at 08:51

score 0 · Accepted Answer · answered Jun 20 '19 at 05:08

Turns out I was being silly! pyproj.Proj expects input in lon, lat order not lat, lon :-(

Thus the following works:

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

Why does pyproj.Proj forward projection not seem to account for lat lon origin?

EDIT: Update 1

2 Answers2