[GIS] Why do pyproj and matplotlib basemap produce different results when they convert the lat-long data into a CRS

coordinate systemmatplotlib-basemappyprojpython

pyproj and matplotlib basemap seem to produce differing x-y values when they convert my lat-long data into the same CRS (such as Albers Equal Area). These x-y values tend to be off by hundreds of meters. I must be doing something wrong but I cannot figure out what.

Here is my code (I also posted this question at StackOverflow). I take a lat/lon and convert it to the defined CRS (Albers) using basemap and then pyproj:

from mpl_toolkits.basemap import Basemap
from pyproj import Proj, transform

lat = 48.57446
lon = 9.066455

# specify the width and height of the map domain in projection coordinates (meters)
map_width_m = 4000 * 1000. #4000 km
map_height_m = 2500 * 1000. #2500 km

# define my projected CRS as albers equal area
albers_crs = {'proj':'aea', 
              'lat_1':35., 
              'lat_2':55., 
              'lon_0':10., 
              'lat_0':45., 
              'x_0':map_width_m/2., 
              'y_0':map_height_m/2.}

# first project the point using mpl basemap
m = Basemap(width=map_width_m, 
            height=map_height_m,
            projection=albers_crs['proj'],
            lat_1=albers_crs['lat_1'], 
            lat_2=albers_crs['lat_2'], 
            lon_0=albers_crs['lon_0'], 
            lat_0=albers_crs['lat_0'])
basemap_x, basemap_y = m(lon, lat)
print 'basemap:', basemap_x, basemap_y

# now project the point using pyproj
pyproj_convert = Proj(albers_crs)
pyproj_x, pyproj_y = pyproj_convert(lon, lat)
print 'pyproj:', pyproj_x, pyproj_y

# what's the difference between the two, in meters?
print 'difference:', basemap_x - pyproj_x, basemap_y - pyproj_y

Output:

basemap: 1932284.3542 1653858.27802
pyproj: 1932077.41654 1653737.11296
difference: 206.937659041 121.165060601

As you can see, my x's are 207 meters off each other and my y's are 121 meters off each other. The same holds for other projections, such as Lambert Conformal Conic. I have several hundred spatial data points in this dataset and these other lat/lon inputs can produce x-y values from basemap and pyproj that are 1000s of meters off each other. If I set the latitude and longitude of false origin in my CRS definition equal to my lat and lon variables (above), basemap and pyproj produce the same x and y values. But I need to set the latitude and longitude of false origin to the center of my map (to center it at the center of my entire spatial data set).

Why do basemap and pyproj produce different results? Why aren't basemap and pyproj producing the exact same x and y values given the exact same lat/lon and the exact same CRS definition?

Best Answer

The correct solution is to use m.proj4stringand not directly albers_crs because it is not a "valid" Proj4 string (even if it apparently works: major axis or radius = 0 not given)

print m.proj4string
'+lon_0=10.0 +y_0=1250000.0 +R=6370997.0 +proj=aea +x_0=2000000.0 +units=m +lat_2=55.0 +lat_1=35.0 +lat_0=45.0 '
# so
pyproj_convert = Proj(m.proj4string) 
pyproj_x, pyproj_y = pyproj_convert(lon, lat)
print 'pyproj:', pyproj_x, pyproj_y
pyproj: 1932284.3542 1653858.27802

You can control the result with the Python module GDAL

from osgeo import osr
from_srs = osr.SpatialReference()
from_srs.ImportFromEPSG(4326)
to_srs =  osr.SpatialReference()
to_srs.ImportFromProj4(m.proj4string)
transf = osr.CoordinateTransformation(from_srs,to_srs)
print transf.TransformPoint(lon,lat)[:2]
(1932284.3541970258, 1653858.278018097)
Related Question