Lambert operation

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
+}

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)
+	}
+}