Set spatial reference system to a geodataframe and plot it

Specific Domains Geo
1

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

# 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)

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

2

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

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.

3

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

4

@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  Geospatial joins

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…

5

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

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

2 Likes
6

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

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.

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.

7

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

8

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

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

9

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:

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