# Why do pyproj and matplotlib basemap produce different results when they convert my lat-long data into a CRS?

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?

The correct solution is to usem.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)