1

I have a list of data points given in lon/lat coordinates, which I want to convert to cartesian ones using pyproj. My dataset refers to Italy, so I was considering using either the Lambert Azimuthal Equal Area projection or the Albers Conic Equal Area one. For now I am trying out the former, but since I am not very familiar with coordinate systems, I have run into a basic question with respect to the definition of the datum, ellipsoid and earth radius.

According to what I understand so far, each ellipsoid has a standard radius value. For instance, for the WGS84 ellipsoid the radius is R=6378137. However, when I specify both the ellipsoid and R value in pyproj.proj, I get different coordinates compared to when I only define the ellipsoid.

from pyproj import Proj, transform, Geod
import numpy as np

lambert_aea1 = {'proj': 'laea',
          'lat_0':43.2947, 
          'lon_0':12.1295, 
          'x_0':0., 
          'y_0':0.,
          'ellps': 'WGS84',
          'datum': 'WGS84',
          'R':6378137.0}

lambert_aea2 = {'proj': 'laea',
          'lat_0':43.2947, 
          'lon_0':12.1295, 
          'x_0':0., 
          'y_0':0.,
          'ellps': 'WGS84',
          'datum': 'WGS84'}

xi = [12.1295, 8.4555, 11.1193, 15.8035, 13.1496]
yi = [43.2947, 44.8834, 47.1653, 41.7059, 39.4241]
inProj = Proj(init = 'epsg:4326') 
outProj1 = Proj(lambert_aea1)
outProj2 = Proj(lambert_aea2)
x1,y1 = np.array(transform(inProj,outProj1,xi,yi))
x2,y2 = np.array(transform(inProj,outProj2,xi,yi))

>> x1 = ([0., -290653.25071928, -76769.35932204, 306156.50224268, 88008.51279864])
>> x2 = ([0., -290185.4370178, -76641.501629, 305687.406043, 87878.46410642])
>> y1 = ([-21378.23435213, 161863.38128257, 409910.18924861, -191449.19444042, -451253.91439645])
>> y2 = ([0., 182924.98967511, 430536.02515148, -169835.95115785, -429281.93580092])

Now if I calculate the distances between the first point and all other points, and compare them to the distances I get by using the function pyproj.Geod.inv, I get the following differences:

geod = Geod(ellps='WGS84')
dist1 = []
for jj in range(1,len(xi)):
    dist1.append(geod.inv(xi[0],yi[0],xi[jj],yi[jj])[2])
dist2 = np.sqrt((x1[0]-x1[1:])**2 +(y1[0]-y1[1:])**2)
dist3 = np.sqrt((x2[0]-x2[1:])**2 +(y2[0]-y2[1:])**2) 
print (dist2-dist1)
print (dist3-dist1)

>> [ 517.81025651  660.65516468  485.62015613  537.69336351]
>> [ -46.78192359 -102.50618725  -38.57430363  -70.06941143]

What am I missing here?

Also are these distance "errors" reasonable? This is the first time I deal with coordinate systems and projections, so even though this precision is probably fine for my purposes, I was surprised to find such errors for my relatively small region.


To rephrase my question, I am wondering why does my 'laea' projection using pyproj change, if in addition to the ellipsoid (WGS84), I specify the radius of the earth R = 6378137. From what I understand, this is the radius dictated by the WGS84 ellipsoid anyways. Is pyproj using some default radius despite me defining the ellipsoid, or am I doing something wrong?

I also cannot find any good documentation for pyproj and proj4 with information about the arguments needed for each specific projection type.

closed as unclear what you're asking by BradHards, BERA, John Powell, mgri, AndreJ Jul 12 '17 at 13:36

Please clarify your specific problem or add additional details to highlight exactly what you need. As it's currently written, it’s hard to tell exactly what you're asking. See the How to Ask page for help clarifying this question. If this question can be reworded to fit the rules in the help center, please edit the question.

  • 1
    You'd probably be better off using a defined code (like EPSG:3035) than defining it yourself. You did this for 4326, just do it the same. – BradHards Jul 12 '17 at 9:40
  • Isn't EPSG:3035 centered on Europe though instead of Italy? I have no feel for the magnitude of difference this will have on the accuracy of the coordinates however. – thanp Jul 12 '17 at 9:55
  • It would at least give you an initial indication. You are trying an equal area projection, so it'll preserve area, not distance. By the way, pyproj.Geod.inv takes latitude and longitude pairs, not projected pairs, so the results are random. – BradHards Jul 12 '17 at 10:16
  • Thank you very much for your replies. Yes, I am entering the lon/lat pairs in the Geod.inv function. Will try with the EPSG:3035 and see what I get. The thing is though that I am interested in both areas and distances, so I was hoping that by centering the projection on my dataset I would minimize the error in the distances. – thanp Jul 12 '17 at 10:19
  • 1
    If you specify a radius, you don't have an ellipsoid anymore, but a sphere. So different results are obvious. – AndreJ Jul 12 '17 at 13:35