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

**URL:** <https://discourse.julialang.org/t/how-to-transform-crs-of-entire-shapefile-not-a-single-point/84489>\
**Category:** Geo\
**Tags:** geo, gdal, geodesy\
**Created:** [July 19, 2022, 8:30pm UTC](https://discourse.julialang.org/t/how-to-transform-crs-of-entire-shapefile-not-a-single-point/84489 "2022-07-19T20:30:31Z")\
**Posts on this page:** 7\
**Page:** 1

<div class="post-metadata">

**Author:** ![floswald](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/floswald/32/195_2.png) [@floswald](https://discourse.julialang.org/u/floswald)\
**Post date:** [July 19, 2022, 8:30pm UTC](https://discourse.julialang.org/t/how-to-transform-crs-of-entire-shapefile-not-a-single-point/84489/1 "2022-07-19T20:30:31Z")

</div>

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.

```julia
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](https://stackoverflow.com/questions/63677035/is-there-an-julia-function-to-transform-coordinate-reference-system-e-g-st-tr).

---

<div class="post-metadata">

**Author:** ![visr](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/visr/32/17204_2.png) [@visr](https://discourse.julialang.org/u/visr)\
**Post date:** [July 19, 2022, 10:22pm UTC](https://discourse.julialang.org/t/how-to-transform-crs-of-entire-shapefile-not-a-single-point/84489/2 "2022-07-19T22:22:33Z")

</div>

This will do the trick:

```julia
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](https://yeesian.com/ArchGDAL.jl/dev/projections/). I give the example using GeoDataFrames since it builds on ArchGDAL to provide probably the simplest API for what you want.

---

<div class="post-metadata">

**Author:** ![floswald](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/floswald/32/195_2.png) [@floswald](https://discourse.julialang.org/u/floswald)\
**Post date:** [July 21, 2022, 9:58am UTC](https://discourse.julialang.org/t/how-to-transform-crs-of-entire-shapefile-not-a-single-point/84489/3 "2022-07-21T09:58:10Z")

</div>

> [@visr](#):
>
> ```julia
> import GeoDataFrames as GDF
> import GeoFormatTypes as GFT
> 
> ```

awesome, thanks!

---

<div class="post-metadata">

**Author:** ![Amine\_o](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/amine_o/32/209936_2.png) [@Amine\_o](https://discourse.julialang.org/u/Amine_o)\
**Post date:** [June 12, 2024, 11:09pm UTC](https://discourse.julialang.org/t/how-to-transform-crs-of-entire-shapefile-not-a-single-point/84489/4 "2024-06-12T23:09:23Z")

</div>

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.

---

<div class="post-metadata">

**Author:** ![asinghvi17](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/asinghvi17/32/8272_2.png) [@asinghvi17](https://discourse.julialang.org/u/asinghvi17)\
**Post date:** [June 16, 2024, 1:55pm UTC](https://discourse.julialang.org/t/how-to-transform-crs-of-entire-shapefile-not-a-single-point/84489/5 "2024-06-16T13:55:35Z")

</div>

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

```julia
DataFrames.metadata(df, "crs")

```

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

---

<div class="post-metadata">

**Author:** ![Raf](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/raf/32/3383_2.png) [@Raf](https://discourse.julialang.org/u/Raf)\
**Post date:** [June 16, 2024, 2:08pm UTC](https://discourse.julialang.org/t/how-to-transform-crs-of-entire-shapefile-not-a-single-point/84489/6 "2024-06-16T14:08:13Z")

</div>

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

```julia
GeoInterface.crs(x) = DataAPI.metadata(x, "crs")

```

---

<div class="post-metadata">

**Author:** ![asinghvi17](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/asinghvi17/32/8272_2.png) [@asinghvi17](https://discourse.julialang.org/u/asinghvi17)\
**Post date:** [June 18, 2024, 6:45pm UTC](https://discourse.julialang.org/t/how-to-transform-crs-of-entire-shapefile-not-a-single-point/84489/7 "2024-06-18T18:45:04Z")

</div>

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.