Shapefiles read with v0.4.0 can no longer use common GeoInterface methods

[niclasmattsson]

I'm not sure if the diagnosis in my title is correct, feel free to edit it when you figure out what the problem is. Anyway, my code utterly broke when I upgraded to v0.4.0. The code reads and manipulates geodata in various formats, both geopackages and shapefiles. The parts that operate exclusively on geopackages work as before in v0.3.13, but everything that interacts with shapefiles now produce method errors as shown below. Basic methods like getgeom work fine, but optional methods like area, contains or intersects don't work anymore.

Here's an example using a shapefile of Austria from GADM level 0.

v0.4.0

julia> df = GDF.read("C:/Users/niclas/Downloads/gadm41_AUT_shp/gadm41_AUT_0.shp")
1×3 DataFrame
 Row │ geometry                GID_0    COUNTRY
     │ Polygon?                String?  String?
─────┼──────────────────────────────────────────
   1 │ Polygon(249273 Points)  AUT      Austria

julia> GI.area(df.geometry[1])
ERROR: MethodError: no method matching area(::GeoInterface.MultiPolygonTrait, ::Shapefile.Polygon)
The function `area` exists, but no method is defined for this combination of argument types.

Closest candidates are:
  area(::Any)
   @ GeoInterface C:\Users\niclas\.julia\packages\GeoInterface\4tyIo\src\interface.jl:244
  area(::Union{GeoInterface.CircularStringTrait, GeoInterface.CompoundCurveTrait, GeoInterface.CurvePolygonTrait, GeoInterface.GeometryCollectionTrait, GeoInterface.LineStringTrait, GeoInterface.LinearRingTrait, GeoInterface.MultiLineStringTrait, GeoInterface.MultiPointTrait, GeoInterface.MultiPolygonTrait, GeoInterface.MultiSurfaceTrait, GeoInterface.PointTrait, GeoInterface.PolygonTrait, GeoInterface.PolyhedralSurfaceTrait, GeoInterface.TINTrait, GeoInterface.TriangleTrait}, ::ArchGDAL.AbstractGeometry)
   @ ArchGDAL C:\Users\niclas\.julia\packages\ArchGDAL\ujstt\src\geointerface.jl:310

Stacktrace:
 [1] area(geom::Shapefile.Polygon)
   @ GeoInterface C:\Users\niclas\.julia\packages\GeoInterface\4tyIo\src\interface.jl:244
 [2] top-level scope
   @ REPL[178]:1

v0.3.13

julia> df = GDF.read("C:/Users/niclas/Downloads/gadm41_AUT_shp/gadm41_AUT_0.shp")
1×3 DataFrame
 Row │ geometry                   GID_0   COUNTRY
     │ IGeometr…                  String  String
─────┼────────────────────────────────────────────
   1 │ Geometry: wkbMultiPolygon  AUT     Austria

julia> GI.area(df.geometry[1])
10.038532695098086

[asinghvi17]

That was the point of 0.4.0 iirc. Now everything is forwarded to GeometryOps or package specific methods. I don't remember if we added a changelog but this should be in docs/tutorials. Will forward some links when I find them.

[niclasmattsson]

Thank you, that comment led me on the right track. The missing methods are available by prefixing with GeometryOps, e.g. GDF.GeometryOps.area. I wasn't aware that this reshuffling of methods into GeometryOps was taking place. I noticed the "Now exports GeometryOps" line in the release notes but didn't understand its significance.

I think an explicit line about this would be helpful in the GeoDataFrames README, but otherwise this issue can be closed.

[niclasmattsson]

I still think there are interoperability problems though. When I tried calculating the intersection of multipolygons from different formats I needed to use a wrapper on one argument because the geopackage import ended up in ArchGDAL format while the shapefile became a Julia polygon. I can't tell if this difference is by design.

julia> GO.intersection(df_gpkg.geometry[1], df_shp.geometry[250]; target=GI.MultiPolygonTrait())
ERROR: MethodError: Cannot `convert` an object of type
  GeoInterface.Wrappers.Polygon{false,false,Array{GeoInterface.Wrappers.LineString{false, false, Vector{Tuple{Float64, Float64}}, Nothing, GeoFormatTypes.WellKnownText{GeoFormatTypes.CRS}},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, Extents.Extent}, C})(::Any, ::Any, ::Any)
   @ GeoInterface C:\Users\niclas\.julia\packages\GeoInterface\4tyIo\src\wrappers.jl:166
  convert(::Type{T}, ::T) where T
   @ Base Base.jl:126

julia> GO.intersection(GI.MultiPolygon(df_gpkg.geometry[1]), df_shp.geometry[250]; target=GI.MultiPolygonTrait())
GeoInterface.Wrappers.MultiPolygon{false, false}([GeoInterface.Wrappers.Polygon([GeoInterface.Wrappers.LinearRing([(537000.0, 6.234e6), … (3) … , (537000.0, 6.234e6)])])])

Is this intended?

[asinghvi17]

Ah, that's a long standing GeometryOps bug with CRS in geometries. I do need to clean that up.

But for now you could also use ArchGDAL methods (which is what the GeoInterface methods forwarded to previously) in the same way as you used the GeoInterface methods before. Just import ArchGDAL as AG and then AG.* instead of GI.*.

[niclasmattsson]

Thanks, I already solved my problem though - see the line below the error in my last post. But I tried your suggestion yesterday and ran into similar interoperability issues. I experimented with several variations of it but never got it to work.

julia> AG.intersection(df_gpkg.geometry[1], df_shp.geometry[250])
ERROR: MethodError: no method matching intersection(::ArchGDAL.IGeometry{ArchGDAL.wkbPolygon}, ::Shapefile.Polygon)
The function `intersection` exists, but no method is defined for this combination of argument types.

Closest candidates are:
  intersection(::Function, ::Any...; kwargs...)
   @ ArchGDAL C:\Users\niclas\.julia\packages\ArchGDAL\ujstt\src\context.jl:265
  intersection(::ArchGDAL.AbstractGeometry, ::ArchGDAL.AbstractGeometry)
   @ ArchGDAL C:\Users\niclas\.julia\packages\ArchGDAL\ujstt\src\ogr\geometry.jl:806

Problems aside, I appreciate that GeoDataFrames is moving towards pure Julia representations of the data structures. It's so much easier to figure out how things are supposed to work when you can dump an object to see its internals instead of immediately running into C pointers.

[technocrat]

I had the same troubles with shapefiles after an update. While GeoDataFrames.read("my_shapefile.shp") brings the geometry in as typeof IGeometry, plotting with CairoMakie/GeoMakie fails with

caused by: MethodError: no method matching convert_arguments(::PointBased, ::Vector{Union{Missing, Shapefile.Polygon}}) The functionconvert_arguments exists, but no method is defined for this combination of argument types.

I've ended up with a workaround for anyone who needs one while things settle out.

using CSV, CairoMakie, Colors, DataFrames, Dates, GeoDataFrames, GeoMakie, 
GeoInterface, GeometryBasics

include("src/constants.jl")

tigerline_file = "data/2024_shp/cb_2024_us_state_500k.shp"
geo = GeoDataFrames.read(tigerline_file)

f = Figure(size = (1200, 600))
ga = GeoAxis(f[1, 1], dest = conus_crs)

# Comprehensive conversion function for all geometry types
function convert_to_makie_polygons(geom)
    try
        coords = GeoInterface.coordinates(geom)
        geom_type = GeoInterface.geomtrait(geom)
        
        polygons = GeometryBasics.Polygon{2, Float32}[]
        
        if geom_type == GeoInterface.PolygonTrait()
            # Single polygon (possibly with holes)
            poly = convert_single_polygon(coords)
            if poly !== nothing
                push!(polygons, poly)
            end
            
        elseif geom_type == GeoInterface.MultiPolygonTrait()
            # Multiple polygons
            for poly_coords in coords
                poly = convert_single_polygon(poly_coords)
                if poly !== nothing
                    push!(polygons, poly)
                end
            end
        end
        
        return polygons
        
    catch e
        println("Error converting geometry: $e")
        return GeometryBasics.Polygon{2, Float32}[]
    end
end

function convert_single_polygon(coords)
    try
        if isempty(coords)
            return nothing
        end
        
        # Get exterior ring (first ring)
        exterior_coords = coords[1]
        
        # Convert to Point2f
        points = Point2f[]
        for coord in exterior_coords
            if length(coord) >= 2 && coord[1] isa Real && coord[2] isa Real
                push!(points, Point2f(Float32(coord[1]), Float32(coord[2])))
            end
        end
        
        # Remove duplicate last point if it exists
        if length(points) > 1 && points[1] ≈ points[end]
            pop!(points)
        end
        
        # Need at least 3 points for a valid polygon
        if length(points) >= 3
            # Handle holes if they exist
            if length(coords) > 1
                holes = Point2f[][]
                for hole_coords in coords[2:end]
                    hole_points = Point2f[]
                    for coord in hole_coords
                        if length(coord) >= 2 && coord[1] isa Real && coord[2] isa Real
                            push!(hole_points, Point2f(Float32(coord[1]), Float32(coord[2])))
                        end
                    end
                    # Remove duplicate last point
                    if length(hole_points) > 1 && hole_points[1] ≈ hole_points[end]
                        pop!(hole_points)
                    end
                    if length(hole_points) >= 3
                        push!(holes, hole_points)
                    end
                end
                return GeometryBasics.Polygon(points, holes)
            else
                return GeometryBasics.Polygon(points)
            end
        else
            return nothing
        end
        
    catch e
        return nothing
    end
end

# Convert all geometries to individual polygons
all_polygons = GeometryBasics.Polygon{2, Float32}[]
skipped_count = 0

for (i, geom) in enumerate(geo.geometry)
    converted_polys = convert_to_makie_polygons(geom)
    if isempty(converted_polys)
        skipped_count += 1
        println("Skipped geometry $i")
    else
        append!(all_polygons, converted_polys)
    end
end

println("Successfully converted $(length(all_polygons)) polygons from $(length(geo.geometry)) input geometries")
println("Skipped $skipped_count geometries")

# Plot all polygons
if !isempty(all_polygons)
    poly!(ga, all_polygons, color = :white, strokecolor = :black, strokewidth = 0.25)
else
    println("No valid polygons to plot!")
end

display(f)

[asinghvi17]

Ah that's probably a missing dispatch in the makie extension for union{missing, geom}. If you know there are no missing geoms just reset the array type by identity.(geoms)