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

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.

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…

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.

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.

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

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

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?

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: