# Trying to implement Albers projection in Python, what is going wrong?

I have a platform where I can't import Proj4. I only need to convert points to Albers projection, so it seemed simple, but my conversion has some errors:

Projection definition: `+proj=aea +lat_1=52 +lat_2=64 +lat_0=0 +lon_0=105 +x_0=18500000 +y_0=0 +ellps=krass +units=m +towgs84=28,-130,-95,0,0,0,0 +no_defs'`

Here's the code implemented from the formulae on Wikipedia: ``````from math import sin, radians, cos

def albers_proj(lon0, fi0, fi1, fi2, R=6367444.5, x0=0):
lon0, fi0, fi1, fi2 = [radians(i) for i in (lon0, fi0, fi1, fi2)]

n = (sin(fi1) + sin(fi2)) / 2
C = cos(fi1) ** 2 + 2 * n * sin(fi1)
ro0 = R / n * (C - 2 * n * sin(fi0)) ** .5
def transformer(lon, lat):
ro = R / n * (C - 2 * n * sin(lat)) ** .5
theta = n * (lon - lon0)

return ro * sin(theta) + x0, ro0 - ro * cos(theta)

return transformer

albers_siberia = albers_proj(105, 0, 52, 64, R=6367444.5, x0=18500000)
print(albers_siberia(83, 55))
``````

Output:

``````(17127237.118954152, 5777762.666018201)
``````

I found the best R to minimise the error, but still something is not right:

``````    wgs84: 83.0 55.0

r = 6356302.97445488

qgis: 17121825.913615   5767652.8720776
my script: 17129639.13010738 5767652.944553478
error:     7813.21649238       0.072475878
``````

What am I doing wrong?

• These formulas are for the Sphere, but are using an ellipsoid in the Proj4 string : `+ellps=krass`. – FSimardGIS Nov 14 '18 at 15:16
• Formulas for the ellipsoid can be found at page 101 of this working manual on USGS. – FSimardGIS Nov 14 '18 at 15:19
• @FSimardGIS post this as an answer, I'll accept it. Thanks a lot! – culebrón Nov 14 '18 at 15:20

The formulae you have used are for the spherical case. However, your proj4 string indicates that you are using an ellipsoid (`+ellps=krass`). In order to get the proper results on the ellipsoid, you will have to implement the ellipsoidal formulae.