# Set spatial reference system to a geodataframe and plot it

**URL:** <https://discourse.julialang.org/t/set-spatial-reference-system-to-a-geodataframe-and-plot-it/108086>\
**Category:** Geo\
**Created:** [December 27, 2023, 7:01am UTC](https://discourse.julialang.org/t/set-spatial-reference-system-to-a-geodataframe-and-plot-it/108086 "2023-12-27T07:01:48Z")\
**Posts on this page:** 9\
**Page:** 1

<div class="post-metadata">

**Author:** ![anjelinejeline](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/anjelinejeline/32/205552_2.png) [@anjelinejeline](https://discourse.julialang.org/u/anjelinejeline)\
**Post date:** [December 27, 2023, 7:01am UTC](https://discourse.julialang.org/t/set-spatial-reference-system-to-a-geodataframe-and-plot-it/108086/1 "2023-12-27T07:01:48Z")

</div>

Hello,  
I am trying to create a geotable from a csv and plot it along with a shapefile,  
However I am not able to to set the right spatial reference system as I am getting a method error when using the setspatialref! function of ArchGDAL …

Here is my code

```julia
# Load libraries 
using Shapefile, CSV, DataFrames
using GeoDataFrames, GeoFormatTypes
using ArchGDAL, GeoStats, Makie

kg=GeoDataFrames.read("kg.shp")

# Check if the Spatial reference system is the Malloweide 
spatialref=ArchGDAL.getspatialref(kg.geometry[1])
spatialref

Spatial Reference System: +proj=moll +lon_0=0 +x_0=0 +y_0=0 + ... _defs

# Read the df 
df=CSV.read("df.csv", DataFrame)
# Discard the first column
df = select(df, Not(1))

# georeference the table
gf=GeoStats.georef(df,(:x, :y))

# Define the coords reference system 
gf= ArchGDAL.setspatialref!(gf.geometry,spatialref)

{
	"name": "LoadError",
	"message": "MethodError: no method matching setspatialref!(::PointSet{2, Float64}, ::ArchGDAL.ISpatialRef)

Closest candidates are:
  setspatialref!(::ArchGDAL.GeomFieldDefn, ::ArchGDAL.AbstractSpatialRef)
   @ ArchGDAL ~/.julia/packages/ArchGDAL/R3wJR/src/ogr/fielddefn.jl:402
",
	"stack": "MethodError: no method matching setspatialref!(::PointSet{2, Float64}, ::ArchGDAL.ISpatialRef)

Closest candidates are:
  setspatialref!(::ArchGDAL.GeomFieldDefn, ::ArchGDAL.AbstractSpatialRef)
   @ ArchGDAL ~/.julia/packages/ArchGDAL/R3wJR/src/ogr/fielddefn.jl:402

Stacktrace:
 [1] top-level scope
   @ In[14]:2"
}

```

If I overcome this issue I would like to plot them, and join them afterwards (spatial intersection)

```julia
fig = Makie.Figure()
ax = Makie.Axis(fig[1,1], title = "Points and kg",
              xlabel = "longitude", ylabel = "latitude")
viz!(ClimateZones.geometry)
viz!(gf.geometry, color = "teal", pointsize = 10)
fig

```

Thanks for your help

---

<div class="post-metadata">

**Author:** ![juliohm](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/juliohm/32/215266_2.png) [@juliohm](https://discourse.julialang.org/u/juliohm)\
**Post date:** [December 27, 2023, 11:36am UTC](https://discourse.julialang.org/t/set-spatial-reference-system-to-a-geodataframe-and-plot-it/108086/2 "2023-12-27T11:36:32Z")

</div>

Hi @anjelinejeline , welcome to the community.

If you are interested in the GeoStats.jl framework, please read the book

[https://juliaearth.github.io/geospatial-data-science-with-julia](https://juliaearth.github.io/geospatial-data-science-with-julia)

There you will find the recommended usage, and how to interface with other GIS software.

Regarding your specific question, we do not yet support getting/setting the CRS. The new coordinate system infrastructure is planned for the first semester of 2024.

---

<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:** [December 27, 2023, 1:01pm UTC](https://discourse.julialang.org/t/set-spatial-reference-system-to-a-geodataframe-and-plot-it/108086/3 "2023-12-27T13:01:31Z")

</div>

Seems to me this is more of a question for @evetion and GeoDataFrames.jl

---

<div class="post-metadata">

**Author:** ![anjelinejeline](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/anjelinejeline/32/205552_2.png) [@anjelinejeline](https://discourse.julialang.org/u/anjelinejeline)\
**Post date:** [December 27, 2023, 1:07pm UTC](https://discourse.julialang.org/t/set-spatial-reference-system-to-a-geodataframe-and-plot-it/108086/4 "2023-12-27T13:07:05Z")

</div>

@juliohm thanks for your rapid response,  
Can I perform a spatial joint without assigning the crs ?  
I saw that there’s an example [Geospatial Data Science with Julia - 9&nbsp; Geospatial joins](https://juliaearth.github.io/geospatial-data-science-with-julia/09-geojoins.html)

I was not able to plot gf and kg using the viz function that’s why I thought that I needed to define the crs of gf … I got an error related to the receipt that was absent …

Will it be possible to perform a spatial join between gf and kg though the spatial reference system of gf is not recognised? I can see that in the example a join is done but I am not sure how this works…

---

<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:** [December 27, 2023, 1:12pm UTC](https://discourse.julialang.org/t/set-spatial-reference-system-to-a-geodataframe-and-plot-it/108086/5 "2023-12-27T13:12:41Z")

</div>

Spatial intersection you can do with LibGeos.jl `intersection`, you should be able to pass your geometries from GeoDataFrames.jl `geometry` column to it, or at least map over the individual geometries.

Also note you will not need `viz` to plot these as soon as this ArchGDAL PR is merged, they will just plot with `Makie.plot` [Add GeoInterfaceMakie as a weak dependency by felixcremer · Pull Request #404 · yeesian/ArchGDAL.jl · GitHub](https://github.com/yeesian/ArchGDAL.jl/pull/404)

(I will in fact merge and register that now!)

---

<div class="post-metadata">

**Author:** ![juliohm](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/juliohm/32/215266_2.png) [@juliohm](https://discourse.julialang.org/u/juliohm)\
**Post date:** [December 27, 2023, 1:25pm UTC](https://discourse.julialang.org/t/set-spatial-reference-system-to-a-geodataframe-and-plot-it/108086/6 "2023-12-27T13:25:29Z")

</div>

> [@anjelinejeline](#):
>
> Can I perform a spatial joint without assigning the crs ?

That will be possible in the future. So far we don’t have the CRS explicitly modeled in the framework.

> [@anjelinejeline](#):
>
> I was not able to plot gf and kg using the viz function that’s why I thought that I needed to define the crs of gf … I got an error related to the receipt that was absent …

The `gf` and `kg` you are using in your example are from ArchGDAL.jl and GeoDataFrames.jl, which are wrappers to external GIS libraries written in other languages. We do not support these data structures in our framework.

Also, notice that the `viz` function is equivalent to the `Makie.plot` function if you are using `GeoTable` as explained in the book.

> [@anjelinejeline](#):
>
> I can see that in the example a join is done but I am not sure how this works…

Our `geojoin` should work regardless of the CRS, with any kind of geometry, including geometries that are 3D. The problem is that we are still working on the new coordinate system infrastructure to hide these CRS issues from end-users.

My suggestion is to change the CRS of all objects to a single CRS, using whatever method you want, write the results to disk, and then use GeoIO.jl to load the files in the recommended data structure for geospatial data science with GeoStats.jl if you are interested in using the framework at some point.

---

<div class="post-metadata">

**Author:** ![anjelinejeline](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/anjelinejeline/32/205552_2.png) [@anjelinejeline](https://discourse.julialang.org/u/anjelinejeline)\
**Post date:** [December 29, 2023, 3:34pm UTC](https://discourse.julialang.org/t/set-spatial-reference-system-to-a-geodataframe-and-plot-it/108086/7 "2023-12-29T15:34:54Z")

</div>

Hi @Raf thanks for jumping in …

I just give it a try with

`LibGEOS.intersection(gf.geometry,kg.geometry)`

but I got a` MethodError: no method matching geointerface_geomtype(::Nothing)`

I am also wondering whether it wil be possible to retain all columns in the geodataframe given that I only pass the geometry column to the intersection

---

<div class="post-metadata">

**Author:** ![anjelinejeline](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/anjelinejeline/32/205552_2.png) [@anjelinejeline](https://discourse.julialang.org/u/anjelinejeline)\
**Post date:** [December 29, 2023, 4:21pm UTC](https://discourse.julialang.org/t/set-spatial-reference-system-to-a-geodataframe-and-plot-it/108086/8 "2023-12-29T16:21:44Z")

</div>

> [@anjelinejeline](#):
>
> LibGEOS.intersection(gf.geometry,kg.geometry)

R /my message  
If I create a geodataframe using ArchGDAL.createpoint instead of GeoStats.georef and I take a single element in the geometry column the intersection works

```julia
gf=DataFrame(geometry=ArchGDAL.createpoint.(df.x, df.y), year=df.year)
LibGEOS.intersection(gf.geometry[1],kg.geometry[1])

```

So I guess I should iterate over the geometry column of each object and intersect them … which may not be straightforward but how can I keep the year column of gf then?  
Is there any feature that allows to easily intersect two objects while keeping all attributes and not only the geometry?

Thanks for the useful discussion

---

<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:** [December 29, 2023, 6:11pm UTC](https://discourse.julialang.org/t/set-spatial-reference-system-to-a-geodataframe-and-plot-it/108086/9 "2023-12-29T18:11:48Z")

</div>

Yeah looks like you have to iterate over them manually. But being julia its probably trivial, like just broadcast the `intersection` function (note the single added `.`) and assign the result to a column:

```julia
gf.intersection = LibGEOS.intersection.(gf.geometry,kg.geometry)

```

Sometime in the next year we will have native julia versions of all these things and they be a bit smoother

> **[GitHub - asinghvi17/GeometryOps.jl: GeoInterface-based geometry operations](https://github.com/asinghvi17/GeometryOps.jl)**
>
> GeoInterface-based geometry operations. Contribute to asinghvi17/GeometryOps.jl development by creating an account on GitHub.