diff --git a/src/geo_interface.jl b/src/geo_interface.jl index 7e88cb1..6d84518 100644 --- a/src/geo_interface.jl +++ b/src/geo_interface.jl @@ -19,47 +19,59 @@ GeoInterface.geomtrait(::MultiPolygon) = MultiPolygonTrait() GeoInterface.geomtrait(::GeometryCollection) = GeometryCollectionTrait() GeoInterface.geomtrait(geom::PreparedGeometry) = GeoInterface.geomtrait(geom.ownedby) -GeoInterface.ngeom(::AbstractGeometryTrait, geom::AbstractGeometry) = - isEmpty(geom) ? 0 : numGeometries(geom) -GeoInterface.ngeom(::AbstractPointTrait, geom::AbstractGeometry) = 0 +GeoInterface.ngeom(::AbstractGeometryCollectionTrait, geom::AbstractMultiGeometry) = + isEmpty(geom) ? 0 : Int(numGeometries(geom)) +GeoInterface.ngeom(::LineStringTrait, geom::LineString) = Int(numPoints(geom)) +GeoInterface.ngeom(::LinearRingTrait, geom::LinearRing) = Int(numPoints(geom)) +GeoInterface.ngeom(::PolygonTrait, geom::Polygon) = 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 -function GeoInterface.getgeom(::AbstractGeometryTrait, geom::AbstractGeometry, i) +function GeoInterface.getgeom( + ::AbstractGeometryCollectionTrait, + geom::AbstractMultiGeometry, + i, +) getGeometry(geom, i) end - -GeoInterface.getgeom(::AbstractPointTrait, geom::AbstractGeometry, i) = nothing -GeoInterface.ngeom(::AbstractGeometryTrait, geom::Union{LineString,LinearRing}) = - numPoints(geom) -GeoInterface.ngeom(t::AbstractPointTrait, geom::Union{LineString,LinearRing}) = 0 -GeoInterface.getgeom(::AbstractGeometryTrait, geom::Union{LineString,LinearRing}, i) = - Point(getPoint(geom, i)) -GeoInterface.getgeom(::AbstractPointTrait, geom::Union{LineString,LinearRing}, i) = nothing - -GeoInterface.ngeom(::AbstractGeometryTrait, geom::Polygon) = numInteriorRings(geom) + 1 -GeoInterface.ngeom(::AbstractPointTrait, geom::Polygon) = 0 -function GeoInterface.getgeom(::AbstractGeometryTrait, geom::Polygon, i) +GeoInterface.getgeom(::MultiPointTrait, geom::MultiPoint, i) = getGeometry(geom, i)::Point +GeoInterface.getgeom(::MultiLineStringTrait, geom::MultiLineString, i) = + getGeometry(geom, i)::LineString +GeoInterface.getgeom(::MultiPolygonTrait, geom::MultiPolygon, i) = + getGeometry(geom, i)::Polygon +GeoInterface.getgeom( + ::Union{LineStringTrait,LinearRingTrait}, + geom::Union{LineString,LinearRing}, + i, +) = Point(getPoint(geom, i)) +GeoInterface.getgeom(t::AbstractPointTrait, geom::PreparedGeometry) = nothing +function GeoInterface.getgeom(::PolygonTrait, geom::Polygon, i::Int) if i == 1 LinearRing(exteriorRing(geom)) else LinearRing(interiorRing(geom, i - 1)) end end -GeoInterface.getgeom(::AbstractPointTrait, geom::Polygon, i) = nothing - -GeoInterface.ngeom(t::AbstractGeometryTrait, geom::PreparedGeometry) = - GeoInterface.ngeom(t, geom.ownedby) -GeoInterface.ngeom(t::AbstractPointTrait, geom::PreparedGeometry) = 0 GeoInterface.getgeom(t::AbstractGeometryTrait, geom::PreparedGeometry, i) = GeoInterface.getgeom(t, geom.ownedby, i) GeoInterface.getgeom(t::AbstractPointTrait, geom::PreparedGeometry, i) = 0 -GeoInterface.ncoord(::AbstractGeometryTrait, geom::AbstractGeometry) = - isEmpty(geom) ? 0 : getCoordinateDimension(geom) -GeoInterface.getcoord(::AbstractGeometryTrait, geom::AbstractGeometry, i) = - getCoordinates(getCoordSeq(geom), 1)[i] +GeoInterface.coordinates(t::AbstractPointTrait, geom::Point) = collect(getcoord(t, geom)) +GeoInterface.coordinates(t::AbstractPointTrait, geom::AbstractGeometry) = nothing +GeoInterface.coordinates(t::AbstractGeometryTrait, geom::AbstractGeometry) = + [GeoInterface.coordinates(x) for x in getgeom(t, geom)] +GeoInterface.coordinates(t::AbstractGeometryCollectionTrait, geom::AbstractMultiGeometry) = + [GeoInterface.coordinates(x) for x in getgeom(t, geom)] +GeoInterface.ncoord(::AbstractGeometryTrait, geom::AbstractGeometry) = + isEmpty(geom) ? 0 : Int(getCoordinateDimension(geom)) GeoInterface.ncoord(t::AbstractGeometryTrait, geom::PreparedGeometry) = GeoInterface.ncoord(t, geom.ownedby) + +GeoInterface.getcoord(::AbstractGeometryTrait, geom::AbstractGeometry, i) = + getCoordinates(getCoordSeq(geom), 1)[i] GeoInterface.getcoord(t::AbstractGeometryTrait, geom::PreparedGeometry, i) = GeoInterface.getcoord(t, geom.ownedby, i) @@ -67,7 +79,7 @@ GeoInterface.getcoord(t::AbstractGeometryTrait, geom::PreparedGeometry, i) = function GeoInterface.extent(::AbstractGeometryTrait, geom::AbstractGeometry) # minx, miny, maxx, maxy = getExtent(geom) env = envelope(geom) - return Extent(X = (getXMin(env), getXMax(env)), Y = (getYMin(env), getYMax(env))) + return Extent(; X = (getXMin(env), getXMax(env)), Y = (getYMin(env), getYMax(env))) end GI.convert(::Type{Point}, ::PointTrait, geom::Point; context = nothing) = geom @@ -180,7 +192,6 @@ function _geom_to_coord_seq(geom, context) return seq end - GeoInterface.distance( ::AbstractGeometryTrait, ::AbstractGeometryTrait, @@ -268,7 +279,6 @@ GeoInterface.union( GeoInterfaceRecipes.@enable_geo_plots AbstractGeometry - # ----- # LibGeos operations for any GeoInterface.jl compatible geometries # ----- diff --git a/src/geos_types.jl b/src/geos_types.jl index 9bfa779..07825c7 100644 --- a/src/geos_types.jl +++ b/src/geos_types.jl @@ -1,5 +1,13 @@ +"""All Geometries in LibGEOS are an AbstractGeometry.""" abstract type AbstractGeometry end +""" +All MultiGeometries in LibGEOS are an AbstractMultiGeometry. + +Used to specialize on methods that only work on MultiGeometries, like `ngeom`. +""" +abstract type AbstractMultiGeometry <: AbstractGeometry end + function Base.show(io::IO, geo::AbstractGeometry) compact = get(io, :compact, false) if compact @@ -41,7 +49,7 @@ mutable struct Point <: AbstractGeometry Point(createPoint(x, y, z, context), context) end -mutable struct MultiPoint <: AbstractGeometry +mutable struct MultiPoint <: AbstractMultiGeometry ptr::GEOSGeom context::GEOSContext # create a multipoint from a pointer - only makes sense if it is a pointer to a multipoint @@ -111,7 +119,7 @@ mutable struct LineString <: AbstractGeometry end end -mutable struct MultiLineString <: AbstractGeometry +mutable struct MultiLineString <: AbstractMultiGeometry ptr::GEOSGeom context::GEOSContext # create a multiline string from a multilinestring or a linestring pointer, else error @@ -185,7 +193,6 @@ mutable struct LinearRing <: AbstractGeometry finalizer(destroyGeom, ring) ring end - end mutable struct Polygon <: AbstractGeometry @@ -227,7 +234,7 @@ mutable struct Polygon <: AbstractGeometry ) = Polygon(createPolygon(exterior, holes, context), context) end -mutable struct MultiPolygon <: AbstractGeometry +mutable struct MultiPolygon <: AbstractMultiGeometry ptr::GEOSGeom context::GEOSContext # create multipolygon using a multipolygon or polygon pointer, else error @@ -274,7 +281,7 @@ mutable struct MultiPolygon <: AbstractGeometry ) end -mutable struct GeometryCollection <: AbstractGeometry +mutable struct GeometryCollection <: AbstractMultiGeometry ptr::GEOSGeom context::GEOSContext # create a geometric collection from a pointer to a geometric collection, else error @@ -423,7 +430,7 @@ function compare( ng1 = ngeom(geo1) ng2 = ngeom(geo2) ng1 == ng2 || return false - for i = 1:ng1 + for i in 1:ng1 compare(cmp, getgeom(geo1, i), getgeom(geo2, i), ctx) || return false end end @@ -443,7 +450,7 @@ function compare_coord_seqs(cmp, geo1, geo2, ctx) np1 == np2 || return false coords1 = Vector{Float64}(undef, ncoords1) coords2 = Vector{Float64}(undef, ncoords1) - for i = 1:np1 + for i in 1:np1 coordinates!(coords1, geo1, i, ctx) coordinates!(coords2, geo2, i, ctx) cmp(coords1, coords2) || return false @@ -468,7 +475,7 @@ function compare_coord_seqs(cmp::IsApprox, geo1, geo2, ctx) s1 = 0.0 s2 = 0.0 s12 = 0.0 - for i = 1:np1 + for i in 1:np1 coordinates!(coords1, geo1, i, ctx) coordinates!(coords2, geo2, i, ctx) if ncoords1 == 2 @@ -501,7 +508,7 @@ function Base.hash(geo::AbstractGeometry, h::UInt)::UInt if has_coord_seq(geo) return hash_coord_seq(geo, h) else - for i = 1:ngeom(geo) + for i in 1:ngeom(geo) h = hash(getgeom(geo, i), h) end end @@ -514,14 +521,13 @@ function hash_coord_seq(geo::HasCoordSeq, h::UInt)::UInt end buf = Vector{Float64}(undef, nc) ctx = get_context(geo) - for i = 1:npoints(geo) + for i in 1:npoints(geo) coordinates!(buf, geo, i, ctx) h = hash(buf, h) end return h end - # teach ccall how to get the pointer to pass to libgeos # this way the Julia compiler will track the lifetimes for us Base.unsafe_convert(::Type{Ptr{Cvoid}}, x::AbstractGeometry) = x.ptr diff --git a/test/test_geo_interface.jl b/test/test_geo_interface.jl index 3d65bd0..9f00a9d 100644 --- a/test/test_geo_interface.jl +++ b/test/test_geo_interface.jl @@ -11,7 +11,7 @@ const LG = LibGEOS @test GeoInterface.ncoord(pt) == 2 @test GeoInterface.getcoord(pt, 1) ≈ 1.0 @test GeoInterface.testgeometry(pt) - @test GeoInterface.extent(pt) == Extent(X = (1.0, 1.0), Y = (2.0, 2.0)) + @test GeoInterface.extent(pt) == Extent(X=(1.0, 1.0), Y=(2.0, 2.0)) plot(pt) pt = LibGEOS.Point(1.0, 2.0, 3.0) @@ -32,6 +32,11 @@ const LG = LibGEOS @test GeoInterface.geomtrait(pt) == PointTrait() @test GeoInterface.testgeometry(pt) + @inferred GeoInterface.ncoord(pt) + @inferred GeoInterface.ngeom(pt) + @inferred GeoInterface.getgeom(pt) + @inferred GeoInterface.coordinates(pt) + pt = LibGEOS.readgeom("POINT EMPTY") @test GeoInterface.coordinates(pt) ≈ Float64[] atol = 1e-5 @test GeoInterface.geomtrait(pt) == PointTrait() @@ -49,6 +54,11 @@ const LG = LibGEOS @test GeoInterface.testgeometry(mpt) plot(mpt) + @inferred GeoInterface.ncoord(mpt) + @inferred GeoInterface.ngeom(mpt) + @inferred GeoInterface.getgeom(mpt) + @inferred GeoInterface.coordinates(mpt) + coords = Vector{Float64}[[8, 1], [9, 1], [9, 2], [8, 2]] ls = LibGEOS.LineString(coords) @test GeoInterface.coordinates(ls) == coords @@ -60,6 +70,11 @@ const LG = LibGEOS @test GeoInterface.testgeometry(ls) plot(ls) + @inferred GeoInterface.ncoord(ls) + @inferred GeoInterface.ngeom(ls) + @inferred GeoInterface.getgeom(ls) + @inferred GeoInterface.coordinates(ls) + ls = LibGEOS.readgeom("LINESTRING EMPTY") @test GeoInterface.coordinates(ls) == [] @test GeoInterface.geomtrait(ls) == LineStringTrait() @@ -72,6 +87,11 @@ const LG = LibGEOS @test GeoInterface.testgeometry(mls) plot(mls) + @inferred GeoInterface.ncoord(mls) + @inferred GeoInterface.ngeom(mls) + @inferred GeoInterface.getgeom(mls) + @inferred GeoInterface.coordinates(mls) + coords = Vector{Float64}[[8, 1], [9, 1], [9, 2], [8, 2], [8, 1]] lr = LibGEOS.LinearRing(coords) @test GeoInterface.coordinates(lr) == coords @@ -84,6 +104,11 @@ const LG = LibGEOS # Cannot convert LinearRingTrait to series data for plotting # plot(lr) + @inferred GeoInterface.ncoord(lr) + @inferred GeoInterface.ngeom(lr) + @inferred GeoInterface.getgeom(lr) + @inferred GeoInterface.coordinates(lr) + coords = Vector{Vector{Float64}}[ Vector{Float64}[[0, 0], [10, 0], [10, 10], [0, 10], [0, 0]], Vector{Float64}[[1, 8], [2, 8], [2, 9], [1, 9], [1, 8]], @@ -99,6 +124,11 @@ const LG = LibGEOS @test GeoInterface.testgeometry(polygon) plot(polygon) + @inferred GeoInterface.ncoord(polygon) + @inferred GeoInterface.ngeom(polygon) + @inferred GeoInterface.getgeom(polygon) + @inferred GeoInterface.coordinates(polygon) + polygon = LibGEOS.readgeom("POLYGON EMPTY") @test GeoInterface.coordinates(polygon) == [[]] @test GeoInterface.geomtrait(polygon) == PolygonTrait() @@ -115,13 +145,18 @@ const LG = LibGEOS ]]] @test GeoInterface.geomtrait(multipolygon) == MultiPolygonTrait() @test GeoInterface.testgeometry(multipolygon) - @test GeoInterface.extent(multipolygon) == Extent(X = (0.0, 10.0), Y = (0.0, 10.0)) + @test GeoInterface.extent(multipolygon) == Extent(X=(0.0, 10.0), Y=(0.0, 10.0)) plot(multipolygon) + @inferred GeoInterface.ncoord(multipolygon) + @inferred GeoInterface.ngeom(multipolygon) + @inferred GeoInterface.getgeom(multipolygon) + @inferred GeoInterface.coordinates(multipolygon) + pmultipolygon = LibGEOS.prepareGeom(multipolygon) @test GeoInterface.geomtrait(pmultipolygon) == MultiPolygonTrait() @test GeoInterface.testgeometry(pmultipolygon) - @test GeoInterface.extent(pmultipolygon) == Extent(X = (0.0, 10.0), Y = (0.0, 10.0)) + @test GeoInterface.extent(pmultipolygon) == Extent(X=(0.0, 10.0), Y=(0.0, 10.0)) LibGEOS.destroyGeom(pmultipolygon) geomcollection = LibGEOS.readgeom( @@ -196,6 +231,11 @@ const LG = LibGEOS @test GeoInterface.testgeometry(geomcollection) plot(geomcollection) + @inferred GeoInterface.ncoord(geomcollection) + @inferred GeoInterface.ngeom(geomcollection) + @inferred GeoInterface.getgeom(geomcollection) + # @inferred GeoInterface.coordinates(geomcollection) # can't be inferred + geomcollection = LibGEOS.readgeom( "GEOMETRYCOLLECTION(MULTIPOINT(0 0, 0 0, 1 1),LINESTRING(1 1, 2 2, 2 2, 0 0),POLYGON((5 5, 0 0, 0 2, 2 2, 5 5)))", )