Implemented GeoInterface.

.github/workflows/ci.yml

@@ -3,7 +3,7 @@ on:
   push:
     branches:
       - master
-    tags: '*'
+    tags: "*"
   pull_request:
     branches:
       - master
@@ -15,8 +15,8 @@ jobs:
       fail-fast: false
       matrix:
         version:
-          - '1.3'
-          - '1.7'
+          - "1.6"
+          - "1"
         os:
           - ubuntu-latest
           - macOS-latest
@@ -54,7 +54,7 @@ jobs:
       - uses: actions/checkout@v2
       - uses: julia-actions/setup-julia@v1
         with:
-          version: '1.7'
+          version: "1.7"
       - run: |
           julia --project=docs -e '
             using Pkg

Project.toml

@@ -1,10 +1,11 @@
 name = "GeometryBasics"
 uuid = "5c1252a2-5f33-56bf-86c9-59e7332b4326"
 authors = ["SimonDanisch <[email protected]>"]
-version = "0.4.2"
+version = "0.4.3"
 
 [deps]
 EarCut_jll = "5ae413db-bbd1-5e63-b57d-d24a61df00f5"
+GeoInterface = "cf35fbd7-0cd7-5166-be24-54bfbe79505f"
 IterTools = "c8e1da08-722c-5040-9ed9-7db0dc04731e"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 StaticArrays = "90137ffa-7385-5640-81b9-e52037218182"
@@ -13,17 +14,19 @@ Tables = "bd369af6-aec1-5ad0-b16a-f7cc5008161c"
 
 [compat]
 EarCut_jll = "2"
+GeoInterface = "1.0.1"
 IterTools = "1.3.0"
 StaticArrays = "0.12, 1.0"
 StructArrays = "0.6"
 Tables = "0.2, 1"
-julia = "1.3"
+julia = "1.6"
 
 [extras]
 Aqua = "4c88cf16-eb10-579e-8560-4a9242c79595"
+GeoJSON = "61d90e0f-e114-555e-ac52-39dfb47a3ef9"
 OffsetArrays = "6fe1bfb0-de20-5000-8ca7-80f57d26f881"
 Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 
 [targets]
-test = ["Aqua", "Test", "Random", "OffsetArrays"]
+test = ["Aqua", "GeoJSON", "Test", "Random", "OffsetArrays"]

src/GeometryBasics.jl

@@ -1,6 +1,7 @@
 module GeometryBasics
 
 using StaticArrays, Tables, StructArrays, IterTools, LinearAlgebra
+using GeoInterface
 using EarCut_jll
 
 using Base: @propagate_inbounds
@@ -25,6 +26,7 @@ include("lines.jl")
 include("boundingboxes.jl")
 
 include("deprecated.jl")
+include("geointerface.jl")
 
 export AbstractGeometry, GeometryPrimitive
 export Mat, Point, Vec

src/basic_types.jl

@@ -57,7 +57,6 @@ Fixed Size Polygon, e.g.
 - ...
 """
 struct Ngon{Dim,T<:Real,N,Point<:AbstractPoint{Dim,T}} <: AbstractPolygon{Dim,T}
-
     points::SVector{N,Point}
 end
 
@@ -138,7 +137,6 @@ It applies to infinite dimensions. The structure of this type is designed
 to allow embedding in higher-order spaces by parameterizing on `T`.
 """
 struct Simplex{Dim,T<:Real,N,Point<:AbstractPoint{Dim,T}} <: Polytope{Dim,T}
-
     points::SVector{N,Point}
 end
 
@@ -330,7 +328,6 @@ Base.size(mp::MultiPolygon) = size(mp.polygons)
 
 struct MultiLineString{Dim,T<:Real,Element<:LineString{Dim,T},A<:AbstractVector{Element}} <:
        AbstractVector{Element}
-
     linestrings::A
 end
 
@@ -349,7 +346,6 @@ A collection of points
 """
 struct MultiPoint{Dim,T<:Real,P<:AbstractPoint{Dim,T},A<:AbstractVector{P}} <:
        AbstractVector{P}
-
     points::A
 end
 

src/geointerface.jl

@@ -0,0 +1,115 @@
+# Implementation of trait based interface from https://github.com/JuliaGeo/GeoInterface.jl/
+
+GeoInterface.isgeometry(::Type{<:AbstractGeometry}) = true
+GeoInterface.isgeometry(::Type{<:AbstractFace}) = true
+GeoInterface.isgeometry(::Type{<:AbstractPoint}) = true
+GeoInterface.isgeometry(::Type{<:AbstractVector{<:AbstractGeometry}}) = true
+GeoInterface.isgeometry(::Type{<:AbstractVector{<:AbstractPoint}}) = true
+GeoInterface.isgeometry(::Type{<:AbstractVector{<:LineString}}) = true
+GeoInterface.isgeometry(::Type{<:AbstractVector{<:AbstractPolygon}}) = true
+GeoInterface.isgeometry(::Type{<:AbstractVector{<:AbstractFace}}) = true
+GeoInterface.isgeometry(::Type{<:Mesh}) = true
+
+GeoInterface.geomtrait(::Point) = PointTrait()
+GeoInterface.geomtrait(::Line) = LineTrait()
+GeoInterface.geomtrait(::LineString) = LineStringTrait()
+GeoInterface.geomtrait(::Polygon) = PolygonTrait()
+GeoInterface.geomtrait(::MultiPoint) = MultiPointTrait()
+GeoInterface.geomtrait(::MultiLineString) = MultiLineStringTrait()
+GeoInterface.geomtrait(::MultiPolygon) = MultiPolygonTrait()
+GeoInterface.geomtrait(::Ngon) = PolygonTrait()
+GeoInterface.geomtrait(::AbstractMesh) = PolyhedralSurfaceTrait()
+
+GeoInterface.geomtrait(::Simplex{Dim,T,1}) where {Dim,T} = PointTrait()
+GeoInterface.geomtrait(::Simplex{Dim,T,2}) where {Dim,T} = LineStringTrait()
+GeoInterface.geomtrait(::Simplex{Dim,T,3}) where {Dim,T} = PolygonTrait()
+
+GeoInterface.ncoord(::PointTrait, g::Point) = length(g)
+GeoInterface.getcoord(::PointTrait, g::Point, i::Int) = g[i]
+
+GeoInterface.ngeom(::LineTrait, g::Line) = length(g)
+GeoInterface.getgeom(::LineTrait, g::Line, i::Int) = g[i]
+
+GeoInterface.ngeom(::LineStringTrait, g::LineString) = length(g) + 1  # n line segments + 1
+function GeoInterface.getgeom(::LineStringTrait, g::LineString, i::Int)
+    return GeometryBasics.coordinates(g)[i]
+end
+
+GeoInterface.ngeom(::PolygonTrait, g::Polygon) = length(g.interiors) + 1  # +1 for exterior
+function GeoInterface.getgeom(::PolygonTrait,
+                              g::Polygon,
+                              i::Int)::typeof(g.exterior)
+    return i > 1 ? g.interiors[i - 1] : g.exterior
+end
+
+GeoInterface.ngeom(::MultiPointTrait, g::MultiPoint) = length(g)
+GeoInterface.getgeom(::MultiPointTrait, g::MultiPoint, i::Int) = g[i]
+
+function GeoInterface.ngeom(::MultiLineStringTrait, g::MultiLineString)
+    return length(g)
+end
+function GeoInterface.getgeom(::MultiLineStringTrait, g::MultiLineString,
+                              i::Int)
+    return g[i]
+end
+
+GeoInterface.ngeom(::MultiPolygonTrait, g::MultiPolygon) = length(g)
+GeoInterface.getgeom(::MultiPolygonTrait, g::MultiPolygon, i::Int) = g[i]
+
+function GeoInterface.ncoord(::AbstractGeometryTrait,
+                             ::Simplex{Dim,T,N,P}) where {Dim,T,N,P}
+    return Dim
+end
+function GeoInterface.ncoord(::AbstractGeometryTrait,
+                             ::AbstractGeometry{Dim,T}) where {Dim,T}
+    return Dim
+end
+function GeoInterface.ngeom(::AbstractGeometryTrait,
+                            ::Simplex{Dim,T,N,P}) where {Dim,T,N,P}
+    return N
+end
+GeoInterface.ngeom(::PolygonTrait, ::Ngon) = 1  # can't have any holes
+GeoInterface.getgeom(::PolygonTrait, g::Ngon, _) = LineString(g.points)
+
+function GeoInterface.ncoord(::PolyhedralSurfaceTrait,
+                             ::Mesh{Dim,T,E,V} where {Dim,T,E,V})
+    return Dim
+end
+GeoInterface.ngeom(::PolyhedralSurfaceTrait, g::AbstractMesh) = length(g)
+GeoInterface.getgeom(::PolyhedralSurfaceTrait, g::AbstractMesh, i) = g[i]
+
+function GeoInterface.convert(::Type{Point}, type::PointTrait, geom)
+    dim = Int(ncoord(geom))
+    return Point{dim}(GeoInterface.coordinates(geom))
+end
+
+function GeoInterface.convert(::Type{LineString}, type::LineStringTrait, geom)
+    dim = Int(ncoord(geom))
+    return LineString([Point{dim}(GeoInterface.coordinates(p)) for p in getgeom(geom)])
+end
+
+function GeoInterface.convert(::Type{Polygon}, type::PolygonTrait, geom)
+    t = LineStringTrait()
+    exterior = GeoInterface.convert(LineString, t, GeoInterface.getexterior(geom))
+    if GeoInterface.nhole(geom) == 0
+        return Polygon(exterior)
+    else
+        interiors = GeoInterface.convert.(LineString, Ref(t), GeoInterface.gethole(geom))
+        return Polygon(exterior, interiors)
+    end
+end
+
+function GeoInterface.convert(::Type{MultiPoint}, type::MultiPointTrait, geom)
+    dim = Int(ncoord(geom))
+    return MultiPoint([Point{dim}(GeoInterface.coordinates(p)) for p in getgeom(geom)])
+end
+
+function GeoInterface.convert(::Type{MultiLineString}, type::MultiLineStringTrait, geom)
+    t = LineStringTrait()
+    return MultiLineString([GeoInterface.convert(LineString, t, l) for l in getgeom(geom)])
+end
+
+function GeoInterface.convert(::Type{MultiPolygon}, type::MultiPolygonTrait, geom)
+    t = PolygonTrait()
+    return MultiPolygon([GeoInterface.convert(Polygon, t, poly) for poly in getgeom(geom)])
+end

test/geointerface.jl

@@ -0,0 +1,101 @@
+@testset "Basic types" begin
+    point = Point(2, 3)
+    @test testgeometry(point)
+    @test ncoord(point) == 2
+    @test getcoord(point, 2) == 3
+    @test GeoInterface.coordinates(point) == [2, 3]
+
+    mp = MultiPoint([point, point])
+    @test testgeometry(mp)
+    @test ngeom(mp) == 2
+    @test getgeom(mp, 2) == point
+    @test GeoInterface.coordinates(mp) == [[2, 3], [2, 3]]
+
+    linestring = LineString(Point{2,Int}[(10, 10), (20, 20), (10, 40)])
+    @test testgeometry(linestring)
+    @test ngeom(linestring) == 3
+    @test getgeom(linestring, 1) == Point(10, 10)
+    @test getgeom(linestring, 2) == Point(20, 20)
+    @test getgeom(linestring, 3) == Point(10, 40)
+    @test GeoInterface.coordinates(linestring) == [[10, 10], [20, 20], [10, 40]]
+
+    multilinestring = MultiLineString([linestring, linestring])
+    @test testgeometry(multilinestring)
+    @test GeoInterface.coordinates(multilinestring) ==
+          [[[10, 10], [20, 20], [10, 40]], [[10, 10], [20, 20], [10, 40]]]
+
+    poly = Polygon(rand(Point{2,Float32}, 5), [rand(Point{2,Float32}, 5)])
+    @test testgeometry(poly)
+    @test length(GeoInterface.coordinates(poly)) == 2
+    @test length(GeoInterface.coordinates(poly)[1]) == 5
+
+    triangle = Triangle(point, point, point)
+    @test testgeometry(triangle)
+    @test length(GeoInterface.coordinates(triangle)) == 1
+    @test length(GeoInterface.coordinates(triangle)[1]) == 3
+
+    polys = MultiPolygon([poly, poly])
+    @test testgeometry(polys)
+    @test length(GeoInterface.coordinates(polys)) == 2
+    @test length(GeoInterface.coordinates(polys)[1]) == 2
+    @test length(GeoInterface.coordinates(polys)[1][1]) == 5
+end
+
+@testset "Mesh" begin
+    mesh = triangle_mesh(Sphere(Point3f(0), 1))
+    @test testgeometry(mesh)
+end
+
+@testset "Convert" begin
+    # convert GeoJSON geometry types to GeometryBasics via the GeoInterface
+    point_str = """{"type":"Point","coordinates":[30.1,10.1]}"""
+    point_3d_str = """{"type":"Point","coordinates":[30.1,10.1,5.1]}"""
+    linestring_str = """{"type":"LineString","coordinates":[[30.1,10.1],[10.1,30.1],[40.1,40.1]]}"""
+    polygon_str = """{"type":"Polygon","coordinates":[[[30.1,10.1],[40.1,40.1],[20.1,40.1],[10.1,20.1],[30.1,10.1]]]}"""
+    polygon_hole_str = """{"type":"Polygon","coordinates":[[[35.1,10.1],[45.1,45.1],[15.1,40.1],[10.1,20.1],[35.1,10.1]],[[20.1,30.1],[35.1,35.1],[30.1,20.1],[20.1,30.1]]]}"""
+    multipoint_str = """{"type":"MultiPoint","coordinates":[[10.1,40.1],[40.1,30.1],[20.1,20.1],[30.1,10.1]]}"""
+    multilinestring_str = """{"type":"MultiLineString","coordinates":[[[10.1,10.1],[20.1,20.1],[10.1,40.1]],[[40.1,40.1],[30.1,30.1],[40.1,20.1],[30.1,10.1]]]}"""
+    multipolygon_str = """{"type":"MultiPolygon","coordinates":[[[[30.1,20.1],[45.1,40.1],[10.1,40.1],[30.1,20.1]]],[[[15.1,5.1],[40.1,10.1],[10.1,20.1],[5.1,10.1],[15.1,5.1]]]]}"""
+    multipolygon_hole_str = """{"type":"MultiPolygon","coordinates":[[[[40.1,40.1],[20.1,45.1],[45.1,30.1],[40.1,40.1]]],[[[20.1,35.1],[10.1,30.1],[10.1,10.1],[30.1,5.1],[45.1,20.1],[20.1,35.1]],[[30.1,20.1],[20.1,15.1],[20.1,25.1],[30.1,20.1]]]]}"""
+
+    point_json = GeoJSON.read(point_str)
+    point_3d_json = GeoJSON.read(point_3d_str)
+    linestring_json = GeoJSON.read(linestring_str)
+    polygon_json = GeoJSON.read(polygon_str)
+    polygon_hole_json = GeoJSON.read(polygon_hole_str)
+    multipoint_json = GeoJSON.read(multipoint_str)
+    multilinestring_json = GeoJSON.read(multilinestring_str)
+    multipolygon_json = GeoJSON.read(multipolygon_str)
+    multipolygon_hole_json = GeoJSON.read(multipolygon_hole_str)
+
+    point_gb = GeoInterface.convert(Point, point_json)
+    point_3d_gb = GeoInterface.convert(Point, point_3d_json)
+    linestring_gb = GeoInterface.convert(LineString, linestring_json)
+    polygon_gb = GeoInterface.convert(Polygon, polygon_json)
+    polygon_hole_gb = GeoInterface.convert(Polygon, polygon_hole_json)
+    multipoint_gb = GeoInterface.convert(MultiPoint, multipoint_json)
+    multilinestring_gb = GeoInterface.convert(MultiLineString, multilinestring_json)
+    multipolygon_gb = GeoInterface.convert(MultiPolygon, multipolygon_json)
+    multipolygon_hole_gb = GeoInterface.convert(MultiPolygon, multipolygon_hole_json)
+
+    @test point_gb === Point{2, Float64}(30.1, 10.1)
+    @test point_3d_gb === Point{3, Float64}(30.1, 10.1, 5.1)
+    @test linestring_gb isa LineString
+    @test length(linestring_gb) == 2
+    @test eltype(linestring_gb) == Line{2, Float64}
+    @test polygon_gb isa Polygon
+    @test isempty(polygon_gb.interiors)
+    @test polygon_hole_gb isa Polygon
+    @test length(polygon_hole_gb.interiors) == 1
+    @test multipoint_gb isa MultiPoint
+    @test length(multipoint_gb) == 4
+    @test multipoint_gb[4] === Point{2, Float64}(30.1, 10.1)
+    @test multilinestring_gb isa MultiLineString
+    @test length(multilinestring_gb) == 2
+    @test multipolygon_gb isa MultiPolygon
+    @test length(multipolygon_gb) == 2
+    @test multipolygon_hole_gb isa MultiPolygon
+    @test length(multipolygon_hole_gb) == 2
+    @test length(multipolygon_hole_gb[1].interiors) == 0
+    @test length(multipolygon_hole_gb[2].interiors) == 1
+end

test/runtests.jl

@@ -2,6 +2,8 @@ using Test, Random, StructArrays, Tables, StaticArrays, OffsetArrays
 using GeometryBasics
 using LinearAlgebra
 using GeometryBasics: attributes
+using GeoInterface
+using GeoJSON
 
 @testset "GeometryBasics" begin
 
@@ -715,6 +717,10 @@ end
     include("fixed_arrays.jl")
 end
 
+@testset "GeoInterface" begin
+    include("geointerface.jl")
+end
+
 using Aqua
 # Aqua tests
 # Intervals brings a bunch of ambiquities unfortunately