Intersects between (Multi)Polygons

[evetion]

Based on https://discourse.julialang.org/t/area-interpolation-between-shapefile-country-borders-and-wkt-multi-polygon-elements-in-a-csv-file/133458, and the data mentioned there.

using CSV
using WellKnownGeometry
using GeoFormatTypes
import GeometryOps as GO
import Proj
using GeoDataFrames
import GeoInterface as GI

countries = GeoDataFrames.read(
    "/Users/evetion/Downloads/geoBoundariesCGAZ_ADM0/geoBoundariesCGAZ_ADM0.shp",
)
polygons = CSV.read("test.csv", DataFrame)
polygons.geometry =
    GeoFormatTypes.WellKnownText.((GeoFormatTypes.Geom(),), polygons.Geometry)

polygons.countries_area_km2 = map(polygons.geometry) do geom
    icountries = countries.geometry[GO.intersects.(countries.geometry, (geom,))]
    intersections = GO.intersection.(icountries, (geom,) # FAIL
    sum(GO.area.((Geodesic(),), intersections) / 1e6)  # in km²
end

Without a target kwarg I get

ERROR: MethodError: no method matching TraitTarget(::Nothing)
The type `TraitTarget` exists, but no method is defined for this combination of argument types when trying to construct it.

Closest candidates are:
  TraitTarget(::T) where T<:GeoInterface.AbstractTrait
   @ GeometryOpsCore ~/.julia/packages/GeometryOpsCore/cSZdR/src/types/traittarget.jl:33
  TraitTarget(::GeoInterface.AbstractTrait...)
   @ GeometryOpsCore ~/.julia/packages/GeometryOpsCore/cSZdR/src/types/traittarget.jl:36
  TraitTarget(::Type{<:TraitTarget{T}}) where T
   @ GeometryOpsCore ~/.julia/packages/GeometryOpsCore/cSZdR/src/types/traittarget.jl:35
  ...

Stacktrace:
  [1] intersection(alg::GeometryOps.FosterHormannClipping{…}, geom_a::ArchGDAL.IGeometry{…}, geom_b::GeoInterface.Wrappers.MultiPolygon{…}, ::Type{…}; target::Nothing, kwargs::@Kwargs{})
    @ GeometryOps ~/.julia/packages/GeometryOps/6jExF/src/methods/clipping/intersection.jl:45
  [2] intersection(geom_a::ArchGDAL.IGeometry{…}, geom_b::GeoInterface.Wrappers.MultiPolygon{…}, ::Type{…}; target::Nothing, kwargs::@Kwargs{})
    @ GeometryOps ~/.julia/packages/GeometryOps/6jExF/src/methods/clipping/intersection.jl:55
  [3] intersection(geom_a::ArchGDAL.IGeometry{…}, geom_b::GeoInterface.Wrappers.MultiPolygon{…}, ::Type{…})
    @ GeometryOps ~/.julia/packages/GeometryOps/6jExF/src/methods/clipping/intersection.jl:52
  [4] _broadcast_getindex_evalf

With setting target kwarg to PolygonTrait

ERROR: MethodError: Cannot `convert` an object of type 
  GeoInterface.Wrappers.Polygon{false,false,Array{GeoInterface.Wrappers.LineString{false, false, Vector{Tuple{Float64, Float64}}, Nothing, Nothing},1},Nothing,Nothing} to an object of type 
  GeoInterface.Wrappers.Polygon{false,false,Array{GeoInterface.Wrappers.LinearRing{false, false, Vector{Tuple{Float64, Float64}}, Nothing, Nothing},1},Nothing,Nothing}
The function `convert` exists, but no method is defined for this combination of argument types.

Closest candidates are:
  (::Type{GeoInterface.Wrappers.Polygon{Z, M, T, E, C}} where {Z, M, T, E<:Union{Nothing, Extent}, C})(::Any, ::Any, ::Any)
   @ GeoInterface ~/.julia/packages/GeoInterface/X1jn3/src/wrappers.jl:166
  convert(::Type{T}, ::T) where T
   @ Base Base_compiler.jl:133

Stacktrace:
  [1] push!
    @ ./array.jl:1285 [inlined]
  [2] _intersection(alg::GeometryOps.FosterHormannClipping{…}, ::TraitTarget{…}, ::Type{…}, ::PolygonTrait, poly_a::GeoInterface.Wrappers.Polygon{…}, ::PolygonTrait, poly_b::GeoInterface.Wrappers.Polygon{…}; exact::GeometryOpsCore.True, kwargs::@Kwargs{})
    @ GeometryOps ~/.julia/packages/GeometryOps/6jExF/src/methods/clipping/intersection.jl:89
  [3] _intersection
    @ ~/.julia/packages/GeometryOps/6jExF/src/methods/clipping/intersection.jl:74 [inlined]
  [4] #intersection#130
    @ ~/.julia/packages/GeometryOps/6jExF/src/methods/clipping/intersection.jl:45 [inlined]
  [5] intersection
    @ ~/.julia/packages/GeometryOps/6jExF/src/methods/clipping/intersection.jl:42 [inlined]
  [6] _intersection(alg::GeometryOps.FosterHormannClipping{…}, target::TraitTarget{…}, ::Type{…}, ::PolygonTrait, poly_a::GeoInterface.Wrappers.Polygon{…}, ::MultiPolygonTrait, multipoly_b::ArchGDAL.IGeometry{…}; fix_multipoly::UnionIntersectingPolygons, kwargs::@Kwargs{…})
    @ GeometryOps ~/.julia/packages/GeometryOps/6jExF/src/methods/clipping/intersection.jl:135
  [7] _intersection
    @ ~/.julia/packages/GeometryOps/6jExF/src/methods/clipping/intersection.jl:124 [inlined]
  [8] #intersection#130
    @ ~/.julia/packages/GeometryOps/6jExF/src/methods/clipping/intersection.jl:45 [inlined]
  [9] intersection (repeats 2 times)
    @ ~/.julia/packages/GeometryOps/6jExF/src/methods/clipping/intersection.jl:42 [inlined]
 [10] _intersection
    @ ~/.julia/packages/GeometryOps/6jExF/src/methods/clipping/intersection.jl:143 [inlined]
 [11] #_intersection#138
    @ ~/.julia/packages/GeometryOps/6jExF/src/methods/clipping/intersection.jl:179 [inlined]
 [12] _intersection
    @ ~/.julia/packages/GeometryOps/6jExF/src/methods/clipping/intersection.jl:173 [inlined]
 [13] intersection(alg::GeometryOps.FosterHormannClipping{…}, geom_a::ArchGDAL.IGeometry{…}, geom_b::GeoInterface.Wrappers.Polygon{…}, ::Type{…}; target::MultiPolygonTrait, kwargs::@Kwargs{})
    @ GeometryOps ~/.julia/packages/GeometryOps/6jExF/src/methods/clipping/intersection.jl:45
 [14] intersection (repeats 2 times)
    @ ~/.julia/packages/GeometryOps/6jExF/src/methods/clipping/intersection.jl:42 [inlined]
 [15] #intersection#131
    @ ~/.julia/packages/GeometryOps/6jExF/src/methods/clipping/intersection.jl:55 [inlined]
 [16] intersection
    @ ~/.julia/packages/GeometryOps/6jExF/src/methods/clipping/intersection.jl:52 [inlined]
 [17] (::Base.Broadcast.var"#18#19"{…})(::ArchGDAL.IGeometry{…}, ::Vararg{…})
    @ Base.Broadcast ./broadcast.jl:1333
 [18] _broadcast_getindex_evalf

[asinghvi17]

Yeah this is the same CRS issue again. For a quick fix you can run GO.tuples(geom; crs = nothing) and that should fix it. But we need to figure out how to add correct support for all of this...

[evetion]

So I should set the crs for all wrapper geometries in GeoDataFrames? Is there an issue somewhere with more details? I rather not use GO.tuples in random scripts, it should work out of the box.

[asinghvi17]

No we need to fix it here really, the problem is that you get the wrong type on the array and then it fails to convert. I know how to fix it, just need to find the time...