Add curved geometry definitions to geo_interface.jl

Author: jaakkor2

src/geo_interface.jl

@@ -8,6 +8,11 @@ geointerface_geomtype(::LinearRingTrait) = LinearRing
 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(::LinearRing) = LinearRingTrait()
 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 @@ function GeoInterface.getgeom(::PolygonTrait, geom::Polygon, i::Int)
         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

@@ -1336,7 +1336,10 @@ function normalize!(obj::Geometry, context::GEOSContext = get_context(obj))
 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 +1349,7 @@ end
 
 # 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,14 +1411,46 @@ function interiorRings(obj::Polygon, context::GEOSContext = get_context(obj))
     end
 end
 
+# Polygon rings are of type LinearRingcan, and can be inferred, but CurvePolygon rings can be Union{LinearRing, 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::Union{Polygon,CurvePolygon},
+    context::GEOSContext = get_context(obj),
+)::Union{LinearRing,CompoundCurve}
     result = GEOSGetExteriorRing_r(context, obj)
     if result == C_NULL
         error("LibGEOS: Error in GEOSGetExteriorRing")
     end
-    LinearRing(cloneGeom(result, context), context)
+    geomFromGEOS(cloneGeom(result, context), context)
+    #    LinearRing(cloneGeom(result, context), context)
 end
 
 # Return -1 on exception

test/test_geos_types.jl

@@ -461,4 +461,18 @@ end
 
     @test all(geomLength.(geos) .≈ [pi, 4pi, 4 + 2 * sqrt(2), 4pi, pi + 4, pi + 4, 4pi, 1])
     @test all(area.(geos) .≈ [0, 0, 2, 2pi, 0, pi / 2 + 2, 2pi, 0])
+
+    cp = readgeom(
+        "CURVEPOLYGON (COMPOUNDCURVE (CIRCULARSTRING (0 0, 1 1, 2 0), (2 0, 2 -1, 0 -1, 0 0)))",
+    )
+    @test_broken isequal(cp, CurvePolygon(cp))
+
+    cp = readgeom(
+        "CURVEPOLYGON (COMPOUNDCURVE (CIRCULARSTRING (0 0, 2 2, 4 0), (4 0, 4 -1, 0 -1, 0 0)), LINEARRING (1 0, 3 0, 3 1, 1 1, 1 0))",
+    )
+    @test area(cp) ≈ 2pi + 2
+    @test geomLength(cp) ≈ 2pi + 12
+    @test length(interiorRings(cp)) == 1
+    @test typeof(exteriorRing(cp)) == CompoundCurve
+    @test typeof(interiorRings(cp)[1]) == LinearRing
 end