# Projecting data from 'EPSG:4326' to Lambert Azimuthal Equal Area with PyProj [closed]

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,yi,xi[jj],yi[jj]))
dist2 = np.sqrt((x1-x1[1:])**2 +(y1-y1[1:])**2)
dist3 = np.sqrt((x2-x2[1:])**2 +(y2-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.

• 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
• 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