Status of geodesy packages in Julia (DGGS experiments)

Specific Domains Geo
,
1

Hi all,

I would like to instigate a short review discussion on geodesy-related packages in Julia. My background is GIS, geospatial, and recently Discrete Global Grid Systems (DGGS) in particular. If you have heard of Uber H3, Google S2, DGGRID, or A5 systems, this is the direction. I had some experiences with the Pangeo folks picking up DGGS [0]. But still it feels inconvenient in Python. The reason is, most DGGS-related software is currently implemented in C/C++ and increasingly Rust, with Python-bindings at different levels of maturity and speed. I am increasingly enticed by the “not 2 languages” promise of Julia and had some encouraging success in dabbling with JuliaGeo packages. So I would like to try to do more with DGGS in Julia.

One aspect in underlying geo logic is that most DGGS are implemented on the sphere and using them with ellpsoidal/WGS84 coordinates/geometries requires classic geodesic calculations to make an adjustment of the latitude values (authalic conversion) to be correct.
I am aware of two geodesy packages in Julia (Geodesy.jl and [1] and Geodesics.jl [2]). Geodesy.jl is apparently in maintenance mode, so I reached out to Geodesics.jl. I managed to port an extract of some of the routines for authalic latitude conversion to Julia ([3]), based from a Python package [4], which is based, like most packages, one way or another on Geographiclib by M. Karney [5].

But my question is also if there is interest in general to explore/compose native Julia implementations. There are DGGS-related activities like HEALPIX.jl [6] and DGGS.jl [7].

Thanks for consideration, happy to discuss more.
Alex

Some references:

2

Hi @allixender ,

Take a look into CoordRefSystems.jl and Meshes.jl on top of it. It is a strongly typed system with no equivalent in other languages.

Please reach out if you have questions or would like to contribute.

1 Like
3

Thanks, yes and no :slight_smile: Preferably in DGGS I’d like to go as long as possible without geometries, DGGS are also quite well-structured, so a mesh is a bit of a step backwards. But the package looks very powerful.
And CoordRefSystems look good, too, but if I see correclty it also does not yet implement spherical to ellpisoid. Would you see a way forward to add this (or use my draft implementation as provided above) ?

4

Packages like ArchGDAL and GMT.jl or others that depend on GDAL have access to the PROJ library, which AFAIK is what almost everybody uses in the outside world.

5

You might be surprised if you study the Meshes.jl source code. Dont be fooled by the name. Grids are structured and fully supported with lazy data structures.

It does. Study the source code in detail and feel free to ask questions. If it doesnt have a specific conversion yet, it should be straightforward to add. We already worked out the design to acommodate all use cases.

2 Likes
7

But I am not sure how to use it from Proj, I can’t find documentation that actually show to use authalic conversion, even in Pyproj.

8

There seems to be authalic to geodetic in the code, but not in the docs: # AuthalicLatLon <> GeodeticLatLon

Probably gonna have to try.

1 Like
9

If you want to go to a sphere, you could just specify the +R_A flag to Proj for your target CRS - that will make the target ellipsoid be a sphere with the same surface area as the reference ellipsoid, thus, an authalic transformation.

See the Ellipsoids page of the Proj docs for more information.

1 Like
10

SO, I tested both CoordRefSystems and Proj (“+proj=latlon +R_A” like this?), and both do not transform as I expect them to (reference confirmed by geographiclib, pygeodesy, and another independent implementation. But both CoordRefSystems and Proj implement the proj way and behave basically the same.

The geographic latitude 58.39714590 needs to match with the authalic latitude 58.28252559

This is how tested with CoordRefSystems.

using CoordRefSystems

latlon = LatLon{WGS84Latest}(lat, lon)

check = convert(AuthalicLatLon{WGS84Latest}, latlon )

round_trip_lat = convert(LatLon{WGS84Latest}, check )
11

From the GMT manual (the mapproject module has not yet received a GMT.jl translation)

using GMT

D = mapproject([0 58.39714590], N=:a)

1×2 GMTdataset{Float64, 2}
 Row │ col.1    col.2
─────┼────────────────
   1 │   0.0  58.2825

D.data[2]
58.2825255919146
1 Like
12

@allixender you can omit the datum since it is the default and you dont want to change it during the conversion:

using CoordRefSystems

latlon = LatLon(lat, lon)

check = convert(AuthalicLatLon, latlon )

round_trip_lat = convert(LatLon, check )

Or you could use an older realization of the WGS84 datum that is not the WGS84Latest alias.

CoordRefSystems.jl should match C PROJ to the last digit in double precision. Can you point in the source code what is causing the difference between our definition and yours? We should probably continue this discussion in a GitHub issue.

13

I cannot, because I am not even sure if I used it correctly. Thanks by the way @joa-quim , this looks like I would expect it to look like.

@juliohm If you could confirm the correct use below?

The geographic latitude 58.39714590 needs to match with the authalic latitude 58.28252559

julia> using CoordRefSystems

julia> latlon = LatLon{WGS84Latest}(58.39714590, 0)
GeodeticLatLon{WGS84Latest} coordinates
├─ lat: 58.3971459°
└─ lon: 0.0°

julia> check = convert(AuthalicLatLon{WGS84Latest}, latlon )
AuthalicLatLon{WGS84Latest} coordinates
├─ lat: 58.282525581092806°
└─ lon: 0.0°

julia> round_trip_lat = convert(LatLon{WGS84Latest}, check )
GeodeticLatLon{WGS84Latest} coordinates
├─ lat: 58.39714589834038°
└─ lon: 0.0°

julia> isapprox(round_trip_lat, latlon)
true

Ok, weird (in a good sense), this seems to work now. I wonder what I did wrong?

1 Like
14

Hi @allixender, I would just omit the datum as explained above to make the code more readable. Other than that, I am glad it is producing the result you need.

Regarding the more general discussion on DGGS, try to study the Grid subtypes and the different Topology subtypes in Meshes.jl. They might shed some light on how to develop new grids that are compatible with our stack.

1 Like