From 5410a3d015cef61a5a4d6bc9c39275f108efbe9e Mon Sep 17 00:00:00 2001 From: Patrick Yukman Date: Mon, 6 Nov 2023 14:06:01 -0800 Subject: [PATCH 1/6] minor import cleanup + go mod tidy --- geodesic/direct_test.go | 1 - geodesic/test_cases.go | 1 - geodtest/geod_test.go | 1 - go.mod | 6 +++--- go.sum | 14 ++++++++++---- readme/readme_test.go | 1 - 6 files changed, 13 insertions(+), 11 deletions(-) diff --git a/geodesic/direct_test.go b/geodesic/direct_test.go index dec3637..7592e1c 100644 --- a/geodesic/direct_test.go +++ b/geodesic/direct_test.go @@ -5,7 +5,6 @@ import ( "testing" "github.com/pymaxion/geographiclib-go/geodesic/capabilities" - "github.com/stretchr/testify/assert" ) diff --git a/geodesic/test_cases.go b/geodesic/test_cases.go index 20a83fb..2f30fe2 100644 --- a/geodesic/test_cases.go +++ b/geodesic/test_cases.go @@ -6,7 +6,6 @@ import ( "testing" "github.com/pymaxion/geographiclib-go/geodesic/capabilities" - "github.com/stretchr/testify/assert" "github.com/stretchr/testify/require" ) diff --git a/geodtest/geod_test.go b/geodtest/geod_test.go index 3d89431..6cccac2 100644 --- a/geodtest/geod_test.go +++ b/geodtest/geod_test.go @@ -10,7 +10,6 @@ import ( "testing" "github.com/pymaxion/geographiclib-go/geodesic" - "github.com/stretchr/testify/assert" ) diff --git a/go.mod b/go.mod index e09dfa6..d147141 100644 --- a/go.mod +++ b/go.mod @@ -2,10 +2,10 @@ module github.com/pymaxion/geographiclib-go go 1.17 -require github.com/stretchr/testify v1.7.0 +require github.com/stretchr/testify v1.8.4 require ( - github.com/davecgh/go-spew v1.1.0 // indirect + github.com/davecgh/go-spew v1.1.1 // indirect github.com/pmezard/go-difflib v1.0.0 // indirect - gopkg.in/yaml.v3 v3.0.0-20200313102051-9f266ea9e77c // indirect + gopkg.in/yaml.v3 v3.0.1 // indirect ) diff --git a/go.sum b/go.sum index acb88a4..479781e 100644 --- a/go.sum +++ b/go.sum @@ -1,11 +1,17 @@ -github.com/davecgh/go-spew v1.1.0 h1:ZDRjVQ15GmhC3fiQ8ni8+OwkZQO4DARzQgrnXU1Liz8= github.com/davecgh/go-spew v1.1.0/go.mod h1:J7Y8YcW2NihsgmVo/mv3lAwl/skON4iLHjSsI+c5H38= +github.com/davecgh/go-spew v1.1.1 h1:vj9j/u1bqnvCEfJOwUhtlOARqs3+rkHYY13jYWTU97c= +github.com/davecgh/go-spew v1.1.1/go.mod h1:J7Y8YcW2NihsgmVo/mv3lAwl/skON4iLHjSsI+c5H38= github.com/pmezard/go-difflib v1.0.0 h1:4DBwDE0NGyQoBHbLQYPwSUPoCMWR5BEzIk/f1lZbAQM= github.com/pmezard/go-difflib v1.0.0/go.mod h1:iKH77koFhYxTK1pcRnkKkqfTogsbg7gZNVY4sRDYZ/4= github.com/stretchr/objx v0.1.0/go.mod h1:HFkY916IF+rwdDfMAkV7OtwuqBVzrE8GR6GFx+wExME= -github.com/stretchr/testify v1.7.0 h1:nwc3DEeHmmLAfoZucVR881uASk0Mfjw8xYJ99tb5CcY= -github.com/stretchr/testify v1.7.0/go.mod h1:6Fq8oRcR53rry900zMqJjRRixrwX3KX962/h/Wwjteg= +github.com/stretchr/objx v0.4.0/go.mod h1:YvHI0jy2hoMjB+UWwv71VJQ9isScKT/TqJzVSSt89Yw= +github.com/stretchr/objx v0.5.0/go.mod h1:Yh+to48EsGEfYuaHDzXPcE3xhTkx73EhmCGUpEOglKo= +github.com/stretchr/testify v1.7.1/go.mod h1:6Fq8oRcR53rry900zMqJjRRixrwX3KX962/h/Wwjteg= +github.com/stretchr/testify v1.8.0/go.mod h1:yNjHg4UonilssWZ8iaSj1OCr/vHnekPRkoO+kdMU+MU= +github.com/stretchr/testify v1.8.4 h1:CcVxjf3Q8PM0mHUKJCdn+eZZtm5yQwehR5yeSVQQcUk= +github.com/stretchr/testify v1.8.4/go.mod h1:sz/lmYIOXD/1dqDmKjjqLyZ2RngseejIcXlSw2iwfAo= gopkg.in/check.v1 v0.0.0-20161208181325-20d25e280405 h1:yhCVgyC4o1eVCa2tZl7eS0r+SDo693bJlVdllGtEeKM= gopkg.in/check.v1 v0.0.0-20161208181325-20d25e280405/go.mod h1:Co6ibVJAznAaIkqp8huTwlJQCZ016jof/cbN4VW5Yz0= -gopkg.in/yaml.v3 v3.0.0-20200313102051-9f266ea9e77c h1:dUUwHk2QECo/6vqA44rthZ8ie2QXMNeKRTHCNY2nXvo= gopkg.in/yaml.v3 v3.0.0-20200313102051-9f266ea9e77c/go.mod h1:K4uyk7z7BCEPqu6E+C64Yfv1cQ7kz7rIZviUmN+EgEM= +gopkg.in/yaml.v3 v3.0.1 h1:fxVm/GzAzEWqLHuvctI91KS9hhNmmWOoWu0XTYJS7CA= +gopkg.in/yaml.v3 v3.0.1/go.mod h1:K4uyk7z7BCEPqu6E+C64Yfv1cQ7kz7rIZviUmN+EgEM= diff --git a/readme/readme_test.go b/readme/readme_test.go index 36b08c1..f0a11f4 100644 --- a/readme/readme_test.go +++ b/readme/readme_test.go @@ -7,7 +7,6 @@ import ( "github.com/pymaxion/geographiclib-go/geodesic" "github.com/pymaxion/geographiclib-go/geodesic/capabilities" - "github.com/stretchr/testify/assert" ) From cee74cb82fe23ed1bc7495d4294a7cdf84600b5c Mon Sep 17 00:00:00 2001 From: Patrick Yukman Date: Sat, 11 Nov 2023 14:45:56 -0800 Subject: [PATCH 2/6] implement gnomonic forward and reverse calculations, plus basic tests --- .gitignore | 3 ++ geodesic/gnomonic.go | 101 ++++++++++++++++++++++++++++++++++++++ geodesic/gnomonic_test.go | 29 +++++++++++ geodesic/inverse.go | 2 +- 4 files changed, 134 insertions(+), 1 deletion(-) create mode 100644 geodesic/gnomonic.go create mode 100644 geodesic/gnomonic_test.go diff --git a/.gitignore b/.gitignore index 66fd13c..fffa119 100644 --- a/.gitignore +++ b/.gitignore @@ -13,3 +13,6 @@ # Dependency directories (remove the comment below to include it) # vendor/ + +# Custom ignores +data diff --git a/geodesic/gnomonic.go b/geodesic/gnomonic.go new file mode 100644 index 0000000..a7bc2a8 --- /dev/null +++ b/geodesic/gnomonic.go @@ -0,0 +1,101 @@ +package geodesic + +import ( + "math" + + "github.com/pymaxion/geographiclib-go/geodesic/capabilities" +) + +type Gnomonic struct { + Earth *Geodesic +} + +func (g *Gnomonic) Forward(lat0, lon0, lat, lon float64) GnomonicData { + fwd := newGnomonicData() + fwd.Lat0, fwd.Lon0, fwd.Lat, fwd.Lon = lat0, lon0, lat, lon + + caps := capabilities.Azimuth | capabilities.GeodesicScale | capabilities.ReducedLength + inv := g.Earth.InverseWithCapabilities(lat0, lon0, lat, lon, caps) + fwd.Azi, fwd.Rk = inv.Azi2, inv.M12 + + if inv.M12 > 0 { + rho := inv.M12Reduced / inv.M12 + x, y := sincosd(inv.Azi1) + fwd.X, fwd.Y = rho*x, rho*y + } + return fwd +} + +func (g *Gnomonic) Reverse(lat0, lon0, x, y float64) GnomonicData { + rev := newGnomonicData() + rev.Lat0, rev.Lon0, rev.X, rev.Y = lat0, lon0, x, y + + azi0 := atan2d(x, y) + rho := math.Hypot(x, y) + a := g.Earth.EquatorialRadius() + s := a * math.Atan(rho/a) + + little := rho <= a + if !little { + rho = 1 / rho + } + + caps := capabilities.Latitude | capabilities.Longitude | capabilities.Azimuth | capabilities.DistanceIn | + capabilities.ReducedLength | capabilities.GeodesicScale + line := g.Earth.LineWithCapabilities(lat0, lon0, azi0, caps) + + var pos Data + numIt := 10 + trip := 0 + tripEpsilon := 0.01 * math.Sqrt(epsilon) + + for count := numIt; count > 0; count-- { + pos = line.PositionWithCapabilities(s, caps) + if trip > 0 { + break + } + + var ds float64 + if little { + ds = ((pos.M12Reduced / pos.M12) - rho) * pos.M12 * pos.M12 + } else { + ds = (rho - (pos.M12 / pos.M12Reduced)) * pos.M12Reduced * pos.M12Reduced + } + s -= ds + + if math.Abs(ds) < tripEpsilon*a { + trip++ + } + } + + if trip == 0 { + return rev + } + + rev.Lat, rev.Lon, rev.Azi, rev.Rk = pos.Lat2, pos.Lon2, pos.Azi1, pos.M12 + return rev +} + +type GnomonicData struct { + Lat0 float64 + Lon0 float64 + Lat float64 + Lon float64 + X float64 + Y float64 + Azi float64 + Rk float64 +} + +func newGnomonicData() GnomonicData { + return GnomonicData{ + Lat0: math.NaN(), + Lon0: math.NaN(), + Lat: math.NaN(), + Lon: math.NaN(), + X: math.NaN(), + Y: math.NaN(), + Azi: math.NaN(), + Rk: math.NaN(), + } +} diff --git a/geodesic/gnomonic_test.go b/geodesic/gnomonic_test.go new file mode 100644 index 0000000..d3b356a --- /dev/null +++ b/geodesic/gnomonic_test.go @@ -0,0 +1,29 @@ +package geodesic + +import ( + "testing" + + "github.com/stretchr/testify/assert" +) + +func TestGnomonicForward(t *testing.T) { + t.Run("forward calc", func(t *testing.T) { + lat0, lon0 := 48+50/60.0, 2+20/60.0 // Paris + lat, lon := 50.9, 1.8 // Calais + g := Gnomonic{Earth: WGS84} + r := g.Forward(lat0, lon0, lat, lon) + assert.InDelta(t, -37543.7, r.X, 0.05) + assert.InDelta(t, 230103, r.Y, 0.25) + }) +} + +func TestGnomonicReverse(t *testing.T) { + t.Run("reverse calc", func(t *testing.T) { + lat0, lon0 := 48+50/60.0, 2+20/60.0 // Paris + x, y := -38e3, 230e3 // Calais + g := Gnomonic{Earth: WGS84} + r := g.Reverse(lat0, lon0, x, y) + assert.InDelta(t, 50.899, r.Lat, 0.0005) + assert.InDelta(t, 1.79353, r.Lon, 0.000005) + }) +} diff --git a/geodesic/inverse.go b/geodesic/inverse.go index ef41359..0d702a4 100644 --- a/geodesic/inverse.go +++ b/geodesic/inverse.go @@ -99,7 +99,7 @@ func (s *inverseSolver) genInverse(lat1, lon1, lat2, lon2 float64, caps capabili // |bet2|. Alternatively (cbet1 >= -sbet1), abs(sbet2) + sbet1 is a better // measure. This logic is used in assigning calp2 in Lambda12. Sometimes these // quantities vanish and in that case we force bet2 = +/- bet1 exactly. An - // example where is is necessary is the inverse problem 48.522876735459 0 + // example where this is necessary is the inverse problem 48.522876735459 0 // -48.52287673545898293 179.599720456223079643 which failed with Visual Studio // 10 (Release and Debug) if cbet1 < -sbet1 { From ac5ff0a13439b4188e87e662dd7152e130c7c584 Mon Sep 17 00:00:00 2001 From: Patrick Yukman Date: Sat, 11 Nov 2023 15:10:06 -0800 Subject: [PATCH 3/6] add workflows for build and tests --- .github/workflows/build.yml | 24 ++++++++++++++++++++++++ .github/workflows/tests.yml | 26 ++++++++++++++++++++++++++ README.md | 2 ++ 3 files changed, 52 insertions(+) create mode 100644 .github/workflows/build.yml create mode 100644 .github/workflows/tests.yml diff --git a/.github/workflows/build.yml b/.github/workflows/build.yml new file mode 100644 index 0000000..22dd1aa --- /dev/null +++ b/.github/workflows/build.yml @@ -0,0 +1,24 @@ +# This workflow will build a golang project +# For more information see: https://docs.github.com/en/actions/automating-builds-and-tests/building-and-testing-go + +name: Build + +on: + push: + pull_request: + +jobs: + build: + runs-on: ubuntu-latest + steps: + - uses: actions/checkout@v4 + - name: Set up Go + uses: actions/setup-go@v4 + with: + go-version: '1.17' + - name: Display Go version + run: go version + - name: Install dependencies + run: go get ./... + - name: Build + run: go build -v ./... diff --git a/.github/workflows/tests.yml b/.github/workflows/tests.yml new file mode 100644 index 0000000..6471524 --- /dev/null +++ b/.github/workflows/tests.yml @@ -0,0 +1,26 @@ +# This workflow will build a golang project +# For more information see: https://docs.github.com/en/actions/automating-builds-and-tests/building-and-testing-go + +name: Tests + +on: + push: + pull_request: + +jobs: + build: + runs-on: ubuntu-latest + steps: + - uses: actions/checkout@v4 + - name: Set up Go + uses: actions/setup-go@v4 + with: + go-version: '1.17' + - name: Display Go version + run: go version + - name: Install dependencies + run: go get ./... + - name: Build + run: go build -v ./... + - name: Tests + run: go test -v ./... diff --git a/README.md b/README.md index 6c39b7d..ebbf362 100644 --- a/README.md +++ b/README.md @@ -3,6 +3,8 @@ [![Go Reference](https://pkg.go.dev/badge/github.com/pymaxion/geographiclib-go.svg)](https://pkg.go.dev/github.com/pymaxion/geographiclib-go) [![GitHub go.mod Go version of a Go module](https://img.shields.io/github/go-mod/go-version/pymaxion/geographiclib-go.svg)](https://github.com/pymaxion/geographiclib-go) [![Go Report Card](https://goreportcard.com/badge/github.com/pymaxion/geographiclib-go)](https://goreportcard.com/report/github.com/pymaxion/geographiclib-go) +![Build](https://github.com/pymaxion/geographiclib-go/actions/workflows/build.yml/badge.svg) ![Tests](https://github.com/pymaxion/geographiclib-go/actions/workflows/tests.yml/badge.svg) + # Description This is a Go implementation of the geodesic algorithms from Charles F. F. Karney's [GeographicLib](https://geographiclib.sourceforge.io/). Though not an official implementation of GeographicLib, geographiclib-go has feature parity with the officially-maintained [Java](https://geographiclib.sourceforge.io/html/java/) and [Python](https://geographiclib.sourceforge.io/html/python/) implementations and additionally includes a utility to validate the implementation against the official 500k-line [GeodTest.dat](https://geographiclib.sourceforge.io/html/geodesic.html#testgeod) repository of geodesic test data. From 480c3eee2c04e33b42f8437f9ad48816ad464054 Mon Sep 17 00:00:00 2001 From: Patrick Yukman Date: Mon, 20 Nov 2023 12:04:38 -0800 Subject: [PATCH 4/6] use native github latex rendering for readme --- README.md | 74 +++++++++++++++++++++++++++---------------------------- 1 file changed, 37 insertions(+), 37 deletions(-) diff --git a/README.md b/README.md index ebbf362..19a8e05 100644 --- a/README.md +++ b/README.md @@ -36,9 +36,9 @@ See the `examples` directory for some basic code samples. # Geodesics on an ellipsoid ## Introduction -Consider an ellipsoid of revolution with equatorial radius , polar semi-axis , and flattening . Points on the surface of the ellipsoid are characterized by their latitude and longitude . (Note that latitude here means the *geographical latitude*, the angle between the normal to the ellipsoid and the equatorial plane). +Consider an ellipsoid of revolution with equatorial radius $a$, polar semi-axis $b$, and flattening $f=\frac{a-b}{a}$. Points on the surface of the ellipsoid are characterized by their latitude $\varphi$ and longitude $\lambda$. (Note that latitude here means the *geographical latitude*, the angle between the normal to the ellipsoid and the equatorial plane). -The shortest path between two points on the ellipsoid at and is called the *geodesic*. Its length is and the geodesic from point 1 to point 2 has forward azimuths and at the two end points. In this figure, we have . +The shortest path between two points on the ellipsoid at $(\varphi_1, \lambda_1)$ and $(\varphi_2, \lambda_2)$ is called the *geodesic*. Its length is $s_{12}$ and the geodesic from point 1 to point 2 has forward azimuths $\alpha_1$ and $\alpha_2$ at the two end points. In this figure, we have $\lambda_{12}=\lambda_2-\lambda_1$. diagram @@ -46,33 +46,33 @@ A geodesic can be extended indefinitely by requiring that any sufficiently small Traditionally two geodesic problems are considered: -* the direct problem — given , , , , determine , , and ; this is solved by `Geodesic.Direct()`. -* the inverse problem — given , , , , determine , , and ; this is solved by `Geodesic.Inverse()`. +* the direct problem — given $\varphi_1$, $\lambda_1$, $\alpha_1$, $s_{12}$, determine $\varphi_2$, $\lambda_2$, and $\alpha_2$; this is solved by `Geodesic.Direct()`. +* the inverse problem — given $\varphi_1$, $\lambda_1$, $\varphi_2$, $\lambda_2$, determine $s_{12}$, $\alpha_1$, and $\alpha_2$; this is solved by `Geodesic.Inverse()`. ## Additional properties The functions which solve the geodesic problems also calculate several other quantities of interest: -* is the area between the geodesic from point 1 to point 2 and the equator; i.e., it is the area, measured counter-clockwise, of the quadrilateral with corners , , , and . It is given in meters squared. -* , the reduced length of the geodesic, is defined such that if the initial azimuth is perturbed by (radians) then the second point is displaced by in the direction perpendicular to the geodesic. is given in meters. On a curved surface the reduced length obeys a symmetry relation, . On a flat surface, we have . -* and are geodesic scales. If two geodesics are parallel at point 1 and separated by a small distance , then they are separated by a distance at point 2. is defined similarly (with the geodesics being parallel to one another at point 2). and are dimensionless quantities. On a flat surface, we have . -* is the arc length on the auxiliary sphere. This is a construct for converting the problem to one in spherical trigonometry. The spherical arc length from one equator crossing to the next is always 180°. +* $S_{12}$ is the area between the geodesic from point 1 to point 2 and the equator; i.e., it is the area, measured counter-clockwise, of the quadrilateral with corners $(\varphi_1, \lambda_1)$, $(0,\lambda_1)$, $(0,\lambda_2)$, and $(\varphi_2, \lambda_2)$. It is given in meters squared. +* $m_{12}$, the reduced length of the geodesic, is defined such that if the initial azimuth is perturbed by $d\alpha_1$ (radians) then the second point is displaced by $m_{12}d\alpha_1$ in the direction perpendicular to the geodesic. $m_{12}$ is given in meters. On a curved surface the reduced length obeys a symmetry relation, $m_{12}+m_{21}=0$. On a flat surface, we have $m_{12}=s_{12}$. +* $M_{12}$ and $M_{21}$ are geodesic scales. If two geodesics are parallel at point 1 and separated by a small distance $dt$, then they are separated by a distance $M_{12}dt$ at point 2. $M_{21}$ is defined similarly (with the geodesics being parallel to one another at point 2). $M_{12}$ and $M_{21}$ are dimensionless quantities. On a flat surface, we have $M_{12}=M_{21}=1$. +* $\sigma_{12}$ is the arc length on the auxiliary sphere. This is a construct for converting the problem to one in spherical trigonometry. The spherical arc length from one equator crossing to the next is always 180°. If points 1, 2, and 3 lie on a single geodesic, then the following addition rules hold: -* -* -* -* -* -* +* $s_{13}=s_{12}+s_{23}$ +* $\sigma_{13}=\sigma_{12}+\sigma_{23}$ +* $S_{13}=S_{12}+S_{23}$ +* $m_{13}=m_{12}M_{23}+m_{23}M_{21}$ +* $M_{13}=M_{12}M_{23}-\frac{(1-M_{12}M_{21})m_{23}}{m_{12}}$ +* $M_{31}=M_{32}M_{21}-\frac{(1-M_{23}M_{32})m_{12}}{m_{23}}$ ## Multiple shortest geodesics The shortest distance found by solving the inverse problem is (obviously) uniquely defined. However, in a few special cases there are multiple azimuths which yield the same shortest distance. Here is a catalog of those cases: -* (with neither point at a pole). If , the geodesic is unique. Otherwise there are two geodesics and the second one is obtained by setting , , . (This occurs when the longitude difference is near for oblate ellipsoids.) -* (with neither point at a pole). If or , the geodesic is unique. Otherwise there are two geodesics and the second one is obtained by setting , . (This occurs when is near for prolate ellipsoids.) -* Points 1 and 2 at opposite poles. There are infinitely many geodesics which can be generated by setting , for arbitrary . (For spheres, this prescription applies when points 1 and 2 are antipodal.) -* (coincident points). There are infinitely many geodesics which can be generated by setting , for arbitrary . +* $\varphi_1=-\varphi_2$ (with neither point at a pole). If $\alpha_1=\alpha2$, the geodesic is unique. Otherwise there are two geodesics and the second one is obtained by setting $[\alpha_1,\alpha_2]\leftarrow[\alpha_2,\alpha_1]$, $[M_{12},M_{21}]\leftarrow[M_{21},M_{12}]$, $S_{12}\leftarrow-S_{12}$. (This occurs when the longitude difference is near $\pm180^{\circ}$ for oblate ellipsoids.) +* $\lambda_2=\lambda_1\pm180^{\circ}$ (with neither point at a pole). If $\alpha_1=0^{\circ}$ or $\pm180^{\circ}$, the geodesic is unique. Otherwise there are two geodesics and the second one is obtained by setting $[\alpha_1,\alpha_2]\leftarrow[-\alpha_1,-\alpha_2]$, $S_{12}\leftarrow-S_{12}$. (This occurs when $\varphi_2$ is near $-\varphi_1$ for prolate ellipsoids.) +* Points 1 and 2 at opposite poles. There are infinitely many geodesics which can be generated by setting $[\alpha_1,\alpha_2]\leftarrow[\alpha_1,\alpha_2]+[\delta,-\delta]$, for arbitrary $\delta$. (For spheres, this prescription applies when points 1 and 2 are antipodal.) +* $s_{12}=0$ (coincident points). There are infinitely many geodesics which can be generated by setting $[\alpha_1,\alpha_2]\leftarrow[\alpha_1,\alpha_2]+[\delta,\delta]$, for arbitrary $\delta$. ## Background The algorithms implemented by this package are given in Karney (2013) and are based on Bessel (1825) and Helmert (1880); the algorithm for areas is based on Danielsen (1989). These improve on the work of Vincenty (1975) in the following respects: @@ -83,25 +83,25 @@ The algorithms implemented by this package are given in Karney (2013) and are ba # The library interface ## Units -All angles (latitude, longitude, azimuth, arc length) are measured in degrees with latitudes increasing northwards, longitudes increasing eastwards, and azimuths measured clockwise from north. For a point at a pole, the azimuth is defined by keeping the longitude fixed, writing , and taking the limit . +All angles (latitude, longitude, azimuth, arc length) are measured in degrees with latitudes increasing northwards, longitudes increasing eastwards, and azimuths measured clockwise from north. For a point at a pole, the azimuth is defined by keeping the longitude fixed, writing $\varphi=\pm(90^{\circ}-\epsilon)$, and taking the limit $\epsilon\rightarrow0+$. ## `geodesic.Data` The results returned by `Geodesic.Direct()`, `Geodesic.Inverse()`, `Line.Position()`, etc., return a struct type (`geodesic.Data`) with some of the following 12 fields set: -* `Lat1` = , latitude of point 1 (degrees) -* `Lon1` = , longitude of point 1 (degrees) -* `Azi1` = , azimuth of line at point 1 (degrees) -* `Lat2` = , latitude of point 2 (degrees) -* `Lon2` = , longitude of point 2 (degrees) -* `Azi2` = , (forward) azimuth of line at point 2 (degrees) -* `S12` = , distance from 1 to 2 (meters) -* `A12` = , arc length on auxiliary sphere from 1 to 2 (degrees) -* `M12Reduced` = , reduced length of geodesic (meters) -* `M12` = , geodesic scale at 2 relative to 1 (dimensionless) -* `M21` = , geodesic scale at 1 relative to 2 (dimensionless) -* `S12Area` = , area between geodesic and equator (meters) - -**Note**: because of Go's use of capitalization to denote exported symbols, the names of a few of these fields do not exactly match their symbolic representations. For example, it is not possible to distinguish between and using lowercase vs. capitalized field names as is possible in other languages' implementations of GeographicLib; *must* be capitalized to be accessible from other packages. Thus, where necessary to remove ambiguity, suffixes have been added to the less-commonly-used field names (e.g., `S12`, `S12Area`). +* `Lat1` = $\varphi_1$, latitude of point 1 (degrees) +* `Lon1` = $\lambda_1$, longitude of point 1 (degrees) +* `Azi1` = $\alpha_1$, azimuth of line at point 1 (degrees) +* `Lat2` = $\varphi_2$, latitude of point 2 (degrees) +* `Lon2` = $\lambda_2$, longitude of point 2 (degrees) +* `Azi2` = $\alpha_2$, (forward) azimuth of line at point 2 (degrees) +* `S12` = $s_{12}$, distance from 1 to 2 (meters) +* `A12` = $\sigma_{12}$, arc length on auxiliary sphere from 1 to 2 (degrees) +* `M12Reduced` = $m_{12}$, reduced length of geodesic (meters) +* `M12` = $M_{12}$, geodesic scale at 2 relative to 1 (dimensionless) +* `M21` = $M_{21}$, geodesic scale at 1 relative to 2 (dimensionless) +* `S12Area` = $S_{12}$, area between geodesic and equator (meters$^2$) + +**Note**: because of Go's use of capitalization to denote exported symbols, the names of a few of these fields do not exactly match their symbolic representations. For example, it is not possible to distinguish between $s_{12}$ and $S_{12}$ using lowercase vs. capitalized field names as is possible in other languages' implementations of GeographicLib; $s_{12}$ *must* be capitalized to be accessible from other packages. Thus, where necessary to remove ambiguity, suffixes have been added to the less-commonly-used field names (e.g., $s_{12}\rightarrow$`S12`, $S_{12}\rightarrow$`S12Area`). ## `capabilities.Mask` The standard geodesic functions (`Geodesic.Direct()`, `Geodesic.Inverse()`, etc.) all return the 7 basic quantities: `Lat1`, `Lon1`, `Azi1`, `Lat2`, `Lon2`, `Azi2`, `S12`, together with the arc length `A12`. However, each function has a counterpart `WithCapabilities()` requiring an additional parameter of type `capabilities.Mask`; this additional parameter is used to tailor which quantities to calculate. In addition, when a `Line` is instantiated, it too can be provided with a `capabilities.Mask` parameter specifying what quantities can be returned from the resulting instance. @@ -132,10 +132,10 @@ A custom `capabilities.Mask` parameter can be created by bitwise-or(`|`)’ing t * The distance `S12` is unrestricted. This allows geodesics to wrap around the ellipsoid. Such geodesics are no longer shortest paths; however, they retain the property that they are the straightest curves on the surface. * Similarly, the spherical arc length `A12` is unrestricted. * Longitudes and azimuths are unrestricted; internally these are exactly reduced to the range [-180°, 180°), but see also the `LongUnroll` capability. -* The equatorial radius and the polar semi-axis must both be positive and finite (this implies that ). -* The flattening should satisfy in order to retain full accuracy. This condition holds for most applications in geodesy. +* The equatorial radius $a$ and the polar semi-axis $b$ must both be positive and finite (this implies that $-\infty \lt f \lt 1$). +* The flattening $f$ should satisfy $f\in[-\frac{1}{50},\frac{1}{50}]$ in order to retain full accuracy. This condition holds for most applications in geodesy. -Reasonably accurate results can be obtained for . Here is a table of the approximate maximum error (expressed as a distance) for an ellipsoid with the same equatorial radius as the WGS84 ellipsoid and different values of the flattening. +Reasonably accurate results can be obtained for $-0.2\le f\le0.2$. Here is a table of the approximate maximum error (expressed as a distance) for an ellipsoid with the same equatorial radius as the WGS84 ellipsoid and different values of the flattening. |abs(f)|error| |---|---| @@ -146,7 +146,7 @@ Reasonably accurate results can be obtained for m (not 1 nautical mile!) +Here 1 nm = 1 nanometer = $10^{-9}$ m (not 1 nautical mile!) # Examples ## Initializing From 5db89983287023a832d3245bbd0f69ea9dce62b6 Mon Sep 17 00:00:00 2001 From: Patrick Yukman Date: Mon, 20 Nov 2023 12:36:09 -0800 Subject: [PATCH 5/6] add documentatiuon to gnomonic.go --- geodesic/gnomonic.go | 101 +++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 101 insertions(+) diff --git a/geodesic/gnomonic.go b/geodesic/gnomonic.go index a7bc2a8..cd3576c 100644 --- a/geodesic/gnomonic.go +++ b/geodesic/gnomonic.go @@ -6,10 +6,89 @@ import ( "github.com/pymaxion/geographiclib-go/geodesic/capabilities" ) +/* + Gnomonic is a type representing calculations for gnomonic projections centered + at an arbitrary position C on the ellipsoid. This projection is derived in + Section 8 of C. F. F. Karney, Algorithms for geodesics, J. Geodesy 87, 43–55 + (2013); DOI: 10.1007/s00190-012-0578-z; addenda: + https://geographiclib.sourceforge.io/geod-addenda.html. + + The gnomonic projection of a point P on the ellipsoid is defined as follows: + compute the geodesic line from C to P; compute the reduced length m12, geodesic + scale M12, and ρ = m12/M12; finally, this gives the coordinates x and y of P in + gnomonic projection with x = ρ sin azi1; y = ρ cos azi1, where azi1 is the + azimuth of the geodesic at C. The method Forward(double, double, double, double) + performs the forward projection and Reverse(double, double, double, double) is + the inverse of the projection. The methods also return the azimuth azi of the + geodesic at P and reciprocal scale rk in the azimuthal direction. The scale in + the radial direction is 1/(rk^2). + + For a sphere, ρ reduces to a tan(s12/a), where s12 is the length of the geodesic + from C to P, and the gnomonic projection has the property that all geodesics + appear as straight lines. For an ellipsoid, this property holds only for + geodesics intersecting the centers. However geodesic segments close to the center + are approximately straight. + + Consider a geodesic segment of length l. Let T be the point on the geodesic + (extended if necessary) closest to C, the center of the projection, and t, be + the distance CT. To lowest order, the maximum deviation (as a true distance) of + the corresponding gnomonic line segment (i.e., with the same end points) from + the geodesic is: + + K(T) − K(C)) l^2 t / 32. + + where K is the Gaussian curvature. + + This result applies for any surface. For an ellipsoid of revolution, consider + all geodesics whose end points are within a distance r of C. For a given r, the + deviation is maximum when the latitude of C is 45°, when endpoints are a + distance r away, and when their azimuths from the center are ± 45° or ± 135°. To + lowest order in r and the flattening f, the deviation is f (r/2a)^3 r. + + CAUTION: The definition of this projection for a sphere is standard. However, + there is no standard for how it should be extended to an ellipsoid. The choices + are: + * Declare that the projection is undefined for an ellipsoid. + * Project to a tangent plane from the center of the ellipsoid. This causes great + ellipses to appear as straight lines in the projection; i.e., it generalizes the + spherical great circle to a great ellipse. This was proposed by independently by + Bowring and Williams in 1997. + * Project to the conformal sphere with the constant of integration chosen so + that the values of the latitude match for the center point and perform a central + projection onto the plane tangent to the conformal sphere at the center point. + This causes normal sections through the center point to appear as straight lines + in the projection; i.e., it generalizes the spherical great circle to a normal + section. This was proposed by I. G. Letoval'tsev, Generalization of the gnomonic + projection for a spheroid and the principal geodetic problems involved in the + alignment of surface routes, Geodesy and Aerophotography (5), 271–274 (1963). + * The projection given here. This causes geodesics close to the center point to + appear as straight lines in the projection; i.e., it generalizes the spherical + great circle to a geodesic. + + The algorithms are described in + C. F. F. Karney, Algorithms for geodesics, J. Geodesy 87, 43-55 (2013) + Link: https://doi.org/10.1007/s00190-012-0578-z + Addenda: https://geographiclib.sourceforge.io/geod-addenda.html +*/ type Gnomonic struct { Earth *Geodesic } +/* + Forward performs the forward projection from geographic to gnomonic. + + lat0: latitude of center point of projection (degrees). Should be in the range [-90°, 90°]. + lon0: longitude of center point of projection (degrees). Should be in the range [−540°, 540°). + lat - latitude of point (degrees). Should be in the range [-90°, 90°]. + lon - longitude of point (degrees). Should be in the range [−540°, 540°). + + The scale of the projection is 1/rk2 in the "radial" direction, azi clockwise + from true north, and is 1/rk in the direction perpendicular to this. If the + point lies "over the horizon", i.e., if rk ≤ 0, then NaNs are returned for x and + y (the correct values are returned for azi and rk). A call to Forward followed + by a call to Reverse will return the original (lat, lon) (to within roundoff) + provided the point is not over the horizon. +*/ func (g *Gnomonic) Forward(lat0, lon0, lat, lon float64) GnomonicData { fwd := newGnomonicData() fwd.Lat0, fwd.Lon0, fwd.Lat, fwd.Lon = lat0, lon0, lat, lon @@ -26,6 +105,23 @@ func (g *Gnomonic) Forward(lat0, lon0, lat, lon float64) GnomonicData { return fwd } +/* + Reverse performs the reverse projection from gnomonic to geographic. + + lat0: latitude of center point of projection (degrees). Should be in the range [-90°, 90°]. + lon0: longitude of center point of projection (degrees). Should be in the range [−540°, 540°). + x: easting of point (meters). + y: northing of point (meters). + + The returned values GnomonicData.Lat and GnomonicData.Lon are in the range [-180°, 180°]. + + The scale of the projection is 1/rk2 in the "radial" direction, azi clockwise + from true north, and is 1/rk in the direction perpendicular to this. Even though + all inputs should return a valid lat and lon, it's possible that the procedure + fails to converge for very large x or y; in this case NaNs are returned for all + the output arguments. A call to Reverse followed by a call to Forward will + return the original (x, y) (to roundoff). +*/ func (g *Gnomonic) Reverse(lat0, lon0, x, y float64) GnomonicData { rev := newGnomonicData() rev.Lat0, rev.Lon0, rev.X, rev.Y = lat0, lon0, x, y @@ -76,6 +172,11 @@ func (g *Gnomonic) Reverse(lat0, lon0, x, y float64) GnomonicData { return rev } +// GnomonicData describes the results of gnomonic projection. This is used to +// return the results for a gnomonic projection of a point (lat, lon) given a +// center point of projection (lat0, lon0). The returned GnomonicData structs +// always include the parameters provided to Gnomonic.Forward and +// Gnomonic.Reverse and it always includes the fields x, y, azi. and rk. type GnomonicData struct { Lat0 float64 Lon0 float64 From e1162bf3eb0cbf67ed200c0be56ece9133e758af Mon Sep 17 00:00:00 2001 From: Patrick Yukman Date: Mon, 20 Nov 2023 12:55:24 -0800 Subject: [PATCH 6/6] add version history to README --- README.md | 13 ++++++++++++- 1 file changed, 12 insertions(+), 1 deletion(-) diff --git a/README.md b/README.md index 19a8e05..b5771d5 100644 --- a/README.md +++ b/README.md @@ -9,9 +9,10 @@ This is a Go implementation of the geodesic algorithms from Charles F. F. Karney's [GeographicLib](https://geographiclib.sourceforge.io/). Though not an official implementation of GeographicLib, geographiclib-go has feature parity with the officially-maintained [Java](https://geographiclib.sourceforge.io/html/java/) and [Python](https://geographiclib.sourceforge.io/html/python/) implementations and additionally includes a utility to validate the implementation against the official 500k-line [GeodTest.dat](https://geographiclib.sourceforge.io/html/geodesic.html#testgeod) repository of geodesic test data. More information about the GeographicLib project can be found at https://geographiclib.sourceforge.io. This README lifts heavily from GeographicLib's excellent project documentation, especially the section [Geodesics on an ellipsoid](#geodesics-on-an-ellipsoid), which copies Karney's documentation verbatim to introduce the direct and inverse geodesic problems. - + # Contents 1. [Getting started](#getting-started) + 1. [Version history](#version-history) 1. [Geodesics on an ellipsoid](#geodesics-on-an-ellipsoid) 1. [Introduction](#introduction) 1. [Additional properties](#additional-properties) @@ -34,6 +35,16 @@ More information about the GeographicLib project can be found at https://geograp See the `examples` directory for some basic code samples. +## Version history +See [Releases](https://github.com/pymaxion/geographiclib-go/releases) for full details. + +| Version | Description | +|----------------------------------------------------------------------------|----------------------------------------| +| [v1.1.0](https://github.com/pymaxion/geographiclib-go/releases/tag/v1.1.0) | Add gnomonic projection calculations | +| [v1.0.1](https://github.com/pymaxion/geographiclib-go/releases/tag/v1.0.1) | Fix failing build on ARM architectures | +| [v1.0.0](https://github.com/pymaxion/geographiclib-go/releases/tag/v1.0.0) | Initial release | + + # Geodesics on an ellipsoid ## Introduction Consider an ellipsoid of revolution with equatorial radius $a$, polar semi-axis $b$, and flattening $f=\frac{a-b}{a}$. Points on the surface of the ellipsoid are characterized by their latitude $\varphi$ and longitude $\lambda$. (Note that latitude here means the *geographical latitude*, the angle between the normal to the ellipsoid and the equatorial plane).