feature request: setcrs

[alex-s-gardner]

A simple setcrs utility would be a nice to have. Doing this:

GI.MultiPolygon.(multipoly, crs=GFT.EPSG[3413])
GI.Polygon.(poly, crs=GFT.EPSG[3413])

is a bit awkward. a setcrs utility

GI.setcrs(multipoly, GFT.EPSG[3413])
GI.setcrs(poly, GFT.EPSG[3413])

would help clean things up

[rafaqz]

Nice. Then rasters can extend that method too.

It should work on any geometry but return a GI geometry.

[evetion]

I have it as crs! in GeoArrays, which is more Julian? But yeah, would be good to have as a common definition (like read/write). Not sure about it always wrapping if the object doesn't implement it.

[alex-s-gardner]

I like crs or crs! if it doesn't reallocate

[asinghvi17]

IMO crs!(x) only makes sense if x is mutable and can change its CRS, which is not the case for LibGEOS, Shapefile, GeoJSON, WKB, all other Julia reader geometries, and GI wrappers (I'm probably missing a few here). crs is also used as a getter so I'd tend to prefer setcrs (and maybe optional setcrs! for e.g. ArchGDAL and GeoDataFrames / Tables in general).

[alex-s-gardner]

I guess if we went down this road we might want to introduce:

GI.getcrs and then phase out GI.crs

[rafaqz]

I'm ok with crs and setcrs. That's what Rasters always had

[rafaqz]

setcrs can have the fallback of just wrapping with GI geoms for when the original object can't hold crs

[alex-s-gardner]

Does setcrs need to be defined for all types i.e.:

function setcrs(geom::GeoInterface.Polygon, crs)
    return GeoInterface.Polygon(geom; crs)
end

function setcrs(geom::GeoInterface.MultiPolygon, crs)
    return GeoInterface.MultiPolygon(geom; crs)
end

function setcrs(geom::Missing, crs)
    return Missing
end

or is there a clever way to automate this?

[asinghvi17]

Here's a brief prototype that we can override for known types:

import GeoInterface as GI

setcrs(x, crs) = setcrs(GI.trait(x), x, crs)

function setcrs(trait::Nothing, x, crs)
    error("Trait not found for $(typeof(x))")
    # or, use GeometryOpsCore at the top level
end

function setcrs(trait::GI.AbstractGeometryTrait, x, crs)
    return GI.Wrappers.geointerface_geomtype(trait)(x; crs)
end

# then, usage is

import LibGEOS as LG

p = LG.Point([1., 2.])
setcrs(p, "EPSG:4326")

[asinghvi17]

We can then dispatch on the type of x, or use GeometryOpsCore to make setcrs work on tables etc.

[alex-s-gardner]

@asinghvi17 are you planning to open a PR for the above?

[rafaqz]

This is a one liner with apply in GeometryOpsCore for pretty much anything.

Possibly apply should just live here

[evetion]

I'm hesitant to make this a GeoInterface (1.x) thing that we can't break. With setcrs you can expect one of two behaviours, purely set the crs, or actually transform from the current one to the new one.

That said this seems more a feature request on wrapping a crs for GeoInterface.Wrappers, not a general functionality.

Geopandas seems to have to_crs for the transform and set_crs for this proposal here.

[alex-s-gardner]

purely set the crs, or actually transform from the current one to the new one.

So do we just need to check if a crs already exists, in which case setcrs would throw an error?

[asinghvi17]

That sounds like a good idea

[asinghvi17]

Maybe with a kwarg force that allows you to force the CRS, for example if you've read a GeoJSON (default crs is latlong) that you want to treat as some other coordinate system.