Basic support for curved geometries

src/LibGEOS.jl

@@ -8,12 +8,17 @@ using CEnum
 
 const GI = GeoInterface
 
-export GeometryCollection,
+export CircularString,
+    CompoundCurve,
+    CurvePolygon,
+    GeometryCollection,
     LineString,
     LinearRing,
+    MultiCurve,
     MultiLineString,
     MultiPoint,
     MultiPolygon,
+    MultiSurface,
     Point,
     Polygon,
     STRtree

src/geo_interface.jl

@@ -8,6 +8,11 @@
 geointerface_geomtype(::PolygonTrait) = Polygon
 geointerface_geomtype(::MultiPolygonTrait) = MultiPolygon
 geointerface_geomtype(::GeometryCollectionTrait) = GeometryCollection
+geointerface_geomtype(::CircularStringTrait) = CircularString
+geointerface_geomtype(::CompoundCurveTrait) = CompoundCurve
+geointerface_geomtype(::CurvePolygonTrait) = CurvePolygon
+geointerface_geomtype(::MultiSurfaceTrait) = MultiSurface
+geointerface_geomtype(::MultiCurveTrait) = MultiCurve
 
 GeoInterface.geomtrait(::Point) = PointTrait()
 GeoInterface.geomtrait(::MultiPoint) = MultiPointTrait()
@@ -17,23 +22,33 @@
 GeoInterface.geomtrait(::Polygon) = PolygonTrait()
 GeoInterface.geomtrait(::MultiPolygon) = MultiPolygonTrait()
 GeoInterface.geomtrait(::GeometryCollection) = GeometryCollectionTrait()
+GeoInterface.geomtrait(::CircularString) = CircularStringTrait()
+GeoInterface.geomtrait(::CompoundCurve) = CompoundCurveTrait()
+GeoInterface.geomtrait(::CurvePolygon) = CurvePolygonTrait()
+GeoInterface.geomtrait(::MultiSurface) = MultiSurfaceTrait()
+GeoInterface.geomtrait(::MultiCurve) = MultiCurveTrait()
 GeoInterface.geomtrait(geom::PreparedGeometry) = GeoInterface.geomtrait(geom.ownedby)
 
 GeoInterface.isempty(::AbstractGeometryTrait, geom::AbstractGeometry) = isEmpty(geom)
 
 GeoInterface.ngeom(::AbstractGeometryCollectionTrait, geom::AbstractMultiGeometry) =
     isEmpty(geom) ? 0 : Int(numGeometries(geom))
 GeoInterface.ngeom(::LineStringTrait, geom::LineString) = Int(numPoints(geom))
+GeoInterface.ngeom(::CircularStringTrait, geom::CircularString) = Int(numPoints(geom))
 GeoInterface.ngeom(::LinearRingTrait, geom::LinearRing) = Int(numPoints(geom))
 GeoInterface.ngeom(::PolygonTrait, geom::Polygon) = Int(numInteriorRings(geom)) + 1
+GeoInterface.ngeom(::CurvePolygonTrait, geom::CurvePolygon) =
+    Int(numInteriorRings(geom)) + 1
 GeoInterface.ngeom(t::AbstractGeometryTrait, geom::PreparedGeometry) =
     GeoInterface.ngeom(t, geom.ownedby)
 GeoInterface.ngeom(::AbstractPointTrait, geom::Point) = 0
 GeoInterface.ngeom(::AbstractPointTrait, geom::PreparedGeometry) = 0
 
 GI.is3d(::AbstractGeometryTrait, geom::AbstractGeometry) = hasZ(geom)
-GI.getexterior(::AbstractPolygonTrait, geom::Polygon) = exteriorRing(geom)
-GI.gethole(::AbstractPolygonTrait, geom::Polygon, n) = interiorRing(geom, n)
+GI.getexterior(::AbstractPolygonTrait, geom::Union{Polygon,CurvePolygon}) =
+    exteriorRing(geom)
+GI.gethole(::AbstractPolygonTrait, geom::Union{Polygon,CurvePolygon}, n) =
+    interiorRing(geom, n)
 
 function GeoInterface.getgeom(
     ::AbstractGeometryCollectionTrait,
@@ -60,6 +75,13 @@
         interiorRing(geom, i - 1)
     end
 end
+function GeoInterface.getgeom(::CurvePolygonTrait, geom::CurvePolygon, i::Int)
+    if i == 1
+        exteriorRing(geom)
+    else
+        interiorRing(geom, i - 1)
+    end
+end
 GeoInterface.getgeom(t::AbstractGeometryTrait, geom::PreparedGeometry, i) =
     GeoInterface.getgeom(t, geom.ownedby, i)
 GeoInterface.getgeom(t::AbstractPointTrait, geom::PreparedGeometry, i) = 0

src/geos_functions.jl

@@ -555,6 +555,38 @@
     result
 end
 
+function createCircularString(ptr::GEOSCoordSeq, context::GEOSContext = get_global_context())
+    result = GEOSGeom_createCircularString_r(context, ptr)
+    if result == C_NULL
+        error("LibGEOS: Error in GEOSGeom_createCircularString")
+    end
+    result
+end
+createCircularString(
+    coords::Vector{Vector{Float64}},
+    context::GEOSContext = get_global_context(),
+) = GEOSGeom_createCircularString_r(context, createCoordSeq(coords, context))
+
+# Second argument is an array of GEOSGeometry* objects.
+# The caller remains owner of the array, but pointed-to
+# objects become ownership of the returned GEOSGeometry.
+function createCurvePolygon(
+    shell::Union{LinearRing,CircularString,CompoundCurve,GEOSGeom},
+    holes::AbstractVector,
+    context::GEOSContext = get_global_context(),
+)
+    result = GEOSGeom_createCurvePolygon_r(context, shell, holes, length(holes))
+    if result == C_NULL
+        error("LibGEOS: Error in GEOSGeom_createCurvePolygon")
+    end
+    result
+end
+# convenience function to create polygon without holes
+createCurvePolygon(
+    shell::Union{LinearRing,CircularString,CompoundCurve,GEOSGeom},
+    context::GEOSContext = get_global_context(),
+) = createCurvePolygon(shell, GEOSGeom[], context)
+
 function createEmptyCollection(
     geomtype::GEOSGeomTypes,
     context::GEOSContext = get_global_context(),
@@ -1336,7 +1368,10 @@
 end
 
 # Return -1 on exception
-function numInteriorRings(obj::Polygon, context::GEOSContext = get_context(obj))
+function numInteriorRings(
+    obj::Union{Polygon,CurvePolygon},
+    context::GEOSContext = get_context(obj),
+)
     result = GEOSGetNumInteriorRings_r(context, obj)
     if result == -1
         error("LibGEOS: Error in GEOSGetNumInteriorRings")
@@ -1346,7 +1381,7 @@
 
 # Call only on LINESTRING (returns -1 on exception)
 function numPoints(
-    obj::Union{LineString,LinearRing},
+    obj::Union{LineString,LinearRing,CircularString},
     context::GEOSContext = get_context(obj),
 )
     result = GEOSGeomGetNumPoints_r(context, obj)
@@ -1408,16 +1443,59 @@
     end
 end
 
+# Polygon rings are of type LinearRing, and can be inferred, but CurvePolygon rings can be Union{LinearRing,CircularString,CompoundCurve}
+function interiorRing(
+    obj::CurvePolygon,
+    n::Integer,
+    context::GEOSContext = get_context(obj),
+)
+    if !(0 < n <= numInteriorRings(obj, context))
+        error(
+            "LibGEOS: n=$n is out of bounds for CurvePolygon with $(numInteriorRings(obj, context)) interior ring(s)",
+        )
+    end
+    result = GEOSGetInteriorRingN_r(context, obj, n - 1)
+    if result == C_NULL
+        error("LibGEOS: Error in GEOSGetInteriorRingN")
+    end
+    geomFromGEOS(cloneGeom(result, context), context)
+end
+
+function interiorRings(obj::CurvePolygon, context::GEOSContext = get_context(obj))
+    n = numInteriorRings(obj, context)
+    if n == 0
+        return LinearRing[]
+    else
+        return [interiorRing(obj, i, context) for i = 1:n]
+    end
+end
+
+
 # Return NULL on exception, Geometry must be a Polygon.
 # Returned object is a pointer to internal storage: it must NOT be destroyed directly.
-function exteriorRing(obj::Polygon, context::GEOSContext = get_context(obj))
+function exteriorRing(
+    obj::Polygon,
+    context::GEOSContext = get_context(obj),
+)::LinearRing
     result = GEOSGetExteriorRing_r(context, obj)
     if result == C_NULL
         error("LibGEOS: Error in GEOSGetExteriorRing")
     end
     LinearRing(cloneGeom(result, context), context)
 end
 
+function exteriorRing(
+    obj::CurvePolygon,
+    context::GEOSContext = get_context(obj),
+)::Union{LinearRing,CircularString,CompoundCurve}
+    result = GEOSGetExteriorRing_r(context, obj)
+    if result == C_NULL
+        error("LibGEOS: Error in GEOSGetExteriorRing")
+    end
+    geomFromGEOS(cloneGeom(result, context), context)
+end
+
+
 # Return -1 on exception
 function numCoordinates(obj::Geometry, context::GEOSContext = get_context(obj))
     result = GEOSGetNumCoordinates_r(context, obj)