From b3a18b7318d6f544ac71883f6a7a90c84af8d8bc Mon Sep 17 00:00:00 2001 From: Gautam Dey Date: Thu, 26 Sep 2019 10:40:10 -0700 Subject: [PATCH 1/5] UTM FromLngLat working Added the ability to work with UTM coordinate reference systems --- cmp/cmp.go | 5 + planar/coord/coord.go | 110 ++++++++++ planar/coord/coord_test.go | 260 ++++++++++++++++++++++ planar/coord/utm/utm.go | 408 +++++++++++++++++++++++++++++++++++ planar/coord/utm/utm_test.go | 315 +++++++++++++++++++++++++++ 5 files changed, 1098 insertions(+) create mode 100644 planar/coord/coord.go create mode 100644 planar/coord/coord_test.go create mode 100644 planar/coord/utm/utm.go create mode 100644 planar/coord/utm/utm_test.go diff --git a/cmp/cmp.go b/cmp/cmp.go index ab2f72fb..1a212c88 100644 --- a/cmp/cmp.go +++ b/cmp/cmp.go @@ -24,6 +24,11 @@ var ( NilCollection = (*geom.Collection)(nil) ) +// BitToleranceFor returns the BitToleranceFor the given tolerance +func BitToleranceFor(tol float64) int64 { + return int64(math.Float64bits(1.0+tol) - math.Float64bits(1.0)) +} + // FloatSlice compare two sets of float slices. func FloatSlice(f1, f2 []float64) bool { return Float64Slice(f1, f2, Tolerance, BitTolerance) } diff --git a/planar/coord/coord.go b/planar/coord/coord.go new file mode 100644 index 00000000..f3c1eb9b --- /dev/null +++ b/planar/coord/coord.go @@ -0,0 +1,110 @@ +package coord + +import ( + "fmt" + "log" + "math" +) + +var Precision float64 = 3.0 + + +type Interface interface { + ToLngLat() LngLat +} + +// LngLat describes Latitude,Longitude values +type LngLat struct { + // Longitude in decimal degrees + Lng float64 + // Latitude in decimal degrees + Lat float64 +} + +// NormalizeLng will insure that the longitude value is between 0-360 +func (l LngLat) NormalizeLng() LngLat { + return LngLat{ + Lat: l.Lat, + Lng: l.Lng - float64(int64((l.Lng+180.0)/360.0)*360.0), + } +} + +// LatInRadians returns the Latitude value in radians +func (l LngLat) LatInRadians() float64 { return ToRadian(l.Lat) } + +// LngInRadians returns the Longitude value in radians +func (l LngLat) LngInRadians() float64 { return ToRadian(l.Lng) } + +// LatAsDMS returns the latitude value in Degree Minute Seconds values +func (l LngLat) LatAsDMS() DMS { + log.Printf("lng: %v",l.Lat) + latD, latM, latS := toDMS(l.Lat) + latH := 'N' + if l.Lat < 0 { + latH = 'S' + } + return DMS{ + Degree: latD, + Minute: latM, + Second: latS, + Hemisphere: latH, + } +} + +// LngAsDMS returns the longitude value in Degree Minute Seconds values +func (l LngLat) LngAsDMS() DMS { + log.Printf("lng: %v",l.Lng) + lngD, lngM, lngS := toDMS(l.Lng) + lngH := 'E' + if l.Lng < 0 { + lngH = 'W' + } + return DMS{ + Degree: lngD, + Minute: lngM, + Second: lngS, + Hemisphere: lngH, + } +} + + +// ToRadian will convert a value in degree to radians +func ToRadian(degree float64) float64 { + return degree * math.Pi / 180.000 +} + +// ToDegree will convert a value in radians to degrees +func ToDegree(radian float64) float64 { + return radian * 180.000/math.Pi +} + +// Ellipsoid describes an Ellipsoid +// this may change when we get a proper projection package +type Ellipsoid struct { + Name string + Radius float64 + Eccentricity float64 + NATOCompatible bool +} + +// Convert the given lng or lat value to the degree minute seconds values +func toDMS(v float64) (d int64, m int64, s float64) { + var frac float64 + df, frac := math.Modf(v) + mf, frac := math.Modf(60.0 * frac) + s = 60.0 * frac + return int64(math.Abs(df)), int64(math.Abs(mf)), math.Abs(s) +} + +// DMS is the degree minutes and seconds +type DMS struct { + Degree int64 + Minute int64 + Second float64 + Hemisphere rune +} + +// String returns the string representation. +func (dms DMS) String() string { + return fmt.Sprintf(`%d°%d'%f"%c`, dms.Degree, dms.Minute, dms.Second, dms.Hemisphere) +} diff --git a/planar/coord/coord_test.go b/planar/coord/coord_test.go new file mode 100644 index 00000000..3662af19 --- /dev/null +++ b/planar/coord/coord_test.go @@ -0,0 +1,260 @@ +package coord + +import ( + "github.com/go-spatial/geom/cmp" + "testing" +) + +func tolerance(tol *float64) (float64, int64) { + + t := 0.0001 + if tol != nil { + t = *tol + } + return t, cmp.BitToleranceFor(t) +} + +func TestToRadianDegree(t *testing.T) { + type tcase struct { + Desc string + // Add Additonal Fields here + Degree float64 + Radian float64 + tolerance *float64 + } + + fn := func(tc tcase) func(*testing.T) { + + tol, bitTol := tolerance(tc.tolerance) + return func(t *testing.T) { + + t.Run("ToRadian", func(t *testing.T) { + rad := ToRadian(tc.Degree) + if !cmp.Float64(rad, tc.Radian, tol,bitTol) { + t.Errorf("radian, expect %v, got %v", tc.Radian, rad) + } + }) + + t.Run("ToDegree", func(t *testing.T) { + deg := ToDegree(tc.Radian) + if !cmp.Float64(deg, tc.Degree, tol,bitTol) { + t.Errorf("degree, expect %v, got %v", tc.Degree, deg) + } + }) + + } + } + + tests := []tcase{ + // Subtests + { + Degree: 0.0, + Radian: 0.0, + }, + { + Degree: 1.0, + Radian: 0.017453, + }, + { + Degree: 60.0, + Radian: 1.0472, + }, + { + Degree: 90.0, + Radian: 1.5708, + }, + { + Degree: 180.0, + Radian: 3.14159, + }, + { + Degree: 360.0, + Radian: 6.28319, + }, + { + Degree: 580.0, + Radian: 10.1229, + }, + } + + for i := range tests { + t.Run(tests[i].Desc, fn(tests[i])) + } +} + +func cmpDMS(a, b DMS) bool { + + return !(a.Degree != b.Degree || + a.Minute != b.Minute || + !cmp.Float(a.Second, b.Second) || + a.Hemisphere != b.Hemisphere) + +} + +func TestLngLat_ToDMS(t *testing.T) { + type tcase struct { + Desc string + // Add Additonal Fields here + LngLat LngLat + LngDMS DMS + LatDMS DMS + } + + fn := func(tc tcase) func(*testing.T) { + return func(t *testing.T) { + + t.Run("LatDMS", func(t *testing.T) { + dms := tc.LngLat.LatAsDMS() + if !cmpDMS(dms,tc.LatDMS) { + t.Errorf("dms, expected %v got %v", tc.LatDMS, dms) + } + }) + + t.Run("LngDMS", func(t *testing.T) { + dms := tc.LngLat.LngAsDMS() + if !cmpDMS(dms,tc.LngDMS) { + t.Errorf("dms, expected %v got %v", tc.LngDMS, dms) + } + }) + + } + } + + tests := []tcase{ + // Subtests + { + Desc: "noman's land", + LngLat: LngLat{ + Lng: 0.0, + Lat: 0.0, + }, + LngDMS: DMS{ + Degree: 0, + Minute: 0, + Second: 0.0, + Hemisphere: 'E', + }, + LatDMS: DMS{ + Degree: 0, + Minute: 0, + Second: 0.0, + Hemisphere: 'N', + }, + }, + { + Desc: "india", + LngLat: LngLat{ + Lat: 21.991952, + Lng: 78.873755, + }, + LngDMS: DMS{ + Degree: 78, + Minute: 52, + Second: 25.518, + Hemisphere: 'E', + }, + LatDMS: DMS{ + Degree: 21, + Minute: 59, + Second: 31.0272, + Hemisphere: 'N', + }, + }, + { + Desc: "zambia", + LngLat: LngLat{ + Lat: -14.723885, + Lng: 26.162606, + }, + LatDMS: DMS{ + Degree: 14, + Minute: 43, + Second: 25.986, + Hemisphere: 'S', + }, + LngDMS: DMS{ + Degree: 26, + Minute: 9, + Second: 45.3816, + Hemisphere: 'E', + }, + }, + { + Desc: "brasil", + LngLat: LngLat{ + Lat: -11.126663, + Lng: -49.038633, + }, + LatDMS: DMS{ + Degree: 11, + Minute: 7, + Second: 35.9868, + Hemisphere: 'S', + }, + LngDMS: DMS{ + Degree: 49, + Minute: 2, + Second: 19.0788, + Hemisphere: 'W', + }, + }, + { + Desc: "north canada", + LngLat: LngLat{ + Lat: 66.743373, + Lng: -102.452597, + }, + LatDMS: DMS{ + Degree: 66, + Minute: 44, + Second: 36.1428, + Hemisphere: 'N', + }, + LngDMS: DMS{ + Degree: 102, + Minute: 27, + Second: 9.3492, + Hemisphere: 'W', + }, + }, + } + + for i := range tests { + t.Run(tests[i].Desc, fn(tests[i])) + } +} + +func TestDMS_String(t *testing.T) { + type tcase struct { + Desc string + // Add Additonal Fields here + DMS DMS + Form string + } + + fn := func(tc tcase) func(*testing.T) { + return func(t *testing.T) { + str := tc.DMS.String() + if str == tc.Form { + t.Errorf("str, got %v expected %v",str, tc.Form) + } + } + } + + tests := []tcase{ + // Subtests + { + Form: `66°44'36.1428" N`, + DMS: DMS{ + Degree: 66, + Minute: 44, + Second: 36.1428, + Hemisphere:'N', + }, + }, + } + + for i := range tests { + t.Run(tests[i].Desc, fn(tests[i])) + } +} diff --git a/planar/coord/utm/utm.go b/planar/coord/utm/utm.go new file mode 100644 index 00000000..d528fe49 --- /dev/null +++ b/planar/coord/utm/utm.go @@ -0,0 +1,408 @@ +/* +Package utm provides the ability to work with UTM grids + +Ref: https://stevedutch.net/FieldMethods/UTMSystem.htm + +*/ +package utm + +import ( + "fmt" + "github.com/go-spatial/geom/planar/coord" + "log" + "math" +) + +const ( + k0 = 0.9996 // k0 - 0.9996 for UTM + l0 = 0 // λ0 = center of the map. == 0 for us. + precision = 0.0001 +) + +// UTM Zone Letters +const ( + ZoneC ZoneLetter = 'C' + ZoneD ZoneLetter = 'D' + ZoneE ZoneLetter = 'E' + ZoneF ZoneLetter = 'F' + ZoneG ZoneLetter = 'G' + ZoneH ZoneLetter = 'H' + ZoneI ZoneLetter = 'I' + ZoneJ ZoneLetter = 'J' + ZoneK ZoneLetter = 'K' + ZoneL ZoneLetter = 'L' + ZoneM ZoneLetter = 'M' + ZoneN ZoneLetter = 'N' + ZoneP ZoneLetter = 'P' + ZoneQ ZoneLetter = 'Q' + ZoneR ZoneLetter = 'R' + ZoneS ZoneLetter = 'S' + ZoneT ZoneLetter = 'T' + ZoneU ZoneLetter = 'U' + ZoneV ZoneLetter = 'V' + ZoneW ZoneLetter = 'W' + ZoneX ZoneLetter = 'X' +) + +// ZoneLetter describes the UTM zone letter +type ZoneLetter byte + +// String implements the stringer interface +func (zl ZoneLetter) String() string { return fmt.Sprintf("%c", zl) } + +// IsNorthern returns if the Zone is in the northern hemisphere +func (zl ZoneLetter) IsNorthern() bool { return zl >= 'N' } + +// IsValid will run validity check on the zone letter and number +func (zl ZoneLetter) IsValid() bool { return zl >= 'C' && zl <= 'X' && zl != 'O' } + +// quick lookup table for central meridian for each zone +// this can be calculated using the formula i := zone > 30 ? zone-31 : 30 - zone +// cmd := 3 + (6*i) +var centralMeridianDegrees = []uint{ + 3, 9, 15, 21, 27, 33, 39, 45, 51, 57, 63, 69, 75, 81, 87, 93, 99, 105, 111, 117, 123, 129, 135, 141, 147, 153, 159, 165, 171, 177, +} + +func CentralMeridian(zone Zone) (degree int, err error) { + + if !zone.IsValid() { + return 0, fmt.Errorf("zone is invalid") + } + + if zone.Number <= 30 { + idx := 30 - zone.Number + return -1 * int(centralMeridianDegrees[idx]), nil + } + idx := zone.Number - 31 + return int(centralMeridianDegrees[idx]), nil + +} + +var digraphZones = [...][2][4]rune{ + [2][4]rune{ + [4]rune{'V', 'U', 'T', 'S'}, + [4]rune{'W', 'X', 'Y', 'X'}, + }, + [2][4]rune{ + [4]rune{'D', 'C', 'B', 'A'}, + [4]rune{'E', 'F', 'G', 'H'}, + }, + [2][4]rune{ + [4]rune{'M', 'L', 'K', 'J'}, + [4]rune{'N', 'P', 'Q', 'R'}, + }, +} + +var latDigraphZones = [...]rune{'A', 'B', 'C', 'D', 'E', 'F', 'G', 'H', 'J', 'K', 'L', 'M', 'N', 'P', 'Q', 'R', 'S', 'T', 'U', 'V', 'A', 'B', 'C', 'D','E'} + +type Digraph [2]rune + +func (d Digraph) String() string { return string(d[:]) } + +func newDigraph(zone Zone, lnglat coord.LngLat) (Digraph, error) { + if !zone.IsValid() { + return Digraph{}, fmt.Errorf("zone is invalid") + } + // Howmany times have we cycled through the zones. + dZone := digraphZones[zone.Number%3] + // We already know the zone is valid + cm, _ := CentralMeridian(zone) + degreeDiff := float64(cm) - lnglat.Lng + // This is the distance in km from the centeral meridian for the zone + kmDist := int(111 * math.Cos(lnglat.LatInRadians()) * degreeDiff) + // each square is 100 km wide, so we need to see how many of these do we cross. + letterIdx := int(math.Abs(float64(kmDist / 100))) + sideSelect := 0 + if degreeDiff < 0 { + sideSelect = 1 + } + log.Printf("kmDist: %v : degreeDiff %v", kmDist, degreeDiff) + lngLetter := dZone[sideSelect][letterIdx] + + + kmDistLat := math.Abs(111.0 * lnglat.Lat) + + // if zone is even start at F + // even zones are f-z + // odd zones are a-v + offset := -1 + if zone.Number%2 == 0 { + offset = 4 // start at f + } + + log.Printf("kmDistLat: %v lat: %v", kmDistLat, lnglat.Lat) + log.Printf("km%%2000 %v",int(kmDistLat)%2000 ) + // there are 2000km per set of 20 100km blocks from the eq to the pole which is defined to be 40,000km + idx := int(math.Abs(math.Ceil( + float64(int(kmDistLat)%2000) / 100.0, + ))) + // Southern hemishere the values are laid out from the pole to the equator + if !zone.IsNorthern() { + idx = 21 - idx + } + + log.Printf("idx: %v offset: %v", idx, offset) + + letterIdx = offset + idx + latLetter := latDigraphZones[letterIdx] + + return Digraph{lngLetter, latLetter}, nil +} + +// Zone describes an UTM zone +type Zone struct { + Number int + Letter ZoneLetter +} + +// String implements the stringer interface +func (z Zone) String() string { return fmt.Sprintf("%v%v", z.Number, z.Letter) } + +// IsNorthern returns if the Zone is in the northern hemisphere +func (z Zone) IsNorthern() bool { return z.Letter.IsNorthern() } + +// IsValid will run validity check on the zone letter and number +func (z Zone) IsValid() bool { return z.Letter.IsValid() && z.Number >= 1 && z.Number <= 60 } + +// zoneNumberFromLatLng get the lat zone given the two values. +// The returned value will be from 1-60, if 0 is returned +// it means that the lat,lng value was in the polar region +// and UPS should be used. +// Transcribed from: https://github.com/gdey/GDGeoCocoa/blob/master/GDGeoCoordConv.m +func ZoneNumberFromLngLat(lnglat coord.LngLat) int { + lng, lat := lnglat.Lng, lnglat.Lat + if (lat > 84.0 && lat < 90.0) || // North Pole + (lat > -80.0 && lat < -90.0) { // South Pole + return 0 + } + + // Adjust for projects. + switch { + case lat >= 56.0 && lat < 64.0 && lng >= 3.0 && lng < 12.0: // Exceptions around Norway + return 32 + + case lat >= 72.0 && lat < 84.0: // Exceptions around Svalbard + switch { + case lng >= 0.0 && lng < 9.0: + return 31 + case lng >= 9.0 && lng < 21.0: + return 33 + case lng >= 21.0 && lng < 33.0: + return 35 + case lng >= 33.0 && lng < 42.0: + return 37 + } + } + // Recast from [-180,180) to [0,360). + // the w<-> is then divided into 60 zones from 1-60. + return int((lng+180)/6) + 1 +} + +// ZoneLetterForLat returns the UTM zone letter for the given latitude value +func ZoneLetterForLat(lat float64) (ZoneLetter, error) { + switch { + case 84 >= lat && lat >= 72: + return ZoneX, nil + + case 72 > lat && lat >= 64: + return ZoneW, nil + + case 64 > lat && lat >= 56: + return ZoneV, nil + + case 56 > lat && lat >= 48: + return ZoneU, nil + + case 48 > lat && lat >= 40: + return ZoneT, nil + + case 40 > lat && lat >= 32: + return ZoneS, nil + + case 32 > lat && lat >= 24: + return ZoneR, nil + + case 24 > lat && lat >= 16: + return ZoneQ, nil + + case 16 > lat && lat >= 8: + return ZoneP, nil + + case 8 > lat && lat >= 0: + return ZoneN, nil + + case 0 > lat && lat >= -8: + return ZoneM, nil + + case -8 > lat && lat >= -16: + return ZoneL, nil + + case -16 > lat && lat >= -24: + return ZoneK, nil + + case -24 > lat && lat >= -32: + return ZoneJ, nil + + case -32 > lat && lat >= -40: + return ZoneH, nil + + case -40 > lat && lat >= -48: + return ZoneG, nil + + case -48 > lat && lat >= -56: + return ZoneF, nil + + case -56 > lat && lat >= -64: + return ZoneE, nil + + case -64 > lat && lat >= -72: + return ZoneD, nil + + case -72 > lat && lat >= -80: + return ZoneC, nil + + default: + return 0, fmt.Errorf("latitude out of range") + } +} + +// NewZone returns the UTM zone for the given LngLat value. If there is an error, the Zone value is not valid +func NewZone(lnglat coord.LngLat) (Zone, error) { + number := ZoneNumberFromLngLat(lnglat) + letter, err := ZoneLetterForLat(lnglat.Lat) + return Zone{ + Number: number, + Letter: letter, + }, err +} + +// Coord defins an UTM coordinate +type Coord struct { + Northing float64 + Easting float64 + Zone Zone + Digraph Digraph +} + +// ScalarFactor will calculate the correct k scalar value for the given lnglat and eccentricity +func ScalarFactor(lnglat coord.LngLat, e float64) float64 { + lngRad, latRad := lnglat.LngInRadians(), lnglat.LatInRadians() + + dl := lngRad - l0 + dl2 := dl * dl + dl4 := dl2 * dl2 + dl6 := dl4 * dl2 + + e2 := (e * e) / (1 - (e * e)) + e4 := e2 * e2 + e6 := e4 * e2 + + c := math.Cos(latRad) + c2 := c * c + c4 := c2 * c2 + c6 := c4 * c2 + + t := math.Tan(latRad) + t2 := t * t + t4 := t2 * t2 + + T26 := c2 / 2 * (1 + e2*c2) + T27 := (c4 / 24) * (5 - (4 * t2) + + (24 * e2 * c2) + + (13 * e4 * c4) - + (28 * t2 * e2 * c2) + + (4 * e6 * c6) - + (48 * t2 * e4 * c4) - + (24 * t2 * e6 * c6)) + + T28 := (c6 / 720) * (61 - (148 * t2) + (16 * t4)) + + return k0 * (1 + (dl2 * T26) + (dl4 * T27) + (dl6 * T28)) +} + +// fromLngLat does the majority of the work. It assumes a valid zone and ellipsoid values +func fromLngLat(lnglat coord.LngLat, zone Zone, ellips coord.Ellipsoid) Coord { + + eccentricity, radius, nato := ellips.Eccentricity, ellips.Radius, ellips.NATOCompatible + latRad, lngRad := lnglat.LatInRadians(), lnglat.LngInRadians() + + lngOrigin := float64((zone.Number-1)*6 - 180 + 3) + lngOriginRad := coord.ToRadian(lngOrigin) + eccentPrime := eccentricity / (1 - eccentricity) + + scale := k0 + + sinLatRad := math.Sin(latRad) + + n := radius / math.Sqrt(1-eccentricity*sinLatRad*sinLatRad) + t0 := 0.0 + if latRad != 0.0 { + t0 = math.Tan(latRad) + } + cosLatRad := math.Cos(latRad) + + t := math.Pow(t0, 2.0) + c := math.Pow(eccentPrime, 2.0) * math.Pow(cosLatRad, 2.0) + a := (lngRad - lngOriginRad) * cosLatRad + + t2 := math.Pow(t, 2.0) + t3 := math.Pow(t, 3.0) + + c2 := math.Pow(c, 2.0) + + a2 := math.Pow(a, 2.0) + a3 := math.Pow(a, 3.0) + a4 := math.Pow(a, 4.0) + a5 := math.Pow(a, 5.0) + a6 := math.Pow(a, 6.0) + + e2 := math.Pow(eccentricity, 2.0) + e3 := math.Pow(eccentricity, 3.0) + + m0101 := (eccentricity / 4.0) + m0102 := (3.0 / 64.0 * e2) + m0103 := (5.0 / 256.0 * e3) + + m01 := (1 - m0101 - m0102 - m0103) * latRad + m02 := ((3.0 / 8.0 * eccentricity) + (3.0 / 32.0 * e2) + (45.0 / 1024.0 * e3)) * math.Sin(latRad*2.0) + + m0301 := math.Sin(latRad * 4.0) + m03 := ((15.0 / 256.0 * e2) + (45.0 / 1024.0 * e3)) * m0301 + + m0401 := math.Sin(latRad * 6.0) + m04 := (35.0 / 3072.0 * e3) * m0401 + + m := radius * (m01 - m02 + m03 - m04) + + easting := scale*n*(a+(1.0-t+c)*a3/6.0+(5.0-(10.0*t3)+(72.0*c)-(58.0*eccentPrime))*a5/120.0) + 500000.0 + northing := scale * (m + n*t0*(a2/2.0+ + (5.0-t+(9.0*c)+(4.0*c2))* + a4/24.0+ + (61.0-(58.0*t)+t2+(600.0*c)-(330.0*eccentPrime))* + a6/720.0)) + + if lnglat.Lat < 0.0 { + northing += 10000000.0 + } + var digraph Digraph + // only compute digraph for nato + if nato { + // we know the zone is good + digraph, _ = newDigraph(zone, lnglat) + } + return Coord{ + Northing: math.Round(northing), + Easting: math.Round(easting), + Zone: zone, + Digraph: digraph, + } +} + +// FromLngLat returns a new utm coordinate based on the provided longitude and latitude values. +func FromLngLat(lnglat coord.LngLat, ellips coord.Ellipsoid) (Coord, error) { + zone, err := NewZone(lnglat) + if err != nil { + return Coord{}, err + } + return fromLngLat(lnglat.NormalizeLng(), zone, ellips), nil +} diff --git a/planar/coord/utm/utm_test.go b/planar/coord/utm/utm_test.go new file mode 100644 index 00000000..ea93e434 --- /dev/null +++ b/planar/coord/utm/utm_test.go @@ -0,0 +1,315 @@ +package utm + +import ( + "fmt" + "github.com/go-spatial/geom/cmp" + "github.com/go-spatial/geom/planar/coord" + "strings" + "testing" + "unicode" +) + +// datum is the ellipsoid Structure for various datums. +// this is here to reduce the dependency tree. Don't count on these +// to be valid or accuract, better to use the official values +var knownEllipsoids = []coord.Ellipsoid{ + { + Name: normalizeName("Airy"), + Radius: 6377563, + Eccentricity: 0.00667054, + }, + { + Name: normalizeName("Clarke_1866"), + Radius: 6378206, + Eccentricity: 0.006768658, + }, + { + Name: normalizeName("WGS_84"), + Radius: 6378137, + Eccentricity: 0.00669438, + NATOCompatible: true, + }, +} + +func tolerance(tol *float64) (float64, int64) { + if tol != nil { + return *tol, cmp.BitToleranceFor(*tol) + } + return cmp.Tolerance, int64(cmp.BitTolerance) +} + +// normalizeName will modify the value a bit; remove trailing spaces, collapsing and transform spaces to '_' and uppercase everything else +func normalizeName(s string) string { + var str strings.Builder + s = strings.TrimSpace(s) + lastIsSpace := false + for _, r := range s { + if unicode.IsSpace(r) { + if !lastIsSpace { + lastIsSpace = true + str.WriteRune('_') + } + continue + } else { + lastIsSpace = false + } + str.WriteRune(unicode.ToUpper(r)) + } + return str.String() +} + +func getEllipsoidByName(name string) coord.Ellipsoid { + if name == "" { + name = "WGS 84" + } + name = normalizeName(name) + for _, ellp := range knownEllipsoids { + if ellp.Name == name { + return ellp + } + } + panic(fmt.Sprintf("Unknown ellipsoid: %v", name)) +} + +func TestFromLngLat(t *testing.T) { + type tcase struct { + Desc string + LngLat coord.LngLat + EllipsoidName string + Tolerance *float64 + UTM Coord + } + + fn := func(tc tcase) func(*testing.T) { + return func(t *testing.T) { + + tol, bitTol := tolerance(tc.Tolerance) + ellips := getEllipsoidByName(tc.EllipsoidName) + utm, err := FromLngLat(tc.LngLat, ellips) + // TODO(gdey): add support for tests that error + if err != nil { + t.Errorf("error, expected nil, got: %v", err) + return + } + if !cmp.Float64(utm.Northing, tc.UTM.Northing, tol, bitTol) { + t.Errorf("northing, expected %v, got: %v", tc.UTM.Northing, utm.Northing) + } + if !cmp.Float64(utm.Easting, tc.UTM.Easting, tol, bitTol) { + t.Errorf("easting, expected %v, got: %v", tc.UTM.Easting, utm.Easting) + } + if utm.Zone != tc.UTM.Zone { + t.Errorf("zone, expected %v, got: %v", tc.UTM.Zone, utm.Zone) + } + if utm.Digraph[0] != tc.UTM.Digraph[0] || utm.Digraph[1] != tc.UTM.Digraph[1] { + t.Errorf("Digraph, expected %v, got: %v", tc.UTM.Digraph, utm.Digraph) + } + } + } + + tests := []tcase{ + // Subtests + { + Desc: "Kabul", + LngLat: coord.LngLat{ + Lng: 69.1503666510912, + Lat: 34.52518357633554, + }, + UTM: Coord{ + Northing: 3820400.0, + Easting: 513800.0, + Zone: Zone{ + Letter: ZoneS, + Number: 42, + }, + Digraph: Digraph{'W', 'D'}, + }, + }, + { + Desc: "Brasil", + LngLat: coord.LngLat{ + Lng: -49.463803, + Lat: -11.126665, + }, + UTM: Coord{ + Northing: 8769581, + Easting: 667767, + Zone: Zone{ + Letter: ZoneL, + Number: 22, + }, + Digraph: Digraph{'F', 'N'}, + }, + }, + { + //https://metacpan.org/pod/Geo::Coordinates::UTM" + Desc: "perl example", + LngLat: coord.LngLat{ + Lng: -2.788951667, + Lat: 57.803055556, + }, + EllipsoidName: "Clarke_1866", + UTM: Coord{ + Northing: 6406592, + Easting: 512544, + Zone: Zone{ + Letter: ZoneV, + Number: 30, + }, + }, + }, + } + + for i := range tests { + t.Run(tests[i].Desc, fn(tests[i])) + } +} + +func TestLngLat_NormalizeLng(t *testing.T) { + type tcase struct { + Desc string + LngLat coord.LngLat + Lng float64 + } + + fn := func(tc tcase) func(*testing.T) { + + return func(t *testing.T) { + lng := tc.LngLat.NormalizeLng() + if !cmp.Float(lng.Lng, tc.Lng) { + t.Errorf("normalized lng, expected %v, got %v", tc.Lng, lng.Lng) + } + } + } + + tests := []tcase{ + // Subtests + {}, + } + + for i := range tests { + t.Run(tests[i].Desc, fn(tests[i])) + } +} + +func TestLngLat_InRadians(t *testing.T) { + type tcase struct { + Desc string + LngLat coord.LngLat + LngRadians float64 + LatRadians float64 + Tolerance *float64 + } + + fn := func(tc tcase) func(*testing.T) { + tol, bitTol := tolerance(tc.Tolerance) + return func(t *testing.T) { + + t.Run("LatInRadians", func(t *testing.T) { + lat := tc.LngLat.LatInRadians() + if !cmp.Float64(lat, tc.LatRadians, tol, bitTol) { + t.Errorf("radians, expected %v, got %v", tc.LatRadians, lat) + } + }) + + t.Run("LngInRadians", func(t *testing.T) { + lng := tc.LngLat.LngInRadians() + if !cmp.Float64(lng, tc.LngRadians, tol, bitTol) { + t.Errorf("radians, expected %v, got %v", tc.LngRadians, lng) + } + }) + + } + } + + tests := []tcase{ + // Subtests + { + LngLat: coord.LngLat{ + Lng: 69.1503666510912, + Lat: 34.5251835763355, + }, + LatRadians: 0.602578128262526, + LngRadians: 1.20690157702283, + }, + } + + for i := range tests { + t.Run(tests[i].Desc, fn(tests[i])) + } +} + +func TestNewDigraph(t *testing.T) { + type tcase struct { + Desc string + // Add Additonal Fields here + LngLat coord.LngLat + Zone *Zone + Digraph Digraph + } + + fn := func(tc tcase) func(*testing.T) { + return func(t *testing.T) { + + var ( + zone Zone + err error + ) + + if tc.Zone == nil { + zone, err = NewZone(tc.LngLat) + if err != nil { + panic("Error test LatLng not producing good zone") + } + + } else { + zone = *tc.Zone + } + + digraph, _ := newDigraph(zone, tc.LngLat) + if digraph[0] != tc.Digraph[0] || digraph[1] != tc.Digraph[1] { + t.Errorf("digraph, expected %v got %v", tc.Digraph, digraph) + } + + } + } + + tests := []tcase{ + // Subtests + { + Desc: "Green Bay", + LngLat: coord.LngLat{ + Lat: 44.438486, + Lng: -88.0, + }, + Digraph: Digraph{'D', 'Q'}, + }, + { + Desc: "Kabul", + LngLat: coord.LngLat{ + Lng: 69.1503666510912, + Lat: 34.52518357633554, + }, + Digraph: Digraph{'W', 'D'}, + }, + { + Desc: "Brasil even zone", + LngLat: coord.LngLat{ + Lng: -49.463803, + Lat: -11.126665, + }, + Digraph: Digraph{'F', 'N'}, + }, + { + Desc: "Brasil odd zone", + LngLat: coord.LngLat{ + Lat: -11.126665015021864, + Lng: -43.46380056756961 , + }, + Digraph: Digraph{'P', 'H'}, + }, + } + + for i := range tests { + t.Run(tests[i].Desc, fn(tests[i])) + } +} From e0454989bf4734a649befd3a4ba19e44e3bed517 Mon Sep 17 00:00:00 2001 From: Gautam Dey Date: Fri, 27 Sep 2019 16:02:07 -0700 Subject: [PATCH 2/5] Add ToLngLat conversion function to UTM Can not convert a UTM value to LngLat value --- planar/coord/utm/utm.go | 94 +++++++++++++++++++++++++++++++++--- planar/coord/utm/utm_test.go | 68 ++++++++++++++++++++++++-- 2 files changed, 153 insertions(+), 9 deletions(-) diff --git a/planar/coord/utm/utm.go b/planar/coord/utm/utm.go index d528fe49..9a6f9c37 100644 --- a/planar/coord/utm/utm.go +++ b/planar/coord/utm/utm.go @@ -2,15 +2,17 @@ Package utm provides the ability to work with UTM grids Ref: https://stevedutch.net/FieldMethods/UTMSystem.htm +Ref: https://gisgeography.com/central-meridian/ */ package utm import ( "fmt" - "github.com/go-spatial/geom/planar/coord" "log" "math" + + "github.com/go-spatial/geom/planar/coord" ) const ( @@ -93,7 +95,7 @@ var digraphZones = [...][2][4]rune{ }, } -var latDigraphZones = [...]rune{'A', 'B', 'C', 'D', 'E', 'F', 'G', 'H', 'J', 'K', 'L', 'M', 'N', 'P', 'Q', 'R', 'S', 'T', 'U', 'V', 'A', 'B', 'C', 'D','E'} +var latDigraphZones = [...]rune{'A', 'B', 'C', 'D', 'E', 'F', 'G', 'H', 'J', 'K', 'L', 'M', 'N', 'P', 'Q', 'R', 'S', 'T', 'U', 'V', 'A', 'B', 'C', 'D', 'E'} type Digraph [2]rune @@ -119,23 +121,22 @@ func newDigraph(zone Zone, lnglat coord.LngLat) (Digraph, error) { log.Printf("kmDist: %v : degreeDiff %v", kmDist, degreeDiff) lngLetter := dZone[sideSelect][letterIdx] - kmDistLat := math.Abs(111.0 * lnglat.Lat) // if zone is even start at F // even zones are f-z // odd zones are a-v offset := -1 - if zone.Number%2 == 0 { + if zone.Number%2 == 0 { offset = 4 // start at f } log.Printf("kmDistLat: %v lat: %v", kmDistLat, lnglat.Lat) - log.Printf("km%%2000 %v",int(kmDistLat)%2000 ) + log.Printf("km%%2000 %v", int(kmDistLat)%2000) // there are 2000km per set of 20 100km blocks from the eq to the pole which is defined to be 40,000km idx := int(math.Abs(math.Ceil( float64(int(kmDistLat)%2000) / 100.0, - ))) + ))) // Southern hemishere the values are laid out from the pole to the equator if !zone.IsNorthern() { idx = 21 - idx @@ -406,3 +407,84 @@ func FromLngLat(lnglat coord.LngLat, ellips coord.Ellipsoid) (Coord, error) { } return fromLngLat(lnglat.NormalizeLng(), zone, ellips), nil } + + +// ToLngLat transforms the utm Coord to it's Lat Lng representation based on the given datum +func (c Coord) ToLngLat(ellips coord.Ellipsoid) (coord.LngLat, error) { + + if !c.Zone.IsValid() { + return coord.LngLat{}, fmt.Errorf("invalid zone") + } + + radius, ecc := ellips.Radius, ellips.Eccentricity + x := c.Easting - 500000.0 // remove longitude offset + y := c.Northing + + if !c.Zone.IsNorthern() { + // remove Southern offset + y -= 10000000.0 + } + + ecc2 := math.Pow(ecc, 2.0) + ecc3 := math.Pow(ecc, 3.0) + + lngOrigin := float64((c.Zone.Number-1)*6 - 180 + 3) + eccPrimeSqr := (ecc / (1.0 - ecc)) + m := y / k0 + mu := m / (radius * (1.0 - (ecc / 4.0) - (3.0 / 64.0 * ecc2) - (5.0 / 256.0 * ecc3))) + + e_1 := 1.0 - ecc + e1 := (1.0 - math.Sqrt(e_1)) / (1.0 + math.Sqrt(e_1)) + e12 := math.Pow(e1, 2.0) + e13 := math.Pow(e1, 3.0) + e14 := math.Pow(e1, 4.0) + p1 := 3.0 / 2.0 * e1 + p2 := 27.0 / 32.0 * e13 + p3 := math.Sin(mu * 2.0) + p4 := 21.0 / 16.0 * e12 + p5 := 55.0 / 32.0 * e14 + p6 := math.Sin(mu * 4.0) + p7 := 151.0 / 96.0 * e13 + p8 := math.Sin(mu * 6.0) + + phi1Rad := mu + (p1-p2)*p3 + (p4-p5)*p6 + (p7)*p8 + + phi1Tan := math.Tan(phi1Rad) + phi1Sin := math.Sin(phi1Rad) + phi1Cos := math.Cos(phi1Rad) + + a := 1 - (ecc * phi1Sin * phi1Sin) + + n1 := radius / math.Sqrt(a) + t1 := math.Pow(phi1Tan, 2) + t12 := math.Pow(t1, 2) + c1 := ecc * math.Pow(phi1Cos, 2) + c12 := math.Pow(c1, 2) + c12_3 := 3 * c12 + r1 := radius * e_1 / math.Pow(a, 1.5) + d := x / (n1 * k0) + + latRad := phi1Rad - + (n1*phi1Tan/r1)* + ((math.Pow(d, 2)/2)- + (5+3*t1+10*c1-4*c12-9*eccPrimeSqr)* + math.Pow(d, 4)/24+ + (61+90*t1+298*c1+45*t12-252*eccPrimeSqr-c12_3)* + math.Pow(d, 6)*720) + + lngRad := (d - + (1+2*t1+c1)* + math.Pow(d, 3)/6 + + (5-2*c1+28*t1-c12_3+8*eccPrimeSqr+24*t12)* + math.Pow(d, 5)/120) / + phi1Cos + + return coord.LngLat{ + Lng: lngOrigin + coord.ToDegree(lngRad), + Lat: coord.ToDegree(latRad), + }, nil +} + +var x50 = int(math.Pow(10,5)) +func (utm Coord) NatoEasting() int { return int(utm.Easting)%x50 } +func (utm Coord) NatoNorthing() int { return int(utm.Northing)%x50 } diff --git a/planar/coord/utm/utm_test.go b/planar/coord/utm/utm_test.go index ea93e434..3cb0fa84 100644 --- a/planar/coord/utm/utm_test.go +++ b/planar/coord/utm/utm_test.go @@ -2,11 +2,12 @@ package utm import ( "fmt" - "github.com/go-spatial/geom/cmp" - "github.com/go-spatial/geom/planar/coord" "strings" "testing" "unicode" + + "github.com/go-spatial/geom/cmp" + "github.com/go-spatial/geom/planar/coord" ) // datum is the ellipsoid Structure for various datums. @@ -303,7 +304,7 @@ func TestNewDigraph(t *testing.T) { Desc: "Brasil odd zone", LngLat: coord.LngLat{ Lat: -11.126665015021864, - Lng: -43.46380056756961 , + Lng: -43.46380056756961, }, Digraph: Digraph{'P', 'H'}, }, @@ -313,3 +314,64 @@ func TestNewDigraph(t *testing.T) { t.Run(tests[i].Desc, fn(tests[i])) } } + +func TestCoord_ToLngLat(t *testing.T) { + type tcase struct { + Desc string + Datum string + UTM Coord + LngLat coord.LngLat + Err error + } + + fn := func(tc tcase) func(*testing.T) { + return func(t *testing.T) { + ellp := getEllipsoidByName(tc.Datum) + got, _ := tc.UTM.ToLngLat(ellp) + if !cmp.Float(got.Lng, tc.LngLat.Lng) { + t.Errorf("lng, expected %v got %v", tc.LngLat.Lng, got.Lng) + } + if !cmp.Float(got.Lat, tc.LngLat.Lat) { + t.Errorf("lat, expected %v got %v", tc.LngLat.Lat, got.Lat) + } + } + } + + tests := []tcase{ + // Subtests + { + Desc: "Kabul", + LngLat: coord.LngLat{ + Lng: 69.1503666510912, + Lat: 34.52518357633554, + }, + UTM: Coord{ + Northing: 3820400.0, + Easting: 513800.0, + Zone: Zone{ + Letter: ZoneS, + Number: 42, + }, + }, + }, + { + Desc: "Brazil", + LngLat: coord.LngLat{ + Lat: -11.126489480072872, + Lng: -43.46380056756961, + }, + UTM: Coord{ + Northing: 8769581.0, + Easting: 667767.0, + Zone: Zone{ + Letter: ZoneL, + Number: 23, + }, + }, + }, + } + + for i := range tests { + t.Run(tests[i].Desc, fn(tests[i])) + } +} From de497651b2b41abb4ccb8424f3d8e6af3aef8f58 Mon Sep 17 00:00:00 2001 From: Gautam Dey Date: Fri, 27 Sep 2019 16:08:04 -0700 Subject: [PATCH 3/5] Remove extra logging that is not needed anymore --- planar/coord/utm/utm.go | 14 ++++---------- 1 file changed, 4 insertions(+), 10 deletions(-) diff --git a/planar/coord/utm/utm.go b/planar/coord/utm/utm.go index 9a6f9c37..d5d5cb2e 100644 --- a/planar/coord/utm/utm.go +++ b/planar/coord/utm/utm.go @@ -9,7 +9,6 @@ package utm import ( "fmt" - "log" "math" "github.com/go-spatial/geom/planar/coord" @@ -118,7 +117,6 @@ func newDigraph(zone Zone, lnglat coord.LngLat) (Digraph, error) { if degreeDiff < 0 { sideSelect = 1 } - log.Printf("kmDist: %v : degreeDiff %v", kmDist, degreeDiff) lngLetter := dZone[sideSelect][letterIdx] kmDistLat := math.Abs(111.0 * lnglat.Lat) @@ -131,8 +129,6 @@ func newDigraph(zone Zone, lnglat coord.LngLat) (Digraph, error) { offset = 4 // start at f } - log.Printf("kmDistLat: %v lat: %v", kmDistLat, lnglat.Lat) - log.Printf("km%%2000 %v", int(kmDistLat)%2000) // there are 2000km per set of 20 100km blocks from the eq to the pole which is defined to be 40,000km idx := int(math.Abs(math.Ceil( float64(int(kmDistLat)%2000) / 100.0, @@ -142,8 +138,6 @@ func newDigraph(zone Zone, lnglat coord.LngLat) (Digraph, error) { idx = 21 - idx } - log.Printf("idx: %v offset: %v", idx, offset) - letterIdx = offset + idx latLetter := latDigraphZones[letterIdx] @@ -408,7 +402,6 @@ func FromLngLat(lnglat coord.LngLat, ellips coord.Ellipsoid) (Coord, error) { return fromLngLat(lnglat.NormalizeLng(), zone, ellips), nil } - // ToLngLat transforms the utm Coord to it's Lat Lng representation based on the given datum func (c Coord) ToLngLat(ellips coord.Ellipsoid) (coord.LngLat, error) { @@ -485,6 +478,7 @@ func (c Coord) ToLngLat(ellips coord.Ellipsoid) (coord.LngLat, error) { }, nil } -var x50 = int(math.Pow(10,5)) -func (utm Coord) NatoEasting() int { return int(utm.Easting)%x50 } -func (utm Coord) NatoNorthing() int { return int(utm.Northing)%x50 } +var x50 = int(math.Pow(10, 5)) + +func (utm Coord) NatoEasting() int { return int(utm.Easting) % x50 } +func (utm Coord) NatoNorthing() int { return int(utm.Northing) % x50 } From 3a770f5435942b8087d61d378f5638d8e7fe99da Mon Sep 17 00:00:00 2001 From: Gautam Dey Date: Tue, 1 Oct 2019 12:21:03 -0700 Subject: [PATCH 4/5] CodeReview Updates * Moved Tests to appropriate package * Fixed spelling errors * Improved code comments * Added missing tests --- planar/coord/coord.go | 20 +++--- planar/coord/coord_test.go | 118 ++++++++++++++++++++++++++++++----- planar/coord/utm/errors.go | 10 +++ planar/coord/utm/utm.go | 90 ++++++++++++++++++-------- planar/coord/utm/utm_test.go | 78 +---------------------- 5 files changed, 189 insertions(+), 127 deletions(-) create mode 100644 planar/coord/utm/errors.go diff --git a/planar/coord/coord.go b/planar/coord/coord.go index f3c1eb9b..de721ea9 100644 --- a/planar/coord/coord.go +++ b/planar/coord/coord.go @@ -2,13 +2,14 @@ package coord import ( "fmt" - "log" "math" ) -var Precision float64 = 3.0 - - +// Interface is an interface that wraps a ToLngLat methods +// +// ToLngLat should returns the a LngLat value in degrees that +// represents that given value as closely as possible. It is +// understood that this value may not be excate type Interface interface { ToLngLat() LngLat } @@ -37,7 +38,6 @@ func (l LngLat) LngInRadians() float64 { return ToRadian(l.Lng) } // LatAsDMS returns the latitude value in Degree Minute Seconds values func (l LngLat) LatAsDMS() DMS { - log.Printf("lng: %v",l.Lat) latD, latM, latS := toDMS(l.Lat) latH := 'N' if l.Lat < 0 { @@ -53,7 +53,6 @@ func (l LngLat) LatAsDMS() DMS { // LngAsDMS returns the longitude value in Degree Minute Seconds values func (l LngLat) LngAsDMS() DMS { - log.Printf("lng: %v",l.Lng) lngD, lngM, lngS := toDMS(l.Lng) lngH := 'E' if l.Lng < 0 { @@ -67,7 +66,6 @@ func (l LngLat) LngAsDMS() DMS { } } - // ToRadian will convert a value in degree to radians func ToRadian(degree float64) float64 { return degree * math.Pi / 180.000 @@ -75,15 +73,15 @@ func ToRadian(degree float64) float64 { // ToDegree will convert a value in radians to degrees func ToDegree(radian float64) float64 { - return radian * 180.000/math.Pi + return radian * 180.000 / math.Pi } // Ellipsoid describes an Ellipsoid // this may change when we get a proper projection package type Ellipsoid struct { - Name string - Radius float64 - Eccentricity float64 + Name string + Radius float64 + Eccentricity float64 NATOCompatible bool } diff --git a/planar/coord/coord_test.go b/planar/coord/coord_test.go index 3662af19..4c4cba29 100644 --- a/planar/coord/coord_test.go +++ b/planar/coord/coord_test.go @@ -1,8 +1,9 @@ package coord import ( - "github.com/go-spatial/geom/cmp" "testing" + + "github.com/go-spatial/geom/cmp" ) func tolerance(tol *float64) (float64, int64) { @@ -30,14 +31,14 @@ func TestToRadianDegree(t *testing.T) { t.Run("ToRadian", func(t *testing.T) { rad := ToRadian(tc.Degree) - if !cmp.Float64(rad, tc.Radian, tol,bitTol) { + if !cmp.Float64(rad, tc.Radian, tol, bitTol) { t.Errorf("radian, expect %v, got %v", tc.Radian, rad) } }) t.Run("ToDegree", func(t *testing.T) { deg := ToDegree(tc.Radian) - if !cmp.Float64(deg, tc.Degree, tol,bitTol) { + if !cmp.Float64(deg, tc.Degree, tol, bitTol) { t.Errorf("degree, expect %v, got %v", tc.Degree, deg) } }) @@ -105,14 +106,14 @@ func TestLngLat_ToDMS(t *testing.T) { t.Run("LatDMS", func(t *testing.T) { dms := tc.LngLat.LatAsDMS() - if !cmpDMS(dms,tc.LatDMS) { + if !cmpDMS(dms, tc.LatDMS) { t.Errorf("dms, expected %v got %v", tc.LatDMS, dms) } }) t.Run("LngDMS", func(t *testing.T) { dms := tc.LngLat.LngAsDMS() - if !cmpDMS(dms,tc.LngDMS) { + if !cmpDMS(dms, tc.LngDMS) { t.Errorf("dms, expected %v got %v", tc.LngDMS, dms) } }) @@ -175,7 +176,7 @@ func TestLngLat_ToDMS(t *testing.T) { LngDMS: DMS{ Degree: 26, Minute: 9, - Second: 45.3816, + Second: 45.3816, Hemisphere: 'E', }, }, @@ -194,7 +195,7 @@ func TestLngLat_ToDMS(t *testing.T) { LngDMS: DMS{ Degree: 49, Minute: 2, - Second: 19.0788, + Second: 19.0788, Hemisphere: 'W', }, }, @@ -213,7 +214,7 @@ func TestLngLat_ToDMS(t *testing.T) { LngDMS: DMS{ Degree: 102, Minute: 27, - Second: 9.3492, + Second: 9.3492, Hemisphere: 'W', }, }, @@ -228,7 +229,7 @@ func TestDMS_String(t *testing.T) { type tcase struct { Desc string // Add Additonal Fields here - DMS DMS + DMS DMS Form string } @@ -236,7 +237,7 @@ func TestDMS_String(t *testing.T) { return func(t *testing.T) { str := tc.DMS.String() if str == tc.Form { - t.Errorf("str, got %v expected %v",str, tc.Form) + t.Errorf("str, got %v expected %v", str, tc.Form) } } } @@ -246,11 +247,100 @@ func TestDMS_String(t *testing.T) { { Form: `66°44'36.1428" N`, DMS: DMS{ - Degree: 66, - Minute: 44, - Second: 36.1428, - Hemisphere:'N', + Degree: 66, + Minute: 44, + Second: 36.1428, + Hemisphere: 'N', + }, + }, + } + + for i := range tests { + t.Run(tests[i].Desc, fn(tests[i])) + } +} + +func TestLngLat_InRadians(t *testing.T) { + type tcase struct { + Desc string + LngLat LngLat + LngRadians float64 + LatRadians float64 + Tolerance *float64 + } + + fn := func(tc tcase) func(*testing.T) { + tol, bitTol := tolerance(tc.Tolerance) + return func(t *testing.T) { + + t.Run("LatInRadians", func(t *testing.T) { + lat := tc.LngLat.LatInRadians() + if !cmp.Float64(lat, tc.LatRadians, tol, bitTol) { + t.Errorf("radians, expected %v, got %v", tc.LatRadians, lat) + } + }) + + t.Run("LngInRadians", func(t *testing.T) { + lng := tc.LngLat.LngInRadians() + if !cmp.Float64(lng, tc.LngRadians, tol, bitTol) { + t.Errorf("radians, expected %v, got %v", tc.LngRadians, lng) + } + }) + + } + } + + tests := []tcase{ + // Subtests + { + LngLat: LngLat{ + Lng: 69.1503666510912, + Lat: 34.5251835763355, + }, + LatRadians: 0.602578128262526, + LngRadians: 1.20690157702283, + }, + } + + for i := range tests { + t.Run(tests[i].Desc, fn(tests[i])) + } +} + +func TestLngLat_NormalizeLng(t *testing.T) { + type tcase struct { + Desc string + LngLat LngLat + Lng float64 + } + + fn := func(tc tcase) func(*testing.T) { + + return func(t *testing.T) { + lng := tc.LngLat.NormalizeLng() + if !cmp.Float(lng.Lng, tc.Lng) { + t.Errorf("normalized lng, expected %v, got %v", tc.Lng, lng.Lng) + } + } + } + + tests := []tcase{ + { + Desc: "Brasil", + LngLat: LngLat{ + Lng: -49.463803, + Lat: -11.126665, + }, + Lng: -49.463803, + }, + { + + Desc: "Kabul", + LngLat: LngLat{ + Lng: 69.1503666510912, + Lat: 34.52518357633554, }, + Lng: 69.1503666510912, }, } diff --git a/planar/coord/utm/errors.go b/planar/coord/utm/errors.go new file mode 100644 index 00000000..0a3245d1 --- /dev/null +++ b/planar/coord/utm/errors.go @@ -0,0 +1,10 @@ +package utm + +import "github.com/gdey/errors" + +const ( + // ErrInvalidZone will be return if the given zone is invalid + ErrInvalidZone = errors.String("zone is invalid") + // ErrLatitudeOutOfRange will be returned if the latitude is not in the correct range of acceptable values + ErrLatitudeOutOfRange = errors.String("latitude out of range") +) diff --git a/planar/coord/utm/utm.go b/planar/coord/utm/utm.go index d5d5cb2e..b6763952 100644 --- a/planar/coord/utm/utm.go +++ b/planar/coord/utm/utm.go @@ -1,8 +1,9 @@ /* -Package utm provides the ability to work with UTM grids +Package utm provides the ability to work with UTM coordinates -Ref: https://stevedutch.net/FieldMethods/UTMSystem.htm -Ref: https://gisgeography.com/central-meridian/ + References: + https://stevedutch.net/FieldMethods/UTMSystem.htm + https://gisgeography.com/central-meridian/ */ package utm @@ -15,9 +16,8 @@ import ( ) const ( - k0 = 0.9996 // k0 - 0.9996 for UTM - l0 = 0 // λ0 = center of the map. == 0 for us. - precision = 0.0001 + k0 = 0.9996 // k0 - 0.9996 for UTM + l0 = 0 // λ0 = center of the map. == 0 for us. ) // UTM Zone Letters @@ -58,16 +58,25 @@ func (zl ZoneLetter) IsNorthern() bool { return zl >= 'N' } func (zl ZoneLetter) IsValid() bool { return zl >= 'C' && zl <= 'X' && zl != 'O' } // quick lookup table for central meridian for each zone -// this can be calculated using the formula i := zone > 30 ? zone-31 : 30 - zone -// cmd := 3 + (6*i) +// this can be calculated using the formula: +// +// ⎧ 30 - zone if zone <= 30 +// i ⎨ zone - 30 if zone > 30 +// ⎩ +// centralMeridianDegrees( zone ) = 3 + 6i +// var centralMeridianDegrees = []uint{ 3, 9, 15, 21, 27, 33, 39, 45, 51, 57, 63, 69, 75, 81, 87, 93, 99, 105, 111, 117, 123, 129, 135, 141, 147, 153, 159, 165, 171, 177, } +// CentralMeridian returns the central meridian degree for the given zone. +// +// Possible errors: +// ErrInvalidZone func CentralMeridian(zone Zone) (degree int, err error) { if !zone.IsValid() { - return 0, fmt.Errorf("zone is invalid") + return 0, ErrInvalidZone } if zone.Number <= 30 { @@ -79,7 +88,9 @@ func CentralMeridian(zone Zone) (degree int, err error) { } -var digraphZones = [...][2][4]rune{ +// lngDigraphZones are the zone labels for the longitudinal values. The labels are split in middle +// so that one can use the central medial to figure grouping to use +var lngDigraphZones = [...][2][4]rune{ [2][4]rune{ [4]rune{'V', 'U', 'T', 'S'}, [4]rune{'W', 'X', 'Y', 'X'}, @@ -94,23 +105,36 @@ var digraphZones = [...][2][4]rune{ }, } +// latDigraphZones are the zone labels for the latitudinal values. var latDigraphZones = [...]rune{'A', 'B', 'C', 'D', 'E', 'F', 'G', 'H', 'J', 'K', 'L', 'M', 'N', 'P', 'Q', 'R', 'S', 'T', 'U', 'V', 'A', 'B', 'C', 'D', 'E'} +// Digraph is the NATO diagraph type Digraph [2]rune +// String returns the string representation of a diagraph func (d Digraph) String() string { return string(d[:]) } +// newDigraph will calculate the digraph based on the provided zone and lnglat values. +// only error returned is ErrInvalidZone +// +// Reference: https://stevedutch.net/FieldMethods/UTMSystem.htm func newDigraph(zone Zone, lnglat coord.LngLat) (Digraph, error) { if !zone.IsValid() { - return Digraph{}, fmt.Errorf("zone is invalid") + return Digraph{}, ErrInvalidZone } - // Howmany times have we cycled through the zones. - dZone := digraphZones[zone.Number%3] + + // How many times have we cycled through the zones. + // There are only 3 zones total + dZone := lngDigraphZones[zone.Number%3] + // We already know the zone is valid cm, _ := CentralMeridian(zone) degreeDiff := float64(cm) - lnglat.Lng - // This is the distance in km from the centeral meridian for the zone - kmDist := int(111 * math.Cos(lnglat.LatInRadians()) * degreeDiff) + + // This is the distance in km from the central meridian for the zone + // Note by the definition of a zone, we know there are 111km per degree. + kmDist := int(111 * degreeDiff * math.Cos(lnglat.LatInRadians())) + // each square is 100 km wide, so we need to see how many of these do we cross. letterIdx := int(math.Abs(float64(kmDist / 100))) sideSelect := 0 @@ -133,7 +157,7 @@ func newDigraph(zone Zone, lnglat coord.LngLat) (Digraph, error) { idx := int(math.Abs(math.Ceil( float64(int(kmDistLat)%2000) / 100.0, ))) - // Southern hemishere the values are laid out from the pole to the equator + // Southern hemisphere the values are laid out from the pole to the equator if !zone.IsNorthern() { idx = 21 - idx } @@ -159,11 +183,14 @@ func (z Zone) IsNorthern() bool { return z.Letter.IsNorthern() } // IsValid will run validity check on the zone letter and number func (z Zone) IsValid() bool { return z.Letter.IsValid() && z.Number >= 1 && z.Number <= 60 } -// zoneNumberFromLatLng get the lat zone given the two values. -// The returned value will be from 1-60, if 0 is returned -// it means that the lat,lng value was in the polar region -// and UPS should be used. -// Transcribed from: https://github.com/gdey/GDGeoCocoa/blob/master/GDGeoCoordConv.m +// ZoneNumberFromLngLat will get the zone number for the given LngLat value. +// +// The returned value will be from 1-60. +// If 0 is returned it means that the lat,lng value was in the polar region +// and UPS should be used instead. +// +// Transcribed from: +// https://github.com/gdey/GDGeoCocoa/blob/master/GDGeoCoordConv.m func ZoneNumberFromLngLat(lnglat coord.LngLat) int { lng, lat := lnglat.Lng, lnglat.Lat if (lat > 84.0 && lat < 90.0) || // North Pole @@ -194,6 +221,9 @@ func ZoneNumberFromLngLat(lnglat coord.LngLat) int { } // ZoneLetterForLat returns the UTM zone letter for the given latitude value +// +// Possible errors: +// ErrLatitudeOutOfRange func ZoneLetterForLat(lat float64) (ZoneLetter, error) { switch { case 84 >= lat && lat >= 72: @@ -257,11 +287,14 @@ func ZoneLetterForLat(lat float64) (ZoneLetter, error) { return ZoneC, nil default: - return 0, fmt.Errorf("latitude out of range") + return 0, ErrLatitudeOutOfRange } } -// NewZone returns the UTM zone for the given LngLat value. If there is an error, the Zone value is not valid +// NewZone returns the UTM zone for the given LngLat value. +// +// Possible errors: +// ErrLatitudeOutOfRange func NewZone(lnglat coord.LngLat) (Zone, error) { number := ZoneNumberFromLngLat(lnglat) letter, err := ZoneLetterForLat(lnglat.Lat) @@ -315,7 +348,9 @@ func ScalarFactor(lnglat coord.LngLat, e float64) float64 { return k0 * (1 + (dl2 * T26) + (dl4 * T27) + (dl6 * T28)) } -// fromLngLat does the majority of the work. It assumes a valid zone and ellipsoid values +// fromLngLat does the majority of the work. +// +// It assumes a valid zone and ellipsoid values func fromLngLat(lnglat coord.LngLat, zone Zone, ellips coord.Ellipsoid) Coord { eccentricity, radius, nato := ellips.Eccentricity, ellips.Radius, ellips.NATOCompatible @@ -480,5 +515,8 @@ func (c Coord) ToLngLat(ellips coord.Ellipsoid) (coord.LngLat, error) { var x50 = int(math.Pow(10, 5)) -func (utm Coord) NatoEasting() int { return int(utm.Easting) % x50 } -func (utm Coord) NatoNorthing() int { return int(utm.Northing) % x50 } +// NatoEasting returns the easting value for NATO +func (c Coord) NatoEasting() int { return int(c.Easting) % x50 } + +// NatoNorthing returns the northing value for NATO +func (c Coord) NatoNorthing() int { return int(c.Northing) % x50 } diff --git a/planar/coord/utm/utm_test.go b/planar/coord/utm/utm_test.go index 3cb0fa84..3e2a21f8 100644 --- a/planar/coord/utm/utm_test.go +++ b/planar/coord/utm/utm_test.go @@ -12,7 +12,7 @@ import ( // datum is the ellipsoid Structure for various datums. // this is here to reduce the dependency tree. Don't count on these -// to be valid or accuract, better to use the official values +// to be valid or accurate, better to use the official values var knownEllipsoids = []coord.Ellipsoid{ { Name: normalizeName("Airy"), @@ -165,84 +165,10 @@ func TestFromLngLat(t *testing.T) { } } -func TestLngLat_NormalizeLng(t *testing.T) { - type tcase struct { - Desc string - LngLat coord.LngLat - Lng float64 - } - - fn := func(tc tcase) func(*testing.T) { - - return func(t *testing.T) { - lng := tc.LngLat.NormalizeLng() - if !cmp.Float(lng.Lng, tc.Lng) { - t.Errorf("normalized lng, expected %v, got %v", tc.Lng, lng.Lng) - } - } - } - - tests := []tcase{ - // Subtests - {}, - } - - for i := range tests { - t.Run(tests[i].Desc, fn(tests[i])) - } -} - -func TestLngLat_InRadians(t *testing.T) { - type tcase struct { - Desc string - LngLat coord.LngLat - LngRadians float64 - LatRadians float64 - Tolerance *float64 - } - - fn := func(tc tcase) func(*testing.T) { - tol, bitTol := tolerance(tc.Tolerance) - return func(t *testing.T) { - - t.Run("LatInRadians", func(t *testing.T) { - lat := tc.LngLat.LatInRadians() - if !cmp.Float64(lat, tc.LatRadians, tol, bitTol) { - t.Errorf("radians, expected %v, got %v", tc.LatRadians, lat) - } - }) - - t.Run("LngInRadians", func(t *testing.T) { - lng := tc.LngLat.LngInRadians() - if !cmp.Float64(lng, tc.LngRadians, tol, bitTol) { - t.Errorf("radians, expected %v, got %v", tc.LngRadians, lng) - } - }) - - } - } - - tests := []tcase{ - // Subtests - { - LngLat: coord.LngLat{ - Lng: 69.1503666510912, - Lat: 34.5251835763355, - }, - LatRadians: 0.602578128262526, - LngRadians: 1.20690157702283, - }, - } - - for i := range tests { - t.Run(tests[i].Desc, fn(tests[i])) - } -} - func TestNewDigraph(t *testing.T) { type tcase struct { Desc string - // Add Additonal Fields here + // Add Additional Fields here LngLat coord.LngLat Zone *Zone Digraph Digraph From 3387d8262f747c9f60d60f4e702324c1849b5c54 Mon Sep 17 00:00:00 2001 From: Gautam Dey Date: Tue, 1 Oct 2019 12:27:19 -0700 Subject: [PATCH 5/5] Quick comment and spelling error fixes --- planar/coord/coord_test.go | 11 +++++------ planar/coord/utm/utm.go | 4 ++-- 2 files changed, 7 insertions(+), 8 deletions(-) diff --git a/planar/coord/coord_test.go b/planar/coord/coord_test.go index 4c4cba29..35714c66 100644 --- a/planar/coord/coord_test.go +++ b/planar/coord/coord_test.go @@ -18,15 +18,15 @@ func tolerance(tol *float64) (float64, int64) { func TestToRadianDegree(t *testing.T) { type tcase struct { Desc string - // Add Additonal Fields here + // Add Additional Fields here Degree float64 Radian float64 - tolerance *float64 + Tolerance *float64 } fn := func(tc tcase) func(*testing.T) { - tol, bitTol := tolerance(tc.tolerance) + tol, bitTol := tolerance(tc.Tolerance) return func(t *testing.T) { t.Run("ToRadian", func(t *testing.T) { @@ -95,7 +95,7 @@ func cmpDMS(a, b DMS) bool { func TestLngLat_ToDMS(t *testing.T) { type tcase struct { Desc string - // Add Additonal Fields here + // Add Additional Fields here LngLat LngLat LngDMS DMS LatDMS DMS @@ -228,7 +228,7 @@ func TestLngLat_ToDMS(t *testing.T) { func TestDMS_String(t *testing.T) { type tcase struct { Desc string - // Add Additonal Fields here + // Add Additional Fields here DMS DMS Form string } @@ -291,7 +291,6 @@ func TestLngLat_InRadians(t *testing.T) { } tests := []tcase{ - // Subtests { LngLat: LngLat{ Lng: 69.1503666510912, diff --git a/planar/coord/utm/utm.go b/planar/coord/utm/utm.go index b6763952..34cad81e 100644 --- a/planar/coord/utm/utm.go +++ b/planar/coord/utm/utm.go @@ -304,7 +304,7 @@ func NewZone(lnglat coord.LngLat) (Zone, error) { }, err } -// Coord defins an UTM coordinate +// Coord defines an UTM coordinate type Coord struct { Northing float64 Easting float64 @@ -350,7 +350,7 @@ func ScalarFactor(lnglat coord.LngLat, e float64) float64 { // fromLngLat does the majority of the work. // -// It assumes a valid zone and ellipsoid values +// Valid zone and ellipsoid values assumed. func fromLngLat(lnglat coord.LngLat, zone Zone, ellips coord.Ellipsoid) Coord { eccentricity, radius, nato := ellips.Eccentricity, ellips.Radius, ellips.NATOCompatible