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..de721ea9 --- /dev/null +++ b/planar/coord/coord.go @@ -0,0 +1,108 @@ +package coord + +import ( + "fmt" + "math" +) + +// 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 +} + +// 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 { + 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 { + 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..35714c66 --- /dev/null +++ b/planar/coord/coord_test.go @@ -0,0 +1,349 @@ +package coord + +import ( + "testing" + + "github.com/go-spatial/geom/cmp" +) + +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 Additional 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 Additional 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 Additional 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])) + } +} + +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{ + { + 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, + }, + } + + for i := range tests { + t.Run(tests[i].Desc, fn(tests[i])) + } +} 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 new file mode 100644 index 00000000..34cad81e --- /dev/null +++ b/planar/coord/utm/utm.go @@ -0,0 +1,522 @@ +/* +Package utm provides the ability to work with UTM coordinates + + References: + https://stevedutch.net/FieldMethods/UTMSystem.htm + https://gisgeography.com/central-meridian/ + +*/ +package utm + +import ( + "fmt" + "math" + + "github.com/go-spatial/geom/planar/coord" +) + +const ( + k0 = 0.9996 // k0 - 0.9996 for UTM + l0 = 0 // λ0 = center of the map. == 0 for us. +) + +// 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: +// +// ⎧ 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, ErrInvalidZone + } + + if zone.Number <= 30 { + idx := 30 - zone.Number + return -1 * int(centralMeridianDegrees[idx]), nil + } + idx := zone.Number - 31 + return int(centralMeridianDegrees[idx]), nil + +} + +// 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'}, + }, + [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'}, + }, +} + +// 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{}, ErrInvalidZone + } + + // 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 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 + if degreeDiff < 0 { + sideSelect = 1 + } + 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 + } + + // 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 hemisphere the values are laid out from the pole to the equator + if !zone.IsNorthern() { + idx = 21 - idx + } + + 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 } + +// 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 + (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 +// +// Possible errors: +// ErrLatitudeOutOfRange +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, ErrLatitudeOutOfRange + } +} + +// 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) + return Zone{ + Number: number, + Letter: letter, + }, err +} + +// Coord defines 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. +// +// 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 + 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 +} + +// 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)) + +// 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 new file mode 100644 index 00000000..3e2a21f8 --- /dev/null +++ b/planar/coord/utm/utm_test.go @@ -0,0 +1,303 @@ +package utm + +import ( + "fmt" + "strings" + "testing" + "unicode" + + "github.com/go-spatial/geom/cmp" + "github.com/go-spatial/geom/planar/coord" +) + +// 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 accurate, 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 TestNewDigraph(t *testing.T) { + type tcase struct { + Desc string + // Add Additional 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])) + } +} + +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])) + } +}