EPSG:3575 formulas

Question

I have been trying to find out the exact formulas to convert EPSG:3575 to EPSG:4326 (WSG84) and vice-versa using formulas described here: https://mathworld.wolfram.com/LambertAzimuthalEqual-AreaProjection.html

c := 1+math.Sin(sp)*math.Sin(lat)+math.Cos(sp)*math.Cos(lat)*math.Cos(lon-clon)
k := math.Sqrt(2.0/c)
x := radius * k * math.Cos(lat) * math.Sin(lon-clon)
y := radius * k * (math.Cos(sp)*math.Sin(lat) - math.Sin(sp)*math.Cos(lat)*math.Cos(lon-clon))

But came to no luck and I think that either the k or the radius is where I am making the mistake. The projected coordinates are almost-right, but there are small discrepancies.

My result for (45,45):

X:2795304.978972

Y:-3992109.233955

the right result that came from Proj4:

X:2804407.23

Y:-4005108.60

I am calculating the radius by taking the latitiude and calculating the earth's radius there and clon (central longitude) is 10 degrees and sp (standard parallel) is 90 degrees.

Answer 1

The formulas you are using are for the spherical case. You will have to use the ellipsoidal formulas if you want your results to match PROJ and use the WGS84 parameters:

• a = 6378137
• b = 6356752.314245
• e^2 = 1 - b^2 / a^2.

The ellipsoidal formulas are described at page 187 of Snyder's Map Projections Working Manual.

Answer 2

Figured it out, manual debugging of the PROJ library helped me, and I suggest it as the best approach to finding out the formulas, although it is based on altruism of the library and doesn't require deeper knowledge of map projections.

If anyone would want the forward formulas for EPSG:3575, here they are, extracted from the laea.cpp into Go.

lat := ToRadians(45)
lon := ToRadians(45-10)

M_HALFPI := 1.57079632679489661923
q := pj_qsfn(math.Sin(lat), e(), one_es())
qp := pj_qsfn(1., e(), one_es())

b := M_HALFPI + lat;
q = qp - q

var x,y float64

if q >= 1e-15 {
    b = math.Sqrt(q)
    x = b * math.Sin(lon)
    y =  math.Cos(lon) *  -b //b or -b, this is north polar thus -b
} else {
    y = 0
    x = 0
}

fmt.Printf("%F\r\n%f", x*6378137, y*6378137)