0

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?

  • 1
    Make sure both projected CRS are using the base GCS, although I would expect that to only affect the Y values. And what's going on with x_0 and y_0? Are those getting picked up correctly in the Basemap? What happens if you leave them set to 0,0? – mkennedy Nov 5 '15 at 17:46
  • cross posted at stackoverflow.com/questions/33536770/… – Ian Turton Nov 5 '15 at 17:56
  • please see meta.gis.stackexchange.com/questions/4067/… to see a discussion of why we discourage cross posting like this – Ian Turton Nov 5 '15 at 17:57
  • Thanks for the heads up @iant. Didn't realize. I just deleted the post in stackoverflow. – eos Nov 5 '15 at 18:19
  • @mkennedy if I set x_0 and y_0 to 0,0 - I get wildly divergent results from basemap and pyproj. Because then basemap uses width and height to determine northings and eastings but pyproj is just using 0,0. – eos Nov 5 '15 at 19:12
1

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)
  • Ok I feel like I'm getting closer now... So, my pyproj conversion was off because I was using an invalid SRS (ie, that albers_crs dict was incorrectly specified?)... is that right? How can I correctly specify that target SRS for pyproj without using basemap and m.proj4string? So if I want to use only pyproj, how would I correctly specify this target SRS myself to get pyproj to convert lat-long to x-y accurately (ie, yielding the same result x-y values as basemap did, but without using basemap in my code)? – eos Nov 5 '15 at 19:11
  • Add in a radius, 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., 'R' : 6370997.0} – blord-castillo Nov 5 '15 at 22:23
  • Thanks @gene. Can you explain why I need to specify the 'R' parameter, and where does that specific value of 6370997 for the earth's radius come from? I never seem to see 'R' specified for Albers or Lambert etc in Proj4 strings. – eos Nov 6 '15 at 0:17
  • The conceptually clearest solution I have found is to skip the 'R' parameter and instead specify the ellipsoid (ellps) and datum as 'WGS84' in the albers_crs dict that defines the projection. Then I pass the same ellipsoid value into the Basemap function using its ellps argument. This results in pyproj and basemap producing trivially different x and y values: only off each other by 0.0000000296m and 0.0000001288m. Still I don't understand why they'd be different from each other at all. I would think I shouldn't have to specify a radius to get precise results, if I'm specifying an ellipsoid. – eos Nov 6 '15 at 1:26

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.