JuliaGeo / GeoInterface.jl

A Julia Protocol for Geospatial Data
https://juliageo.org/GeoInterface.jl/stable/
MIT License
95 stars 17 forks source link

Add CRS traits #119

Closed evetion closed 2 months ago

evetion commented 8 months ago

This is a design issue, feel free to comment and/or come up with your own thoughts.

I would like to be able to distinguish between spherical and cartesian geometries and coordinate systems. For example, if you currently calculate (using Geomorphology) the slope on a geographic raster the answer is wrong. If you currently calculate the area of a geographic geometry, it is wrong. Libraries like LibGEOS and others assume cartesian coordinates, but don't (and can't) error or warn for geographic geometries. This is something we extensively discussed at https://r-spatial.org/sdsl/.

GeoInterface seems the correct place to implement something to improve the state of the ecosystem. Other libraries distinguish between Geography and Geometry types (PostGIS), others solve it purely using the attached coordinate reference system. I propose to use traits to distinguish between the two modes, forcing libraries to choose what to implement.

abstract type AbstractCRSTrait end
struct ProjectedTrait <: AbstractCRSTrait end
struct GeographicTrait <: AbstractCRSTrait end

GeoInterface.crstrait(geom) = crstrait(crs(geom))  # One of ProjectedTrait() or GeographicTrait()
crstrait(any::Any) = ProjectedTrait()  # by default

# all relevant methods, example for `area`
area(geom) = area(crstrait(geom), geomtrait(geom), geom)

# for current implementations
@deprecate area(::ProjectedTrait, geomtrait(geom), geom) area(geomtrait(geom), geom)

area(::GeographicTrait, geomtrait(geom), geom) = MethodError("Please use package X to enable geographic calculations.")

This design is backwards compatible, defaulting to ProjectedTrait(), but warns if you are not explicit in your implementations. Once a geometry/crs gains the GeographicTrait, you get an error instead of wrong answer.

rafaqz commented 8 months ago

Very good.

My main questions I have are:

  1. Which package checks the crs trait (Proj?) or is that backend deoendent?
  2. Do we store the trait on object construction so we dont have to get it at runtime every time
  3. What do we attach the trait to if we do store it?
  4. Should we have an UnkownTrait option to show the case where we just dont know. Some packages will not have the ability to check with Proj.jl, and sometimes we dont have CRS at all. Then packages can default to treating them as ProjectedTrait but maybe warn about it
evetion commented 8 months ago
1. Which package checks the crs trait (Proj?) or is that backend deoendent?

Proj is the logical candidate yes, and we might need a package extension to GeoFormatTypes for it. But in the end it depends on the backend, so we need to provide good defaults.

2. Do we store the trait on object construction so we dont have to get it at runtime every time

This is backend dependent, in our own Geometries we could/should do so. GDAL can check it on runtime (it stores a crs), and things like LibGEOS don't know about CRS so always use the default (Projected).

3. What do we attach the trait to if we do store it?

Also backend dependent. Could be to the CRS object attached to something.

4. Should we have an `UnkownTrait` option to show the case where we just dont know. Some packages will not have the ability to check with Proj.jl, and sometimes we dont have CRS at all. Then packages can default to treating them as `ProjectedTrait` but maybe warn about it

This goes back to a good default (see https://github.com/evetion/GeoDataFrames.jl/issues/58), which is probably Projected, as we should return an EngineeringCRS thing instead of nothing for crs(geom).

rafaqz commented 8 months ago

For 2. I wonder if GeoFormatTypes can allow storing the trait in CRS?

On 4. to clarify Im looking for a difference between a known value and the unknown default, so they could both work the same but the unknown/default can warn. The ability to warn about lack of information is the key part.

evetion commented 8 months ago

For 2. I wonder if GeoFormatTypes can allow storing the trait in CRS?

We'd have to break the type signature probably, or just accept it's runtime for now. Note that I would not see that as blocking, and here we only need the interface.

On 4. to clarify Im looking for a difference between a known value and the unknown default, so they could both work the same but the unknown/default can warn. The ability to warn about lack of information is the key part.

The above proposal already does a depwarn if you still implement the old method. But maybe we should do a warn here as well crstrait(any::Any) = ProjectedTrait() # by default, which disappear once packages define it for their geometry/crs.

rafaqz commented 8 months ago

Hmm not quite what I mean...

I want a package consuming geometries like GeometryOps.jl to know and throw a warning once at the start of some operation, if the trait is not actually known.

Thats why returning UnknownCRSTrait() is better than returning a default.

evetion commented 8 months ago

It's down to who's responsible. GeometryOps doesn't have geometries/crs and thus doesn't define any Projected/GeographicTraits. It just works with them, why would it warn? I'm not questioning how an object gained some Tables.jl trait, I just use it.

We will probably have a very similar discussion if crs(geom) might change to returning some Engineering WKT string instead of nothing when it's not defined.

But anyway, apart from ProjectedTrait gaining a field with fallback=true or similar, we could also do:

abstract type AbstractCRSTrait end
abstract type AbstractProjectedTrait <: AbstractCRSTrait end
abstract type AbstractGeographicTrait <: AbstractCRSTrait end
struct ProjectedTrait <: AbstractProjectedTrait end
struct UnknownCRSTrait <: AbstractProjectedTrait end
struct GeographicTrait <: AbstractGeographicTrait end
felixcremer commented 8 months ago

I like the solution of having the UnknownCRSTrait, that way LibGeos can fall back to ProjectedTrait and GDAL can use UnknownCRSTrait when the CRS information is missing, but for GDAL data we would always expect that we need to know the CRS to distinguish between projected and geographic.

From what I understand these are possible consumers of this trait: Analysis: Geomorphology.jl
GeometryOps.jl Plotting We could change the plotting behaviour to use GeoMakie or at least warn that the geometry is supposed to be on the sphere.