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

Question

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.