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.

4 Likes
3

awesome, thanks!