How to transform CRS of entire shapefile (not a single point)?

Specific Domains Geo
, ,
1

I could not find a way to transform an entire shapefile into another crs. There are too many options for a first time user. ArchGDAL.jl? Proj4? To me it seems reasonable that the Shapefile package should have a transform method tbh.

using Shapefile, Proj4, GeoInterface

# from https://www.data.gouv.fr/fr/datasets/delimitation-des-aires-geographiques-des-siqo/
aoctable = Shapefile.Table("2021-12-22_delim_aire_geographique_shp.shp")
wine_shapes = Shapefile.shapes(aoctable)[aoctable.type_prod .== "Vins"]

wgs84 = Projection(Proj4.epsg[4326])
lam63 = Projection(Proj4.epsg[2154])

transform(lam63,wgs84, GeoInterface.coordinates.(wine_shapes))
ERROR: MethodError: no method matching transform(::Projection, ::Projection, ::Vector{Vector{Vector{Vector{Vector{Float64}}}}})

This is like the first operation I usually have to do on a shapefile, so I wish it were easier to find how to do this. It does not help that exact use case on SO is not answered so far.

2

This will do the trick:

import GeoDataFrames as GDF
import GeoFormatTypes as GFT

df = GDF.read("2021-12-22_delim_aire_geographique_shp.shp")
df.geometry = GDF.reproject(df.geometry, GFT.EPSG(2154), GFT.EPSG(4326))

Shapefile.jl is really just for reading shapefiles in pure julia. For transforming points between any EPSG code we need to rely on PROJ. PROJ is available directly through Proj.jl, but is more geared towards transforming coordinates rather than datasets. GDAL does both reading and writing as well as reprojection, through PROJ. So that is your best option here. ArchGDAL supports it, there are docs here: Spatial Projections · ArchGDAL.jl. I give the example using GeoDataFrames since it builds on ArchGDAL to provide probably the simplest API for what you want.

5 Likes
3

awesome, thanks!

4

Great solution. Is there a way to get the CRS from the Shapefile instead of hardcoding the EPSG number? I tried

GeoInterface.crs(cbsapg.geometry[1])

and it returns the CRS of the geometry but it doesn’t work with GeoDataFrames.reproject.

5

You can get the CRS from the Shapefile.Table object directly via GeoInterface.crs(table). For GeoDataFrames,

DataFrames.metadata(df, "crs")

should work. We should get this hooked up with GI.crs, though.

6

Maybe with a DataAPI.jl extension (or direct dep?) to GeoInterface.jl that adds a fallback like:

GeoInterface.crs(x) = DataAPI.metadata(x, "crs")
7

Extension wouldn’t work on Julia < 1.9, so we could either declare newer versions of GI only compatible with Julia > 1.9, or add it as a direct dep.

I also just made a PR to GeoDataFrames to use the Apache Arrow namespacing thing (GEOINTERFACE:crs) that @bkamins suggested…my thought was that we can start by looking at the namespace metadata, then fall back to crs if that doesn’t exist.

1 Like