2

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:

enter image description here

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):
        lon, lat = radians(lon), radians(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
  • 1
    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
4

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.

These formulae can be found in this working manual from USGS at page 101. The parameters for different ellispoids (semi-major, semi-minor axis, flattening) can be found in this Wikipedia article.

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.