From d0df086ed9588f43c310da7ffea1517b7ec54501 Mon Sep 17 00:00:00 2001 From: m4rio31 Date: Wed, 14 Jan 2026 18:58:52 +0100 Subject: [PATCH 1/3] add Lambert Conformal Conic (LCC) operation --- operations/Lambert.go | 144 ++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 144 insertions(+) create mode 100644 operations/Lambert.go diff --git a/operations/Lambert.go b/operations/Lambert.go new file mode 100644 index 0000000..a808b39 --- /dev/null +++ b/operations/Lambert.go @@ -0,0 +1,144 @@ +// Copyright (C) 2018, Michael P. Gerlek (Flaxen Consulting) +// +// Portions of this code were derived from the PROJ.4 software +// In keeping with the terms of the PROJ.4 project, this software +// is provided under the MIT-style license in `LICENSE.md` and may +// additionally be subject to the copyrights of the PROJ.4 authors. + +package operations + +import ( + "math" + + "github.com/go-spatial/proj/core" + "github.com/go-spatial/proj/support" +) + +func init() { + core.RegisterConvertLPToXY("lcc", + "Lambert Conic Conformal (LCC)", + "\n\tMisc Sph, no inv.\n\tno_cut lat_b=", + NewLCC, + ) +} + +const LCCIterationEpsilon = 1e-18 + +// LCC implements core.IOperation and core.ConvertLPToXY +type LCC struct { + core.Operation + + n float64 // scale factor of the cone + F float64 // cone constant + rho0 float64 // radius at the origin parallel + lambda0 float64 // longitude of origin + phi0 float64 // latitude of origin + phi1 float64 // first standard parallel + phi2 float64 // second standard parallel + x0 float64 // offset X + y0 float64 // offset Y +} + +// NewLCC returns a new LCC +func NewLCC(system *core.System, desc *core.OperationDescription) (core.IConvertLPToXY, error) { + op := &LCC{} + op.System = system + + err := op.lccSetup(system) + if err != nil { + return nil, err + } + return op, nil +} + +// Forward Operation +func (op *LCC) Forward(lp *core.CoordLP) (*core.CoordXY, error) { + xy := &core.CoordXY{X: 0.0, Y: 0.0} + var rho float64 + + t := support.Tsfn(lp.Phi, math.Sin(lp.Phi), op.System.Ellipsoid.E) + rho = op.F * math.Pow(t, op.n) + + xy.X = rho * math.Sin(op.n*(lp.Lam)) + xy.Y = op.rho0 - (rho * math.Cos(op.n*(lp.Lam))) + + return xy, nil +} + +// Inverse Operation +func (op *LCC) Inverse(xy *core.CoordXY) (*core.CoordLP, error) { + deltaE := xy.X + deltaN := op.rho0 - xy.Y + + rPrime := math.Sqrt(deltaE*deltaE + deltaN*deltaN) + if op.n < 0 { + rPrime = -rPrime + } + + tPrime := math.Pow(rPrime/op.F, 1.0/op.n) + thetaPrime := math.Atan2(deltaE, deltaN) + + lon := (thetaPrime / op.n) - op.lambda0 + + lat := math.Pi/2.0 - 2*math.Atan(tPrime) + for range 10 { // 10 iterations limit for safety + latNew := math.Pi/2.0 - 2*math.Atan(tPrime*math.Pow((1.0-op.System.Ellipsoid.E*math.Sin(lat))/(1.0+op.System.Ellipsoid.E*math.Sin(lat)), op.System.Ellipsoid.E/2.0)) + + if math.Abs(latNew-lat) < LCCIterationEpsilon { + lat = latNew + break + } + lat = latNew + } + + return &core.CoordLP{Phi: lat, Lam: lon}, nil +} + +func (op *LCC) lccSetup(system *core.System) error { + phi0, ok0 := system.ProjString.GetAsFloat("lat_0") + if !ok0 { + phi0 = 0.0 + } + phi1, ok1 := system.ProjString.GetAsFloat("lat_1") + if !ok1 { + phi1 = 0.0 + } + phi2, ok2 := system.ProjString.GetAsFloat("lat_2") + if !ok2 { + phi2 = phi1 + } + lambda0, ok3 := system.ProjString.GetAsFloat("lon_0") + if !ok3 { + lambda0 = 0.0 + } + x0, ok4 := system.ProjString.GetAsFloat("x_0") + if !ok4 { + x0 = 0.0 + } + y0, ok5 := system.ProjString.GetAsFloat("y_0") + if !ok5 { + y0 = 0.0 + } + + op.phi0 = support.DDToR(phi0) + op.phi1 = support.DDToR(phi1) + op.phi2 = support.DDToR(phi2) + op.lambda0 = support.DDToR(lambda0) + op.x0 = x0 + op.y0 = y0 + + PE := system.Ellipsoid + + m1 := support.Msfn(math.Sin(op.phi1), math.Cos(op.phi1), PE.Es) + t1 := support.Tsfn(op.phi1, math.Sin(op.phi1), PE.E) + m2 := support.Msfn(math.Sin(op.phi2), math.Cos(op.phi2), PE.Es) + t2 := support.Tsfn(op.phi2, math.Sin(op.phi2), PE.E) + op.n = math.Log(m1/m2) / math.Log(t1/t2) + + op.F = m1 / (op.n * math.Pow(t1, op.n)) + + t0 := support.Tsfn(op.phi0, math.Sin(op.phi0), PE.E) + op.rho0 = op.F * math.Pow(t0, op.n) + + return nil +} From 56cc78b5499b0064b496f1882965855bd34b3028 Mon Sep 17 00:00:00 2001 From: m4rio31 Date: Wed, 14 Jan 2026 18:59:15 +0100 Subject: [PATCH 2/3] add test for Lambert Conformal Conic (LCC) operation --- operations/operations_test.go | 30 ++++++++++++++++++++++++++++++ 1 file changed, 30 insertions(+) diff --git a/operations/operations_test.go b/operations/operations_test.go index c76dc1b..3a4432f 100644 --- a/operations/operations_test.go +++ b/operations/operations_test.go @@ -117,6 +117,22 @@ var testdata = []data{ {-200, 100, -0.001790493, 0.000895247}, {-200, -100, -0.001790493, -0.000895247}, }, + }, { + // builtins.gie:2251 + proj: "+proj=lcc +ellps=GRS80 +lat_1=0.5 +lat_2=2", + delta: 0.1 * 0.001, + fwd: [][]float64{ + {2, 1, 222588.439735968, 110660.533870800}, + {2, -1, 222756.879700279, -110532.797660827}, + {-2, 1, -222588.439735968, 110660.533870800}, + {-2, -1, -222756.879700279, -110532.797660827}, + }, + inv: [][]float64{ + {200, 100, 0.001796359, 0.000904232}, + {200, -100, 0.001796358, -0.000904233}, + {-200, 100, -0.001796359, 0.000904232}, + {-200, -100, -0.001796358, -0.000904233}, + }, }, } @@ -215,3 +231,17 @@ func BenchmarkConvertEqc(b *testing.B) { _, _ = op.Forward(input) } } + +func BenchmarkConvertLambert(b *testing.B) { + + ps, _ := support.NewProjString("+proj=lcc +ellps=GRS80 +lat_1=0.5 +lat_2=2") + _, opx, _ := core.NewSystem(ps) + op := opx.(core.IConvertLPToXY) + input := &core.CoordLP{Lam: support.DDToR(12.0), Phi: support.DDToR(55.0)} + + b.ResetTimer() + + for i := 0; i < b.N; i++ { + _, _ = op.Forward(input) + } +} \ No newline at end of file From 20e81e72de4201d6e7f74d0fafe3ac927d2986b5 Mon Sep 17 00:00:00 2001 From: m4rio31 Date: Mon, 19 Jan 2026 19:08:10 +0100 Subject: [PATCH 3/3] add Winkel Tripel operation with test --- operations/Wintri.go | 198 ++++++++++++++++++++++++++++++++++ operations/operations_test.go | 27 ++++- 2 files changed, 223 insertions(+), 2 deletions(-) create mode 100644 operations/Wintri.go diff --git a/operations/Wintri.go b/operations/Wintri.go new file mode 100644 index 0000000..011388e --- /dev/null +++ b/operations/Wintri.go @@ -0,0 +1,198 @@ +// Copyright (C) 2018, Michael P. Gerlek (Flaxen Consulting) +// +// Portions of this code were derived from the PROJ.4 software +// In keeping with the terms of the PROJ.4 project, this software +// is provided under the MIT-style license in `LICENSE.md` and may +// additionally be subject to the copyrights of the PROJ.4 authors. + +package operations + +import ( + "math" + + "github.com/go-spatial/proj/core" + "github.com/go-spatial/proj/merror" + "github.com/go-spatial/proj/support" +) + +func init() { + core.RegisterConvertLPToXY("wintri", + "Winkel Tripel", + "\n\tPCyl., Sph.\n\tlat_1= (default: 50.467°)", + NewWintri, + ) +} + +// Wintri implements core.IOperation and core.ConvertLPToXY +type Wintri struct { + core.Operation + lat1 float64 // standard parallel (radiants) + cosLat1 float64 // cosin of standard parallel +} + +// NewWintri returns a new Winkel Tripel projection +func NewWintri(system *core.System, desc *core.OperationDescription) (core.IConvertLPToXY, error) { + op := &Wintri{} + op.System = system + + err := op.wintriSetup(system) + if err != nil { + return nil, err + } + return op, nil +} + +// Forward Operation +func (op *Wintri) Forward(lp *core.CoordLP) (*core.CoordXY, error) { + var xy core.CoordXY + + lat := lp.Phi + lon := lp.Lam + + x1 := lon * op.cosLat1 + y1 := lat + + var x2, y2 float64 + + cosLat := math.Cos(lat) + cosHalfLon := math.Cos(lon * 0.5) + alpha := math.Acos(cosLat * cosHalfLon) + + if alpha < eps10 { + x2 = lon + y2 = lat + } else { + sinAlpha := math.Sin(alpha) + if sinAlpha < eps10 { + x2 = 0.0 + y2 = 0.0 + } else { + factor := alpha / sinAlpha + x2 = 2.0 * cosLat * math.Sin(lon*0.5) * factor + y2 = math.Sin(lat) * factor + } + } + + xy.X = 0.5 * (x1 + x2) + xy.Y = 0.5 * (y1 + y2) + + return &xy, nil +} + +// Inverse Operation +func (op *Wintri) Inverse(xy *core.CoordXY) (*core.CoordLP, error) { + var lp core.CoordLP + + x := xy.X + y := xy.Y + + const maxIter = 30 + const tolerance = 1e-14 + + phi := y + lam := x / op.cosLat1 + + if phi > math.Pi*0.5 { + phi = math.Pi * 0.5 + } else if phi < -math.Pi*0.5 { + phi = -math.Pi * 0.5 + } + + if lam > math.Pi { + lam = math.Pi + } else if lam < -math.Pi { + lam = -math.Pi + } + + for range maxIter { + testLP := core.CoordLP{Phi: phi, Lam: lam} + testXY, err := op.Forward(&testLP) + if err != nil { + return nil, err + } + + dx := testXY.X - x + dy := testXY.Y - y + if math.Abs(dx) < tolerance && math.Abs(dy) < tolerance { + break + } + + if math.Abs(dx) > 10 || math.Abs(dy) > 10 { + phi = y * 0.9 + lam = x * 0.9 / op.cosLat1 + continue + } + + delta := math.Max(1e-8, math.Min(1e-6, math.Max(math.Abs(phi), math.Abs(lam))*1e-8)) + + testLP1 := core.CoordLP{Phi: phi + delta, Lam: lam} + testXY1, err1 := op.Forward(&testLP1) + if err1 != nil { + delta *= 0.5 + continue + } + dxdPhi := (testXY1.X - testXY.X) / delta + dydPhi := (testXY1.Y - testXY.Y) / delta + + testLP2 := core.CoordLP{Phi: phi, Lam: lam + delta} + testXY2, err2 := op.Forward(&testLP2) + if err2 != nil { + delta *= 0.5 + continue + } + dxdLam := (testXY2.X - testXY.X) / delta + dydLam := (testXY2.Y - testXY.Y) / delta + + det := dxdPhi*dydLam - dydPhi*dxdLam + if math.Abs(det) < 1e-15 { + return nil, merror.New(merror.ToleranceCondition, "Jacobian determinant too small in Winkel Tripel inverse") + } + + dphi := (dydLam*dx - dxdLam*dy) / det + dlam := (dxdPhi*dy - dydPhi*dx) / det + + damping := 1.0 + if math.Abs(dphi) > 0.1 || math.Abs(dlam) > 0.1 { + damping = 0.5 + } + + phi -= damping * dphi + lam -= damping * dlam + + if phi > math.Pi*0.5 { + phi = math.Pi * 0.5 + } else if phi < -math.Pi*0.5 { + phi = -math.Pi * 0.5 + } + + for lam > math.Pi { + lam -= 2 * math.Pi + } + for lam < -math.Pi { + lam += 2 * math.Pi + } + } + + lp.Phi = phi + lp.Lam = lam + + return &lp, nil +} + +func (op *Wintri) wintriSetup(system *core.System) error { + system.Ellipsoid.Es = 0.0 + op.lat1 = math.Acos(2.0 / math.Pi) + + if val, ok := system.ProjString.GetAsFloat("lat_1"); ok { + op.lat1 = support.DDToR(val) + } + + op.cosLat1 = math.Cos(op.lat1) + + if math.Abs(op.lat1) > math.Pi*0.5 { + op.lat1 = math.Copysign(math.Pi*0.5, op.lat1) + op.cosLat1 = math.Cos(op.lat1) + } + + return nil +} diff --git a/operations/operations_test.go b/operations/operations_test.go index 3a4432f..710e1ad 100644 --- a/operations/operations_test.go +++ b/operations/operations_test.go @@ -133,6 +133,16 @@ var testdata = []data{ {-200, 100, -0.001796359, 0.000904232}, {-200, -100, -0.001796358, -0.000904233}, }, + }, { + // builtins.gie:5124 + proj: "+proj=wintri +a=6400000 +lat_1=0 +lat_2=2", + delta: 0.1 * 0.001, + fwd: [][]float64{ + {2, 1, 223390.801533485, 111703.907505745}, + {2, -1, 223390.801533485, -111703.907505745}, + {-2, 1, -223390.801533485, 111703.907505745}, + {-2, -1, -223390.801533485, -111703.907505745}, + }, }, } @@ -232,7 +242,7 @@ func BenchmarkConvertEqc(b *testing.B) { } } -func BenchmarkConvertLambert(b *testing.B) { +func BenchmarkConvertLcc(b *testing.B) { ps, _ := support.NewProjString("+proj=lcc +ellps=GRS80 +lat_1=0.5 +lat_2=2") _, opx, _ := core.NewSystem(ps) @@ -244,4 +254,17 @@ func BenchmarkConvertLambert(b *testing.B) { for i := 0; i < b.N; i++ { _, _ = op.Forward(input) } -} \ No newline at end of file +} + +func BenchmarkConvertWintri(b *testing.B) { + ps, _ := support.NewProjString("+proj=wintri +a=6400000 +lat_1=0 +lat_2=2") + _, opx, _ := core.NewSystem(ps) + op := opx.(core.IConvertLPToXY) + input := &core.CoordLP{Lam: support.DDToR(12.0), Phi: support.DDToR(55.0)} + + b.ResetTimer() + + for i := 0; i < b.N; i++ { + _, _ = op.Forward(input) + } +}