Merge pull request #32 from ProminentEdge/add-4087

Add in Inverse convenience method + 4087 proj support

Author: mpgerlek

Convert.go

@@ -23,11 +23,14 @@ type EPSGCode int
 
 // Supported EPSG codes
 const (
-	EPSG3395      EPSGCode = 3395
-	WorldMercator          = EPSG3395
-	EPSG3857               = 3857
-	WebMercator            = EPSG3857
-	EPSG4087               = 4087
+	EPSG3395                    EPSGCode = 3395
+	WorldMercator                        = EPSG3395
+	EPSG3857                             = 3857
+	WebMercator                          = EPSG3857
+	EPSG4087                             = 4087
+	WorldEquidistantCylindrical          = EPSG4087
+	EPSG4326                             = 4326
+	WGS84                                = EPSG4326
 )
 
 // ensure only one person is updating our cache of converters at a time
@@ -54,6 +57,26 @@ func Convert(dest EPSGCode, input []float64) ([]float64, error) {
 	return conv.convert(input)
 }
 
+// Inverse converts from a projected X/Y of a coordinate system to
+// 4326 (lat/lon, 2D).
+//
+// The input is assumed to be an array of x/y points, e.g. [x0, y0,
+// x1, y1, x2, y2, ...]. The length of the array must, therefore, be
+// even.
+//
+// The returned output is a similar array of lon/lat points, e.g. [lon0, lat0, lon1,
+// lat1, lon2, lat2, ...].
+func Inverse(src EPSGCode, input []float64) ([]float64, error) {
+	cacheLock.Lock()
+	conv, err := newConversion(src)
+	cacheLock.Unlock()
+	if err != nil {
+		return nil, err
+	}
+
+	return conv.inverse(input)
+}
+
 //---------------------------------------------------------------------------
 
 // conversion holds the objects needed to perform a conversion
@@ -70,7 +93,7 @@ var conversions = map[EPSGCode]*conversion{}
 var projStrings = map[EPSGCode]string{
 	EPSG3395: "+proj=merc +lon_0=0 +k=1 +x_0=0 +y_0=0 +datum=WGS84",                            // TODO: support +units=m +no_defs
 	EPSG3857: "+proj=merc +a=6378137 +b=6378137 +lat_ts=0.0 +lon_0=0.0 +x_0=0.0 +y_0=0 +k=1.0", // TODO: support +units=m +nadgrids=@null +wktext +no_defs
-	EPSG4087: "+proj=eqc +lat_ts=0 +lat_0=0 +lon_0=0 +x_0=0 +y_0=0 +datum=WGS84",
+	EPSG4087: "+proj=eqc +lat_ts=0 +lat_0=0 +lon_0=0 +x_0=0 +y_0=0 +datum=WGS84",               // TODO: support +units=m +no_defs
 }
 
 // newConversion creates a conversion object for the destination systems. If
@@ -147,3 +170,35 @@ func (conv *conversion) convert(input []float64) ([]float64, error) {
 
 	return output, nil
 }
+
+func (conv *conversion) inverse(input []float64) ([]float64, error) {
+	if conv == nil || conv.converter == nil {
+		return nil, fmt.Errorf("conversion not initialized")
+	}
+
+	if len(input)%2 != 0 {
+		return nil, fmt.Errorf("input array of x/y values must be an even number")
+	}
+
+	output := make([]float64, len(input))
+
+	xy := &core.CoordXY{}
+
+	for i := 0; i < len(input); i += 2 {
+		xy.X = input[i]
+		xy.Y = input[i+1]
+
+		lp, err := conv.converter.Inverse(xy)
+
+		if err != nil {
+			return nil, err
+		}
+
+		l, p := lp.Lam, lp.Phi
+
+		output[i] = support.RToDD(l)
+		output[i+1] = support.RToDD(p)
+	}
+
+	return output, nil
+}

Convert_test.go

@@ -77,6 +77,12 @@ func TestConvert(t *testing.T) {
 		outputB, err := proj.Convert(tc.dest, inputB)
 		assert.NoError(err)
 
+		invA, err := proj.Inverse(tc.dest, tc.expectedA)
+		assert.NoError(err)
+
+		invB, err := proj.Inverse(tc.dest, tc.expectedB)
+		assert.NoError(err)
+
 		const tol = 1.0e-2
 
 		for i := range tc.expectedA {
@@ -89,6 +95,77 @@ func TestConvert(t *testing.T) {
 			assert.InDelta(tc.expectedB[i], outputB[i], tol, tag)
 			assert.InDelta(tc.expectedB[i], outputB[i], tol, tag)
 		}
+
+		for i := range tc.expectedA {
+			tag := fmt.Sprintf("inverse: epsg:%d, input=A.%d", int(tc.dest), i)
+			assert.InDelta(invA[i], inputA[i], tol, tag)
+		}
+
+		for i := range tc.expectedB {
+			tag := fmt.Sprintf("inverse: epsg:%d, input=B.%d", int(tc.dest), i)
+			assert.InDelta(invB[i], inputB[i], tol, tag)
+		}
+	}
+}
+
+func TestEnsureRaisedError(t *testing.T) {
+	type testcase struct {
+		op          string
+		pt          []float64
+		expectedErr string
+		srid        proj.EPSGCode
+	}
+
+	fn := func(tc testcase) func(t *testing.T) {
+		return func(t *testing.T) {
+			var err error
+
+			if tc.op == "convert" {
+				_, err = proj.Convert(proj.EPSGCode(tc.srid), tc.pt)
+			} else {
+				_, err = proj.Inverse(proj.EPSGCode(tc.srid), tc.pt)
+			}
+
+			if err == nil {
+				t.Errorf("didn't get expected error: %v", tc.expectedErr)
+				return
+			}
+
+			if err.Error() != tc.expectedErr {
+				t.Errorf("error: %v not equal to expected error: %v", err.Error(), tc.expectedErr)
+			}
+		}
+	}
+
+	tests := map[string]testcase{
+		"3857 out of bounds WGS84": {
+			op:          "convert",
+			srid:        proj.WebMercator,
+			pt:          []float64{-180.0, 90.0},
+			expectedErr: "tolerance condition error",
+		},
+		"4326 not supported as source srid": {
+			op:          "convert",
+			srid:        proj.EPSG4326,
+			pt:          []float64{0, 0},
+			expectedErr: "epsg code is not a supported projection",
+		},
+		"convert bad point count": {
+			op:          "convert",
+			srid:        proj.WorldMercator,
+			pt:          []float64{-180.0, 90.0, 11.0},
+			expectedErr: "input array of lon/lat values must be an even number",
+		},
+		"inverse bad point count": {
+			op:          "inverse",
+			srid:        proj.WorldMercator,
+			pt:          []float64{-180.0, 90.0, 11.0},
+			expectedErr: "input array of x/y values must be an even number",
+		},
+	}
+
+	for name, tc := range tests {
+		t.Run(name, fn(tc))
 	}
 }