juliaearth / coordrefsystems.jl Goto Github PK

CoordRefSystems currently does not implement a few coordinate conversions that I think it should.

julia> convert(Cartesian, WebMercator(0,0))
ERROR: MethodError: Cannot `convert` an object of type
  WebMercator{WGS84Latest, Quantity{Float64, 𝐋, Unitful.FreeUnits{(m,), 𝐋, nothing}}} to an object of type
  Cartesian

julia> convert(LatLon, LatLonAlt(0,0,0))
ERROR: MethodError: Cannot `convert` an object of type
  GeodeticLatLonAlt{WGS84Latest, Quantity{Float64, NoDims, Unitful.FreeUnits{(°,), NoDims, nothing}}, Quantity{Float64, 𝐋, Unitful.FreeUnits{(m,), 𝐋, nothing}}} to an object of type
  GeodeticLatLon

It would be great if these could be added. Thanks!

It would be great if CoordRefSystems could provide functions cstrip() and cwrap() for converting between coordinate types (e.g. x::LatLon) and coordinate vectors (e.g. SVector(x.lat, xlon)). Such functions are very useful for interfacing with non-CoordRefSystems-aware software (e.g. plotting libraries) and implementing coordinate arithmetic. For example:

function rhumb_line(x::CRS, y::CRS)
    C = Mercator
    return t -> cwrap(C, t * cstrip(C, x) + (1-t) * cstrip(C, y))
end

function mean(coords::AbstractVector{<:CRS})
    C = Cartesian{WGS84latest}
    m = cwrap(C, mean(cstrip(C, x) for x in coords))
    # Don't trust means that are too far underground
    # (i.e. the `coords` are not local enough)
    @assert abs(convert(LatLonAlt, m).alt) < 10
    return m
end

translated_latlon(x::CRS, d::SVector{3}) = cstrip(LatLon, cwrap(ENU(x), d))

I think it should be Proj.jl / Cartography.jl

It would be great if CoordRefSystems could implement the East North Up local coordinate system.

A major design issue when it comes to local coordinates is whether such coordinates should remember their origin. I would like to propose that in CoordRefSystems local coordinates should remember their origin but the cstrip() function (see #47) should return only the offset, i.e.

x = Cartesian(1,0,0)
y = Cartesian(1,1,0)
@assert cstrip(ENU(x), y) == SVector(1,0,0) # y is one unit to the east of x

This would allow local coordinates to represent a unique point in space just like their global counterparts, and at the same time it would allow the user to easily "relocate" a local coordinate using cwrap(ENU(y), cstrip(ENU(x), local_coord)).

In our design, datum are defined with a specific epoch. Sometimes, a dataset is provided years after the definition of the datum, and we need to be able to shift the epoch for dynamic corrections (i.e. Helmert transform).

Based on our investigations, the inverse of the TransverseMercator is not producing the same results of PROJ nor GeographicLib. This needs to be fixed in a second more careful round of review.

The issue shouldn't affect most use cases, but it does affect conversions between TransverseMercator (and UTM) and any other projected CRS since converting between two projected CRS requires calling the inverse of the source CRS.

The LatLon to Robinson conversion is not numerically stable.

julia> reldiff(x,y) = abs(x-y)/max(abs(x), abs(y));

julia> reldiff(
           convert(Robinson, LatLon(prevfloat(30.0,2), 40.0)).x, 
           convert(Robinson, LatLon(          30.0   , 40.0)).x,
       )
2.7474463318253135e-8

This can be an issue e.g. when doing "round trip" conversions.

julia> isapprox(
           convert(Robinson,                                    LatLon(30, 40)), 
           convert(Robinson, convert(LatLon, convert(Cartesian, LatLon(30, 40)))), 
       )
false

I've run into this while addressing #52 (comment). Not sure whether there are any real-world issues related to this, but I thought I'd flag it anyway just in case.

Hey! This looks like a really cool package, and native Julia projections are something I've wanted to introduce to GeoMakie for a while.

Is it possible to apply these transformations in a generic way, without converting between types?

That would simplify the interface considerably for me. For example, treating CRS types in a trait-like way with function signatures like convert(point, source_crs, dest_crs) or maybe a constructed Transformation object that can be called on any pointlike type would make allowing this in GeoMakie pretty easy and fast.

Also, have you considered tying in to https://github.com/JuliaGeo/GeoFormatTypes.jl? They are already integrated in th raster data space so it would be nice to be able to apply these types to other coordinates, even if it's just as an extension package to accept GeoFormatTypes and convert them to Cartography.jl types.

There are a few occasions where it would be handy to have local coordinates that "wrap" around the surface of the earth. For example:

• When putting axis ticks in meters / kilometers on a WebMercator map:

  yticks = [
      convert(
          WebMercator, 
          cwrap(
              WrappingEastNorth(map_origin), 
              (0u"km", y), 
          )
      ).y
     for y in (-5:5)u"km"
  ]

• We work with devices that communicate acoustically under water. In shallow bodies of water, sound is constrained to travel along the surface of the earth rather than tunnel through the seabed, so in this case working in a wrapping coordinate system is more accurate.

I must admit that I am not aware of any reference for such a coordinate system, so this proposal is quite speculative. Curious to hear what other people might think.

isapprox(::Cartesian, ::Cartesian) currently fails for points with lon ≈ ±90° (e.g. Memphis, Tennessee).

julia> r = convert(Cartesian, LatLon(0,0)).x
6.378137e6 m

julia> convert(Cartesian, LatLon(0, 90)) ≈ Cartesian{WGS84Latest}(0u"m", r, 0u"m")
false

The problem is that isapprox(::Cartesian, ::Cartesian) currently compares every entry one by one for approximate equality.

Instead, it should do a vector comparison, isapprox(x::Cartesian, y::Cartesian) = norm(x - y) / max(norm(x), norm(y)), analogous to how Base does it.

This issue is used to trigger TagBot; feel free to unsubscribe.

If you haven't already, you should update your TagBot.yml to include issue comment triggers.
Please see this post on Discourse for instructions and more details.

If you'd like for me to do this for you, comment TagBot fix on this issue.
I'll open a PR within a few hours, please be patient!