Basemap and proj get different projection result

Question

I am struggling to understand the projections using pyproj.

for now my question is to understand the results of projection operations

I have following coordinates that I project on x,y:

from mpl_toolkits.basemap import Basemap
import pyproj

lon  = [3.383789, 5.822754]
lat = [48.920575, 53.72185]

# with Basemap
M = Basemap(projection='merc',ellps = 'WGS84')
q1, q2 = M(lon, lat)
for a,b in zip([x for x in q1], [x for x in q2]):
    print(a,b)
# with pyproj
from pyproj import Proj
p = Proj(proj='merc', ellps='WGS84',errcheck = True)
p1 = Proj(proj='latlong', datum='WGS84',errcheck = True)
print(p(3.383789, 48.920575), p(5.822754, 53.72185))
print(p1(3.383789, 48.920575), p1(5.822754, 53.72185))

20414190.011221122 65799915.8523339 20685694.35308374 66653928.94763097 (376681.6684318804, 6229168.979819128) (648186.0102944968, 7083182.075116195) (0.0590582592427664, 0.8538251057188251) (0.10162622883366988, 0.9376231627625158)

why are the results different while I use the same projections parameters as a newbie in geospatial data processing I apologize in advance for a question which may be trivial

score 0 · Answer 1 · answered Nov 09 '18 at 13:36

For Mercator projection, Basemap uses the origin of the grid coordinate system at lower-left extent of the computable values. With your code, the values can be computed as

M(0, 0, inverse=True)
# output: (-180.0, -89.98999999999992)

If you compute projection coordinates for (long=0, lat=0) and assign the values to (x0, y0). You get the shifts in coordinates (x0, y0) that make its projection coordinates different from standard values (0,0 at center of the map).

lon0, lat0 = 0, 0
x0, y0 = M(lon0, lat0)
# x0=20015077.371242613, y0=59546805.8807

For a test point at (long=3.383789, lat=48.920575),

lon1 = 3.383789
lat1 = 48.920575
x1, y1 = M(lon1, lat1)

with the coordinate shifts applied, the result is

print(x1-x0, y1-y0)
# output: (376259.9924608879, 6254386.049398325)

when compare with the values from pyproj

p0 = Proj(proj='merc', ellps='WGS84', errcheck = True)
print(p0(lon1, lat1))
# output (376681.6684318804, 6229168.979819128)

they are quite agree but not close. For small scale map plotting, you can't see the discrepancies on the maps.

Basemap and proj get different projection result

1 Answers1