From 5f6068b59cee0b8408370581c2658c96f64980dd Mon Sep 17 00:00:00 2001 From: rafaqz Date: Wed, 27 Apr 2022 17:39:23 +0400 Subject: [PATCH 01/16] update for the new GeoInterface --- Project.toml | 4 +- src/Shapefile.jl | 106 +++++++++++++++++++++++++++++++------------ src/geo_interface.jl | 23 +++++----- test/runtests.jl | 43 ++++++++++-------- 4 files changed, 115 insertions(+), 61 deletions(-) diff --git a/Project.toml b/Project.toml index aab02c5..24065ac 100644 --- a/Project.toml +++ b/Project.toml @@ -5,6 +5,7 @@ version = "0.7.4" [deps] DBFTables = "75c7ada1-017a-5fb6-b8c7-2125ff2d6c93" +Extents = "411431e0-e8b7-467b-b5e0-f676ba4f2910" GeoFormatTypes = "68eda718-8dee-11e9-39e7-89f7f65f511f" GeoInterface = "cf35fbd7-0cd7-5166-be24-54bfbe79505f" RecipesBase = "3cdcf5f2-1ef4-517c-9805-6587b60abb01" @@ -13,7 +14,8 @@ Tables = "bd369af6-aec1-5ad0-b16a-f7cc5008161c" [compat] DBFTables = "0.2, 1" GeoFormatTypes = "0.4" -GeoInterface = "0.4, 0.5" +GeoInterface = "1.0" +Extents = "0.1" RecipesBase = "1" Tables = "0.2, 1" julia = "1" diff --git a/src/Shapefile.jl b/src/Shapefile.jl index c050c92..70518d1 100644 --- a/src/Shapefile.jl +++ b/src/Shapefile.jl @@ -1,6 +1,6 @@ module Shapefile -import GeoInterface, DBFTables, Tables, GeoFormatTypes +import GeoFormatTypes, GeoInterface, DBFTables, Extents, Tables using RecipesBase @@ -18,6 +18,13 @@ struct Rect top::Float64 end +abstract type AbstractShape end + +isgeometry(::AbstractShape) = true +GeoInterface.ncoord(::AbstractShape) = 2 # With specific methods when 3 + +GeoInterface.@enable_geo_plots AbstractShape + """ Interval @@ -28,20 +35,26 @@ struct Interval right::Float64 end +abstract type AbstractPoint <: AbstractShape end + +GeoInterface.geomtype(point::Type{<:AbstractPoint}) = GeoInterface.PointTrait() +GeoInterface.x(point::AbstractPoint) = point.x +GeoInterface.y(point::AbstractPoint) = point.y + """ - Point <: GeoInterface.AbstractPoint + Point <: AbstractPoint Point from a shape file. Fields `x`, `y` hold the spatial location. """ -struct Point <: GeoInterface.AbstractPoint +struct Point <: AbstractPoint x::Float64 y::Float64 end """ - PointM <: GeoInterface.AbstractPoint + PointM <: AbstractPoint Point from a shape file. @@ -49,14 +62,16 @@ Fields `x`, `y` hold the spatial location. Includes a measure field `m`, holding a value for the point. """ -struct PointM <: GeoInterface.AbstractPoint +struct PointM <: AbstractPoint x::Float64 y::Float64 m::Float64 # measure end +GeoInterface.m(point::PointM) = point.m + """ - PointZ <: GeoInterface.AbstractPoint + PointZ <: AbstractPoint Three dimensional point, from a shape file. @@ -64,15 +79,23 @@ Fields `x`, `y`, `z` hold the spatial location. Includes a measure field `m`, holding a value for the point. """ -struct PointZ <: GeoInterface.AbstractPoint +struct PointZ <: AbstractPoint x::Float64 y::Float64 z::Float64 m::Float64 # measure end +GeoInterface.m(point::PointZ) = point.m +GeoInterface.z(point::PointZ) = point.z +GeoInterface.ncoord(::PointZ) = 3 + +abstract type AbstractPolyline <: AbstractShape end + +GeoInterface.geomtype(::Type{<:AbstractPolyline}) = GeoInterface.MultiLineStringTrait() + """ - Polyline <: GeoInterface.AbstractMultiLineString + Polyline <: AbstractPolyline Represents a single or multiple polylines from a shape file. @@ -81,14 +104,14 @@ Represents a single or multiple polylines from a shape file. - `parts`: a `Vector` of `Int32` indicating the line each point belongs to. - `MBR`: `nothing` or the known bounding box. Can be retrieved with `GeoInterface.bbox`. """ -struct Polyline <: GeoInterface.AbstractMultiLineString +struct Polyline <: AbstractPolyline MBR::Rect parts::Vector{Int32} points::Vector{Point} end """ - PolylineM <: GeoInterface.AbstractMultiLineString + PolylineM <: AbstractPolyline Polyline from a shape file, with measures. @@ -98,7 +121,7 @@ Polyline from a shape file, with measures. - `MBR`: `nothing` or the known bounding box. Can be retrieved with `GeoInterface.bbox`. - `measures`: holds values from each point. """ -struct PolylineM <: GeoInterface.AbstractMultiLineString +struct PolylineM <: AbstractPolyline MBR::Rect parts::Vector{Int32} points::Vector{Point} @@ -106,7 +129,7 @@ struct PolylineM <: GeoInterface.AbstractMultiLineString end """ - PolylineZ <: GeoInterface.AbstractMultiLineString + PolylineZ <: AbstractPolyline Three dimensional polyline of from a shape file. @@ -117,7 +140,7 @@ Three dimensional polyline of from a shape file. - `MBR`: `nothing` or the known bounding box. Can be retrieved with `GeoInterface.bbox`. - `measures`: holds values from each point. """ -struct PolylineZ <: GeoInterface.AbstractMultiLineString +struct PolylineZ <: AbstractPolyline MBR::Rect parts::Vector{Int32} points::Vector{Point} @@ -125,8 +148,14 @@ struct PolylineZ <: GeoInterface.AbstractMultiLineString measures::Vector{Float64} end +GeoInterface.ncoord(::PolylineZ) = 3 + +abstract type AbstractPolygon <: AbstractShape end + +GeoInterface.geomtype(::Type{<:AbstractPolygon}) = GeoInterface.PolygonTrait() + """ - Polygon <: GeoInterface.AbstractMultiPolygon + Polygon <: AbstractPolygon Represents a polygon from a shape file. @@ -135,7 +164,7 @@ Represents a polygon from a shape file. - `parts`: a `Vector` of `Int32` indicating the polygon each point belongs to. - `MBR`: `nothing` or the known bounding box. Can be retrieved with `GeoInterface.bbox`. """ -struct Polygon <: GeoInterface.AbstractMultiPolygon +struct Polygon <: AbstractPolygon MBR::Rect parts::Vector{Int32} points::Vector{Point} @@ -145,7 +174,7 @@ Base.show(io::IO, p::Polygon) = print(io, "Polygon(", length(p.points), " Points)") """ - PolygonM <: GeoInterface.AbstractMultiPolygon + PolygonM <: AbstractPolygon Represents a polygon from a shape file @@ -155,7 +184,7 @@ Represents a polygon from a shape file - `MBR`: `nothing` or the known bounding box. Can be retrieved with `GeoInterface.bbox`. - `measures`: holds values from each point. """ -struct PolygonM <: GeoInterface.AbstractMultiPolygon +struct PolygonM <: AbstractPolygon MBR::Rect parts::Vector{Int32} points::Vector{Point} @@ -163,7 +192,7 @@ struct PolygonM <: GeoInterface.AbstractMultiPolygon end """ - PolygonZ <: GeoInterface.AbstractMultiPolygon + PolygonZ <: AbstractPolygon A three dimensional polygon from a shape file. @@ -174,7 +203,7 @@ A three dimensional polygon from a shape file. - `MBR`: `nothing` or the known bounding box. Can be retrieved with `GeoInterface.bbox`. - `measures`: holds values from each point. """ -struct PolygonZ <: GeoInterface.AbstractMultiPolygon +struct PolygonZ <: AbstractPolygon MBR::Rect parts::Vector{Int32} points::Vector{Point} @@ -182,8 +211,14 @@ struct PolygonZ <: GeoInterface.AbstractMultiPolygon measures::Vector{Float64} end +GeoInterface.ncoord(::PolygonZ) = 3 + +abstract type AbstractMultiPoint <: AbstractShape end + +GeoInterface.geomtype(::Type{<:AbstractMultiPoint}) = GeoInterface.MultiPointTrait() + """ - MultiPoint <: GeoInterface.AbstractMultiPoint + MultiPoint <: AbstractMultiPoint Collection of points, from a shape file. @@ -191,13 +226,13 @@ Collection of points, from a shape file. - `points`: a `Vector` of [`Point`](@ref). - `MBR`: `nothing` or the known bounding box. Can be retrieved with `GeoInterface.bbox`. """ -struct MultiPoint <: GeoInterface.AbstractMultiPoint +struct MultiPoint <: AbstractMultiPoint MBR::Rect points::Vector{Point} end """ - MultiPointM <: GeoInterface.AbstractMultiPoint + MultiPointM <: AbstractMultiPoint Collection of points, from a shape file. @@ -210,14 +245,20 @@ May have a known bounding box, which can be retrieved with `GeoInterface.bbox`. - `MBR`: `nothing` or the known bounding box. Can be retrieved with `GeoInterface.bbox`. - `measures`: holds values from each point. """ -struct MultiPointM <: GeoInterface.AbstractMultiPoint +struct MultiPointM <: AbstractMultiPoint MBR::Rect points::Vector{Point} measures::Vector{Float64} end """ - MultiPointZ <: GeoInterface.AbstractMultiPoint + MultiPointZ <: AbstractMultiPoint + +Collection of 3d points, from a shape file. + +Includes a `measures` field, holding values from each point. + +May have a known bounding box, which can be retrieved with `GeoInterface.bbox`. # Fields - `points`: a `Vector` of [`Point`](@ref). @@ -225,23 +266,28 @@ end - `MBR`: `nothing` or the known bounding box. Can be retrieved with `GeoInterface.bbox`. - `measures`: holds values from each point. """ -struct MultiPointZ <: GeoInterface.AbstractMultiPoint +struct MultiPointZ <: AbstractMultiPoint MBR::Rect points::Vector{Point} zvalues::Vector{Float64} measures::Vector{Float64} end +GeoInterface.ncoord(::MultiPointZ) = 3 + """ - MultiPatch <: GeoInterface.AbstractGeometry + MultiPatch + +Stores a collection of patch representing the boundary of a 3d object. # Fields - `points`: a `Vector` of [`Point`](@ref) represents a one or multiple spatial objects. - `parts`: a `Vector` of `Int32` indicating the object each point belongs to. - `parttypes`: a `Vector` of `Int32` indicating the type of object each point belongs to. +- `zvalues`: a `Vector` of `Float64` indicating absolute or relative heights. - `MBR`: `nothing` or the known bounding box. Can be retrieved with `GeoInterface.bbox`. """ -struct MultiPatch <: GeoInterface.AbstractGeometry +struct MultiPatch <: AbstractShape MBR::Rect parts::Vector{Int32} parttypes::Vector{Int32} @@ -250,6 +296,10 @@ struct MultiPatch <: GeoInterface.AbstractGeometry # measures::Vector{Float64} # (optional) end +GeoInterface.geomtype(::Type{<:MultiPatch}) = GeoInterface.MultiPointTrait() + +GeoInterface.ncoord(::MultiPatch) = 3 + const SHAPETYPE = Dict{Int32,DataType}( 0 => Missing, 1 => Point, @@ -279,7 +329,7 @@ The Vector of shape object can be accessed with `shapes(handle)`. `Handle` may have a known bounding box, which can be retrieved with `GeoInterface.bbox`. """ -mutable struct Handle{T<:Union{<:GeoInterface.AbstractGeometry,Missing}} +mutable struct Handle{T<:Union{<:AbstractShape,Missing}} code::Int32 length::Int32 version::Int32 diff --git a/src/geo_interface.jl b/src/geo_interface.jl index 9fedc27..4463b37 100644 --- a/src/geo_interface.jl +++ b/src/geo_interface.jl @@ -1,17 +1,17 @@ -GeoInterface.coordinates(obj::Union{Point,PointM}) = Float64[obj.x, obj.y] -GeoInterface.coordinates(obj::PointZ) = Float64[obj.x, obj.y, obj.z] +GeoInterface.coordinates(::GeoInterface.PointTrait, obj::Union{Point,PointM}) = [obj.x, obj.y] +GeoInterface.coordinates(::GeoInterface.PointTrait, obj::PointZ) = [obj.x, obj.y, obj.z] GeoInterface.coordinates(obj::Union{MultiPoint,MultiPointM,MultiPointZ}) = Vector{Float64}[GeoInterface.coordinates(p) for p in obj.points] -function GeoInterface.coordinates(obj::Union{Polyline,PolylineM,PolylineZ}) +function GeoInterface.coordinates(::GeoInterface.MultiLineStringTrait, obj::Union{Polyline,PolylineM,PolylineZ}) npoints = length(obj.points) nparts = length(obj.parts) @assert obj.parts[nparts] <= npoints push!(obj.parts, npoints) coords = Vector{Vector{Float64}}[Vector{Float64}[GeoInterface.coordinates(obj.points[j+1]) for j = obj.parts[i]:(obj.parts[i+1]-1)] for i = 1:nparts] pop!(obj.parts) - coords + return coords end # Only supports 2D geometries for now @@ -26,11 +26,11 @@ function inring(pt::Vector{Float64}, ring::Vector{Vector{Float64}}) for k = 2:length(ring) isinside = intersect(ring[k], ring[k-1]) ? !isinside : isinside end - isinside + return isinside end # ported from https://github.com/Esri/terraformer-arcgis-parser/blob/master/terraformer-arcgis-parser.js#L168-L253 -function GeoInterface.coordinates(obj::Union{Polygon,PolygonZ,PolygonM}) +function GeoInterface.coordinates(::GeoInterface.PolygonTrait, obj::Union{Polygon,PolygonZ,PolygonM}) npoints = length(obj.points) nparts = length(obj.parts) @assert obj.parts[end] <= npoints @@ -68,17 +68,16 @@ function GeoInterface.coordinates(obj::Union{Polygon,PolygonZ,PolygonM}) end end if !found - push!(coords, Vector{Vector{Float64}}[hole]) # hole is not inside any ring; make it a polygon + push!(coords, [hole]) # hole is not inside any ring; make it a polygon end end - coords + return coords end # bbox -const HasMBR = Union{Polyline,PolylineM, PolylineZ, Polygon, PolygonM, PolygonZ, - MultiPoint, MultiPointM, MultiPointZ, MultiPatch, Handle} +const HasMBR = Union{AbstractPolyline,AbstractPolygon,AbstractMultiPoint,MultiPatch,Handle} -function GeoInterface.bbox(x::HasMBR) +function GeoInterface.extent(x::HasMBR) rect = x.MBR - return rect.left, rect.bottom, rect.right, rect.top + return Extents.Extent(X=(rect.left, rect.right), Y=(rect.bottom, rect.top)) end diff --git a/test/runtests.jl b/test/runtests.jl index 0dc75b0..9dcb85a 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -1,82 +1,84 @@ -using Shapefile, GeoInterface, Plots +using Shapefile, GeoInterface, Plots, Extents using Test +plot(Shapefile.Point(1, 2)) + test_tuples = [ ( - path=joinpath("shapelib_testcases", "test.shp"), + path=joinpath("shapelib_testcases", "test.shp"), geomtype=Shapefile.Polygon, coordinates=Array{Array{Array{Array{Float64,1},1},1},1}[Array{Array{Array{Float64,1},1},1}[Array{Array{Float64,1},1}[Array{Float64,1}[[20.0,20.0],[20.0,30.0],[30.0,30.0],[20.0,20.0]]],Array{Array{Float64,1},1}[Array{Float64,1}[[0.0,0.0],[100.0,0.0],[100.0,100.0],[0.0,100.0],[0.0,0.0]]]],Array{Array{Array{Float64,1},1},1}[Array{Array{Float64,1},1}[Array{Float64,1}[[150.0,150.0],[160.0,150.0],[180.0,170.0],[150.0,150.0]]]],Array{Array{Array{Float64,1},1},1}[Array{Array{Float64,1},1}[Array{Float64,1}[[150.0,150.0],[160.0,150.0],[180.0,170.0],[150.0,150.0]]]]], - bbox=(0.0, 0.0, 180.0, 170.0), + extent=Extent(X=(0.0, 180.0), Y=(0.0, 170.0)), ),( path=joinpath("shapelib_testcases", "test0.shp"), geomtype=Missing, coordinates=nothing, - bbox=(0.0, 0.0, 10.0, 20.0), + extent=Extent(X=(0.0, 10.0), Y=(0.0, 20.0)), ),( path=joinpath("shapelib_testcases", "test1.shp"), geomtype=Shapefile.Point, coordinates=Array{Float64,1}[[1.0,2.0],[10.0,20.0]], - bbox=(1.0, 2.0, 10.0, 20.0), + extent=Extent(X=(1.0, 10.0), Y=(2.0, 20.0)), ),( path=joinpath("shapelib_testcases", "test2.shp"), geomtype=Shapefile.PointZ, coordinates=Array{Float64,1}[[1.0,2.0,3.0],[10.0,20.0,30.0]], - bbox=(1.0, 2.0, 10.0, 20.0), + extent=Extent(X=(1.0, 10.0), Y=(2.0, 20.0)), ),( path=joinpath("shapelib_testcases", "test3.shp"), geomtype=Shapefile.PointM, coordinates=Array{Float64,1}[[1.0,2.0],[10.0,20.0]], - bbox=(1.0, 2.0, 10.0, 20.0), + extent=Extent(X=(1.0, 10.0), Y=(2.0, 20.0)), ),( path=joinpath("shapelib_testcases", "test4.shp"), geomtype=Shapefile.MultiPoint, coordinates=Array{Array{Float64,1},1}[Array{Float64,1}[[1.15,2.25],[2.15,3.25],[3.15,4.25],[4.15,5.25]],Array{Float64,1}[[11.15,12.25],[12.15,13.25],[13.15,14.25],[14.15,15.25]],Array{Float64,1}[[21.15,22.25],[22.15,23.25],[23.15,24.25],[24.15,25.25]]], - bbox=(1.15, 2.25, 24.15, 25.25), + extent=Extent(X=(1.15, 24.15), Y=(2.25, 25.25)), ),( path=joinpath("shapelib_testcases", "test5.shp"), geomtype=Shapefile.MultiPointZ, coordinates=Array{Array{Float64,1},1}[Array{Float64,1}[[1.15,2.25],[2.15,3.25],[3.15,4.25],[4.15,5.25]],Array{Float64,1}[[11.15,12.25],[12.15,13.25],[13.15,14.25],[14.15,15.25]],Array{Float64,1}[[21.15,22.25],[22.15,23.25],[23.15,24.25],[24.15,25.25]]], - bbox=(1.15, 2.25, 24.15, 25.25), + extent=Extent(X=(1.15, 24.15), Y=(2.25, 25.25)), ),( path=joinpath("shapelib_testcases", "test6.shp"), geomtype=Shapefile.MultiPointM, coordinates=Array{Array{Float64,1},1}[Array{Float64,1}[[1.15,2.25],[2.15,3.25],[3.15,4.25],[4.15,5.25]],Array{Float64,1}[[11.15,12.25],[12.15,13.25],[13.15,14.25],[14.15,15.25]],Array{Float64,1}[[21.15,22.25],[22.15,23.25],[23.15,24.25],[24.15,25.25]]], - bbox=(1.15, 2.25, 24.15, 25.25), + extent=Extent(X=(1.15, 24.15), Y=(2.25, 25.25)), ),( path=joinpath("shapelib_testcases", "test7.shp"), geomtype=Shapefile.Polyline, coordinates=Array{Array{Array{Float64,1},1},1}[Array{Array{Float64,1},1}[Array{Float64,1}[[1.0,1.0],[2.0,1.0],[2.0,2.0],[1.0,2.0],[1.0,1.0]]],Array{Array{Float64,1},1}[Array{Float64,1}[[1.0,4.0],[2.0,4.0],[2.0,5.0],[1.0,5.0],[1.0,4.0]]],Array{Array{Float64,1},1}[Array{Float64,1}[[1.0,7.0],[2.0,7.0],[2.0,8.0],[1.0,8.0],[1.0,7.0]]],Array{Array{Float64,1},1}[Array{Float64,1}[[0.0,0.0],[0.0,100.0],[100.0,100.0],[100.0,0.0],[0.0,0.0]],Array{Float64,1}[[10.0,20.0],[30.0,20.0],[30.0,40.0],[10.0,40.0],[10.0,20.0]],Array{Float64,1}[[60.0,20.0],[90.0,20.0],[90.0,40.0],[60.0,40.0],[60.0,20.0]]]], - bbox=(0.0, 0.0, 100.0, 100.0), + extent=Extent(X=(0.0, 100.0), Y=(0.0, 100.0)), ),( path=joinpath("shapelib_testcases", "test8.shp"), geomtype=Shapefile.PolylineZ, coordinates=Array{Array{Array{Float64,1},1},1}[Array{Array{Float64,1},1}[Array{Float64,1}[[1.0,1.0],[2.0,1.0],[2.0,2.0],[1.0,2.0],[1.0,1.0]]],Array{Array{Float64,1},1}[Array{Float64,1}[[1.0,4.0],[2.0,4.0],[2.0,5.0],[1.0,5.0],[1.0,4.0]]],Array{Array{Float64,1},1}[Array{Float64,1}[[1.0,7.0],[2.0,7.0],[2.0,8.0],[1.0,8.0],[1.0,7.0]]],Array{Array{Float64,1},1}[Array{Float64,1}[[0.0,0.0],[0.0,100.0],[100.0,100.0],[100.0,0.0],[0.0,0.0]],Array{Float64,1}[[10.0,20.0],[30.0,20.0],[30.0,40.0],[10.0,40.0],[10.0,20.0]],Array{Float64,1}[[60.0,20.0],[90.0,20.0],[90.0,40.0],[60.0,40.0],[60.0,20.0]]]], - bbox=(0.0, 0.0, 100.0, 100.0), + extent=Extent(X=(0.0, 100.0), Y=(0.0, 100.0)), ),( path=joinpath("shapelib_testcases", "test9.shp"), geomtype=Shapefile.PolylineM, coordinates=Array{Array{Array{Float64,1},1},1}[Array{Array{Float64,1},1}[Array{Float64,1}[[1.0,1.0],[2.0,1.0],[2.0,2.0],[1.0,2.0],[1.0,1.0]]],Array{Array{Float64,1},1}[Array{Float64,1}[[1.0,4.0],[2.0,4.0],[2.0,5.0],[1.0,5.0],[1.0,4.0]]],Array{Array{Float64,1},1}[Array{Float64,1}[[1.0,7.0],[2.0,7.0],[2.0,8.0],[1.0,8.0],[1.0,7.0]]],Array{Array{Float64,1},1}[Array{Float64,1}[[0.0,0.0],[0.0,100.0],[100.0,100.0],[100.0,0.0],[0.0,0.0]],Array{Float64,1}[[10.0,20.0],[30.0,20.0],[30.0,40.0],[10.0,40.0],[10.0,20.0]],Array{Float64,1}[[60.0,20.0],[90.0,20.0],[90.0,40.0],[60.0,40.0],[60.0,20.0]]]], - bbox=(0.0, 0.0, 100.0, 100.0), + extent=Extent(X=(0.0, 100.0), Y=(0.0, 100.0)), ),( path=joinpath("shapelib_testcases", "test10.shp"), geomtype=Shapefile.Polygon, coordinates=[[[[[1.0, 1.0], [2.0, 1.0], [2.0, 2.0], [1.0, 2.0], [1.0, 1.0]]]], [[[[1.0, 4.0], [2.0, 4.0], [2.0, 5.0], [1.0, 5.0], [1.0, 4.0]]]], [[[[1.0, 7.0], [2.0, 7.0], [2.0, 8.0], [1.0, 8.0], [1.0, 7.0]]]], [[[[0.0, 0.0], [0.0, 100.0], [100.0, 100.0], [100.0, 0.0], [0.0, 0.0]], [[10.0, 20.0], [30.0, 20.0], [30.0, 40.0], [10.0, 40.0], [10.0, 20.0]], [[60.0, 20.0], [90.0, 20.0], [90.0, 40.0], [60.0, 40.0], [60.0, 20.0]]]]], - bbox=(0.0, 0.0, 100.0, 100.0), + extent=Extent(X=(0.0, 100.0), Y=(0.0, 100.0)), ),( path=joinpath("shapelib_testcases", "test11.shp"), geomtype=Shapefile.PolygonZ, coordinates=[[[[[1.0, 1.0], [2.0, 1.0], [2.0, 2.0], [1.0, 2.0], [1.0, 1.0]]]], [[[[1.0, 4.0], [2.0, 4.0], [2.0, 5.0], [1.0, 5.0], [1.0, 4.0]]]], [[[[1.0, 7.0], [2.0, 7.0], [2.0, 8.0], [1.0, 8.0], [1.0, 7.0]]]], [[[[0.0, 0.0], [0.0, 100.0], [100.0, 100.0], [100.0, 0.0], [0.0, 0.0]], [[10.0, 20.0], [30.0, 20.0], [30.0, 40.0], [10.0, 40.0], [10.0, 20.0]], [[60.0, 20.0], [90.0, 20.0], [90.0, 40.0], [60.0, 40.0], [60.0, 20.0]]]]], - bbox=(0.0, 0.0, 100.0, 100.0), + extent=Extent(X=(0.0, 100.0), Y=(0.0, 100.0)), ),( path=joinpath("shapelib_testcases", "test12.shp"), geomtype=Shapefile.PolygonM, coordinates=[[[[[1.0, 1.0], [2.0, 1.0], [2.0, 2.0], [1.0, 2.0], [1.0, 1.0]]]], [[[[1.0, 4.0], [2.0, 4.0], [2.0, 5.0], [1.0, 5.0], [1.0, 4.0]]]], [[[[1.0, 7.0], [2.0, 7.0], [2.0, 8.0], [1.0, 8.0], [1.0, 7.0]]]], [[[[0.0, 0.0], [0.0, 100.0], [100.0, 100.0], [100.0, 0.0], [0.0, 0.0]], [[10.0, 20.0], [30.0, 20.0], [30.0, 40.0], [10.0, 40.0], [10.0, 20.0]], [[60.0, 20.0], [90.0, 20.0], [90.0, 40.0], [60.0, 40.0], [60.0, 20.0]]]]], - bbox=(0.0, 0.0, 100.0, 100.0), + extent=Extent(X=(0.0, 100.0), Y=(0.0, 100.0)), ),( path=joinpath("shapelib_testcases/test13.shp"), geomtype=Shapefile.MultiPatch, coordinates=nothing, - bbox=(0.0, 0.0, 100.0, 100.0), + extent=Extent(X=(0.0, 100.0), Y=(0.0, 100.0)), ) ] @@ -105,9 +107,9 @@ for test in test_tuples if !(test.geomtype <: Union{Missing,Shapefile.MultiPatch}) @test GeoInterface.coordinates.(shp.shapes) == test.coordinates end - @test shp.MBR == Shapefile.Rect(test.bbox...) - @test GeoInterface.bbox(shp) == test.bbox - + ext = test.extent + @test shp.MBR == Shapefile.Rect(ext.X[1], ext.Y[1], ext.X[2], ext.Y[2]) + @test GeoInterface.extent(shp) == test.extent # Multipatch can't be plotted, but it's obscure anyway if !(test.geomtype == Shapefile.MultiPatch) plot(shp) # Just test that it actually plots @@ -117,6 +119,7 @@ end # Test all .shx files; the values in .shx must match the .shp offsets +test = test_tuples[1] for test in test_tuples offsets = Int32[] From d26bcb59db0e0702a564b452cf9dde6fad43802b Mon Sep 17 00:00:00 2001 From: rafaqz Date: Mon, 9 May 2022 20:05:09 +0400 Subject: [PATCH 02/16] huge refactor and test overhaul --- Project.toml | 2 + src/Shapefile.jl | 587 +------------------------------------------ src/common.jl | 46 ++++ src/extent.jl | 14 ++ src/geo_interface.jl | 83 ------ src/handle.jl | 94 +++++++ src/multipatch.jl | 44 ++++ src/multipoints.jl | 107 ++++++++ src/plotrecipes.jl | 2 + src/points.jl | 84 +++++++ src/polygons.jl | 233 +++++++++++++++++ src/polylines.jl | 147 +++++++++++ src/table.jl | 1 + test/runtests.jl | 66 ++--- 14 files changed, 824 insertions(+), 686 deletions(-) create mode 100644 src/common.jl create mode 100644 src/extent.jl delete mode 100644 src/geo_interface.jl create mode 100644 src/handle.jl create mode 100644 src/multipatch.jl create mode 100644 src/multipoints.jl create mode 100644 src/points.jl create mode 100644 src/polygons.jl create mode 100644 src/polylines.jl diff --git a/Project.toml b/Project.toml index 24065ac..b5382f7 100644 --- a/Project.toml +++ b/Project.toml @@ -8,6 +8,7 @@ DBFTables = "75c7ada1-017a-5fb6-b8c7-2125ff2d6c93" Extents = "411431e0-e8b7-467b-b5e0-f676ba4f2910" GeoFormatTypes = "68eda718-8dee-11e9-39e7-89f7f65f511f" GeoInterface = "cf35fbd7-0cd7-5166-be24-54bfbe79505f" +GeoInterfaceRecipes = "0329782f-3d07-4b52-b9f6-d3137cf03c7a" RecipesBase = "3cdcf5f2-1ef4-517c-9805-6587b60abb01" Tables = "bd369af6-aec1-5ad0-b16a-f7cc5008161c" @@ -15,6 +16,7 @@ Tables = "bd369af6-aec1-5ad0-b16a-f7cc5008161c" DBFTables = "0.2, 1" GeoFormatTypes = "0.4" GeoInterface = "1.0" +GeoInterfaceRecipes = "1.0" Extents = "0.1" RecipesBase = "1" Tables = "0.2, 1" diff --git a/src/Shapefile.jl b/src/Shapefile.jl index 70518d1..06fea61 100644 --- a/src/Shapefile.jl +++ b/src/Shapefile.jl @@ -1,304 +1,19 @@ module Shapefile -import GeoFormatTypes, GeoInterface, DBFTables, Extents, Tables +import GeoFormatTypes, GeoInterface, GeoInterfaceRecipes, DBFTables, Extents, Tables -using RecipesBase - -# Types: could we move these to a shapes.jl file? - -""" - Rect - -A rectangle object to represent the bounding box for other shape file shapes. -""" -struct Rect - left::Float64 - bottom::Float64 - right::Float64 - top::Float64 -end - -abstract type AbstractShape end - -isgeometry(::AbstractShape) = true -GeoInterface.ncoord(::AbstractShape) = 2 # With specific methods when 3 - -GeoInterface.@enable_geo_plots AbstractShape - -""" - Interval - -Represents the range of measures or Z dimension, in a shape file. -""" -struct Interval - left::Float64 - right::Float64 -end - -abstract type AbstractPoint <: AbstractShape end - -GeoInterface.geomtype(point::Type{<:AbstractPoint}) = GeoInterface.PointTrait() -GeoInterface.x(point::AbstractPoint) = point.x -GeoInterface.y(point::AbstractPoint) = point.y - -""" - Point <: AbstractPoint - -Point from a shape file. - -Fields `x`, `y` hold the spatial location. -""" -struct Point <: AbstractPoint - x::Float64 - y::Float64 -end - -""" - PointM <: AbstractPoint - -Point from a shape file. - -Fields `x`, `y` hold the spatial location. - -Includes a measure field `m`, holding a value for the point. -""" -struct PointM <: AbstractPoint - x::Float64 - y::Float64 - m::Float64 # measure -end - -GeoInterface.m(point::PointM) = point.m - -""" - PointZ <: AbstractPoint - -Three dimensional point, from a shape file. - -Fields `x`, `y`, `z` hold the spatial location. - -Includes a measure field `m`, holding a value for the point. -""" -struct PointZ <: AbstractPoint - x::Float64 - y::Float64 - z::Float64 - m::Float64 # measure -end - -GeoInterface.m(point::PointZ) = point.m -GeoInterface.z(point::PointZ) = point.z -GeoInterface.ncoord(::PointZ) = 3 - -abstract type AbstractPolyline <: AbstractShape end - -GeoInterface.geomtype(::Type{<:AbstractPolyline}) = GeoInterface.MultiLineStringTrait() - -""" - Polyline <: AbstractPolyline - -Represents a single or multiple polylines from a shape file. - -# Fields -- `points`: a `Vector` of [`Point`](@ref) represents a one or multiple lines. -- `parts`: a `Vector` of `Int32` indicating the line each point belongs to. -- `MBR`: `nothing` or the known bounding box. Can be retrieved with `GeoInterface.bbox`. -""" -struct Polyline <: AbstractPolyline - MBR::Rect - parts::Vector{Int32} - points::Vector{Point} -end - -""" - PolylineM <: AbstractPolyline - -Polyline from a shape file, with measures. - -# Fields -- `points`: a `Vector` of [`Point`](@ref) represents a one or multiple lines. -- `parts`: a `Vector` of `Int32` indicating the line each point belongs to. -- `MBR`: `nothing` or the known bounding box. Can be retrieved with `GeoInterface.bbox`. -- `measures`: holds values from each point. -""" -struct PolylineM <: AbstractPolyline - MBR::Rect - parts::Vector{Int32} - points::Vector{Point} - measures::Vector{Float64} -end - -""" - PolylineZ <: AbstractPolyline - -Three dimensional polyline of from a shape file. - -# Fields -- `points`: a `Vector` of [`Point`](@ref) represents a one or multiple lines. -- `parts`: a `Vector` of `Int32` indicating the line each point belongs to. -- `zvalues`: a `Vector` of `Float64` representing the z dimension values. -- `MBR`: `nothing` or the known bounding box. Can be retrieved with `GeoInterface.bbox`. -- `measures`: holds values from each point. -""" -struct PolylineZ <: AbstractPolyline - MBR::Rect - parts::Vector{Int32} - points::Vector{Point} - zvalues::Vector{Float64} - measures::Vector{Float64} -end - -GeoInterface.ncoord(::PolylineZ) = 3 - -abstract type AbstractPolygon <: AbstractShape end - -GeoInterface.geomtype(::Type{<:AbstractPolygon}) = GeoInterface.PolygonTrait() - -""" - Polygon <: AbstractPolygon - -Represents a polygon from a shape file. - -# Fields -- `points`: a `Vector` of [`Point`](@ref) represents a one or multiple closed areas. -- `parts`: a `Vector` of `Int32` indicating the polygon each point belongs to. -- `MBR`: `nothing` or the known bounding box. Can be retrieved with `GeoInterface.bbox`. -""" -struct Polygon <: AbstractPolygon - MBR::Rect - parts::Vector{Int32} - points::Vector{Point} -end - -Base.show(io::IO, p::Polygon) = - print(io, "Polygon(", length(p.points), " Points)") - -""" - PolygonM <: AbstractPolygon - -Represents a polygon from a shape file - -# Fields -- `points`: a `Vector` of [`Point`](@ref) represents a one or multiple closed areas. -- `parts`: a `Vector` of `Int32` indicating the polygon each point belongs to. -- `MBR`: `nothing` or the known bounding box. Can be retrieved with `GeoInterface.bbox`. -- `measures`: holds values from each point. -""" -struct PolygonM <: AbstractPolygon - MBR::Rect - parts::Vector{Int32} - points::Vector{Point} - measures::Vector{Float64} -end - -""" - PolygonZ <: AbstractPolygon - -A three dimensional polygon from a shape file. - -# Fields -- `points`: a `Vector` of [`Point`](@ref) represents a one or multiple closed areas. -- `parts`: a `Vector` of `Int32` indicating the polygon each point belongs to. -- `zvalues`: a `Vector` of `Float64` representing the z dimension values. -- `MBR`: `nothing` or the known bounding box. Can be retrieved with `GeoInterface.bbox`. -- `measures`: holds values from each point. -""" -struct PolygonZ <: AbstractPolygon - MBR::Rect - parts::Vector{Int32} - points::Vector{Point} - zvalues::Vector{Float64} - measures::Vector{Float64} -end +const GI = GeoInterface -GeoInterface.ncoord(::PolygonZ) = 3 - -abstract type AbstractMultiPoint <: AbstractShape end - -GeoInterface.geomtype(::Type{<:AbstractMultiPoint}) = GeoInterface.MultiPointTrait() - -""" - MultiPoint <: AbstractMultiPoint - -Collection of points, from a shape file. - -# Fields -- `points`: a `Vector` of [`Point`](@ref). -- `MBR`: `nothing` or the known bounding box. Can be retrieved with `GeoInterface.bbox`. -""" -struct MultiPoint <: AbstractMultiPoint - MBR::Rect - points::Vector{Point} -end - -""" - MultiPointM <: AbstractMultiPoint - -Collection of points, from a shape file. - -Includes a `measures` field, holding values from each point. - -May have a known bounding box, which can be retrieved with `GeoInterface.bbox`. - -# Fields -- `points`: a `Vector` of [`Point`](@ref). -- `MBR`: `nothing` or the known bounding box. Can be retrieved with `GeoInterface.bbox`. -- `measures`: holds values from each point. -""" -struct MultiPointM <: AbstractMultiPoint - MBR::Rect - points::Vector{Point} - measures::Vector{Float64} -end - -""" - MultiPointZ <: AbstractMultiPoint - -Collection of 3d points, from a shape file. - -Includes a `measures` field, holding values from each point. - -May have a known bounding box, which can be retrieved with `GeoInterface.bbox`. - -# Fields -- `points`: a `Vector` of [`Point`](@ref). -- `zvalues`: a `Vector` of `Float64` representing the z dimension values. -- `MBR`: `nothing` or the known bounding box. Can be retrieved with `GeoInterface.bbox`. -- `measures`: holds values from each point. -""" -struct MultiPointZ <: AbstractMultiPoint - MBR::Rect - points::Vector{Point} - zvalues::Vector{Float64} - measures::Vector{Float64} -end - -GeoInterface.ncoord(::MultiPointZ) = 3 - -""" - MultiPatch - -Stores a collection of patch representing the boundary of a 3d object. - -# Fields -- `points`: a `Vector` of [`Point`](@ref) represents a one or multiple spatial objects. -- `parts`: a `Vector` of `Int32` indicating the object each point belongs to. -- `parttypes`: a `Vector` of `Int32` indicating the type of object each point belongs to. -- `zvalues`: a `Vector` of `Float64` indicating absolute or relative heights. -- `MBR`: `nothing` or the known bounding box. Can be retrieved with `GeoInterface.bbox`. -""" -struct MultiPatch <: AbstractShape - MBR::Rect - parts::Vector{Int32} - parttypes::Vector{Int32} - points::Vector{Point} - zvalues::Vector{Float64} - # measures::Vector{Float64} # (optional) -end +using RecipesBase -GeoInterface.geomtype(::Type{<:MultiPatch}) = GeoInterface.MultiPointTrait() +include("common.jl") +include("points.jl") +include("polygons.jl") +include("polylines.jl") +include("multipoints.jl") +include("multipatch.jl") -GeoInterface.ncoord(::MultiPatch) = 3 +# Handle/Table const SHAPETYPE = Dict{Int32,DataType}( 0 => Missing, @@ -317,288 +32,10 @@ const SHAPETYPE = Dict{Int32,DataType}( 31 => MultiPatch, ) -""" - Handle - - Handle(path::AbstractString, [indexpath::AbstractString]) - -Load a shapefile into GeoInterface compatible objects. This can be plotted -with Plots.jl `plot`. - -The Vector of shape object can be accessed with `shapes(handle)`. - -`Handle` may have a known bounding box, which can be retrieved with `GeoInterface.bbox`. -""" -mutable struct Handle{T<:Union{<:AbstractShape,Missing}} - code::Int32 - length::Int32 - version::Int32 - shapeType::Int32 - MBR::Rect - zrange::Interval - mrange::Interval - shapes::Vector{T} - crs::Union{Nothing,GeoFormatTypes.ESRIWellKnownText} -end -function Handle(path::AbstractString, index=nothing) - open(path) do io - read(io, Handle, index; path = path) - end -end - -shapes(h::Handle) = h.shapes - -GeoInterface.crs(h::Handle) = h.crs - -Base.length(shp::Handle) = length(shapes(shp)) - -function Base.read(io::IO, ::Type{Rect}) - minx = read(io, Float64) - miny = read(io, Float64) - maxx = read(io, Float64) - maxy = read(io, Float64) - Rect(minx, miny, maxx, maxy) -end - -function Base.read(io::IO, ::Type{Point}) - x = read(io, Float64) - y = read(io, Float64) - Point(x, y) -end - -function Base.read(io::IO, ::Type{PointM}) - x = read(io, Float64) - y = read(io, Float64) - m = read(io, Float64) - PointM(x, y, m) -end - -function Base.read(io::IO, ::Type{PointZ}) - x = read(io, Float64) - y = read(io, Float64) - z = read(io, Float64) - m = read(io, Float64) - PointZ(x, y, z, m) -end - -function Base.read(io::IO, ::Type{Polyline}) - box = read(io, Rect) - numparts = read(io, Int32) - numpoints = read(io, Int32) - parts = Vector{Int32}(undef, numparts) - read!(io, parts) - points = Vector{Point}(undef, numpoints) - read!(io, points) - Polyline(box, parts, points) -end - -function Base.read(io::IO, ::Type{PolylineM}) - box = read(io, Rect) - numparts = read(io, Int32) - numpoints = read(io, Int32) - parts = Vector{Int32}(undef, numparts) - read!(io, parts) - points = Vector{Point}(undef, numpoints) - read!(io, points) - mrange = Vector{Float64}(undef, 2) - read!(io, mrange) - measures = Vector{Float64}(undef, numpoints) - read!(io, measures) - PolylineM(box, parts, points, measures) -end - -function Base.read(io::IO, ::Type{PolylineZ}) - box = read(io, Rect) - numparts = read(io, Int32) - numpoints = read(io, Int32) - parts = Vector{Int32}(undef, numparts) - read!(io, parts) - points = Vector{Point}(undef, numpoints) - read!(io, points) - zrange = Vector{Float64}(undef, 2) - read!(io, zrange) - zvalues = Vector{Float64}(undef, numpoints) - read!(io, zvalues) - mrange = Vector{Float64}(undef, 2) - read!(io, mrange) - measures = Vector{Float64}(undef, numpoints) - read!(io, measures) - PolylineZ(box, parts, points, zvalues, measures) -end - -function Base.read(io::IO, ::Type{Polygon}) - box = read(io, Rect) - numparts = read(io, Int32) - numpoints = read(io, Int32) - parts = Vector{Int32}(undef, numparts) - read!(io, parts) - points = Vector{Point}(undef, numpoints) - read!(io, points) - Polygon(box, parts, points) -end - -function Base.read(io::IO, ::Type{PolygonM}) - box = read(io, Rect) - numparts = read(io, Int32) - numpoints = read(io, Int32) - parts = Vector{Int32}(undef, numparts) - read!(io, parts) - points = Vector{Point}(undef, numpoints) - read!(io, points) - mrange = Vector{Float64}(undef, 2) - read!(io, mrange) - measures = Vector{Float64}(undef, numpoints) - read!(io, measures) - PolygonM(box, parts, points, measures) -end - -function Base.read(io::IO, ::Type{PolygonZ}) - box = read(io, Rect) - numparts = read(io, Int32) - numpoints = read(io, Int32) - parts = Vector{Int32}(undef, numparts) - read!(io, parts) - points = Vector{Point}(undef, numpoints) - read!(io, points) - zrange = Vector{Float64}(undef, 2) - read!(io, zrange) - zvalues = Vector{Float64}(undef, numpoints) - read!(io, zvalues) - mrange = Vector{Float64}(undef, 2) - read!(io, mrange) - measures = Vector{Float64}(undef, numpoints) - read!(io, measures) - PolygonZ(box, parts, points, zvalues, measures) -end - -function Base.read(io::IO, ::Type{MultiPoint}) - box = read(io, Rect) - numpoints = read(io, Int32) - points = Vector{Point}(undef, numpoints) - read!(io, points) - MultiPoint(box, points) -end - -function Base.read(io::IO, ::Type{MultiPointM}) - box = read(io, Rect) - numpoints = read(io, Int32) - points = Vector{Point}(undef, numpoints) - read!(io, points) - mrange = Vector{Float64}(undef, 2) - read!(io, mrange) - measures = Vector{Float64}(undef, numpoints) - read!(io, measures) - MultiPointM(box, points, measures) -end - -function Base.read(io::IO, ::Type{MultiPointZ}) - box = read(io, Rect) - numpoints = read(io, Int32) - points = Vector{Point}(undef, numpoints) - read!(io, points) - zrange = Vector{Float64}(undef, 2) - read!(io, zrange) - zvalues = Vector{Float64}(undef, numpoints) - read!(io, zvalues) - mrange = Vector{Float64}(undef, 2) - read!(io, mrange) - measures = Vector{Float64}(undef, numpoints) - read!(io, measures) - MultiPointZ(box, points, zvalues, measures) -end - -function Base.read(io::IO, ::Type{MultiPatch}) - box = read(io, Rect) - numparts = read(io, Int32) - numpoints = read(io, Int32) - parts = Vector{Int32}(undef, numparts) - read!(io, parts) - parttypes = Vector{Int32}(undef, numparts) - read!(io, parttypes) - points = Vector{Point}(undef, numpoints) - read!(io, points) - zrange = Vector{Float64}(undef, 2) - read!(io, zrange) - zvalues = Vector{Float64}(undef, numpoints) - read!(io, zvalues) - # mrange = Vector{Float64}(2) - # read!(io, mrange) - # measures = Vector{Float64}(numpoints) - # read!(io, measures) - MultiPatch(box, parts, parttypes, points, zvalues) #,measures) -end - -function Base.read(io::IO, ::Type{Handle}, index = nothing; path = nothing) - code = bswap(read(io, Int32)) - read!(io, Vector{Int32}(undef, 5)) - fileSize = bswap(read(io, Int32)) - version = read(io, Int32) - shapeType = read(io, Int32) - MBR = read(io, Rect) - zmin = read(io, Float64) - zmax = read(io, Float64) - mmin = read(io, Float64) - mmax = read(io, Float64) - jltype = SHAPETYPE[shapeType] - shapes = Vector{Union{jltype,Missing}}(undef, 0) - crs = nothing - if path !== nothing - prjfile = string(splitext(path)[1], ".prj") - if isfile(prjfile) - try - crs = GeoFormatTypes.ESRIWellKnownText(read(open(prjfile), String)) - catch - @warn "Projection file $prjfile appears to be corrupted. `nothing` used for `crs`" - end - end - end - file = Handle( - code, - fileSize, - version, - shapeType, - MBR, - Interval(zmin, zmax), - Interval(mmin, mmax), - shapes, - crs, - ) - num = Int32(0) - while (!eof(io)) - seeknext(io, num, index) - num = bswap(read(io, Int32)) - rlength = bswap(read(io, Int32)) - shapeType = read(io, Int32) - if shapeType === Int32(0) - push!(shapes, missing) - else - push!(shapes, read(io, jltype)) - end - end - file -end - include("shx.jl") +include("handle.jl") include("table.jl") -include("geo_interface.jl") +include("extent.jl") include("plotrecipes.jl") -seeknext(io, num, ::Nothing) = nothing - -function seeknext(io, num, index::IndexHandle) - seek(io, index.indices[num+1].offset * 2) -end - -function Handle(path::AbstractString, indexpath::AbstractString) - index = open(indexpath) do io - read(io, IndexHandle) - end - Handle(path, index) -end - -function Base.:(==)(a::Rect, b::Rect) - a.left == b.left && - a.bottom == b.bottom && a.right == b.right && a.top == b.top -end - end # module diff --git a/src/common.jl b/src/common.jl new file mode 100644 index 0000000..e52c75a --- /dev/null +++ b/src/common.jl @@ -0,0 +1,46 @@ + +abstract type AbstractShape end + +isgeometry(::AbstractShape) = true +GeoInterface.ncoord(::GI.AbstractGeometryTrait, ::AbstractShape) = 2 # With specific methods when 3 + +""" + Rect + +A rectangle object to represent the bounding box for other shape file shapes. +""" +struct Rect + left::Float64 + bottom::Float64 + right::Float64 + top::Float64 +end + +function Base.read(io::IO, ::Type{Rect}) + minx = read(io, Float64) + miny = read(io, Float64) + maxx = read(io, Float64) + maxy = read(io, Float64) + Rect(minx, miny, maxx, maxy) +end + +function Base.:(==)(a::Rect, b::Rect) + a.left == b.left && + a.bottom == b.bottom && a.right == b.right && a.top == b.top +end + +""" + Interval + +Represents the range of measures or Z dimension, in a shape file. +""" +struct Interval + left::Float64 + right::Float64 +end + +function Base.read(io::IO, ::Type{Interval}) + left = read(io, Float64) + right = read(io, Float64) + Interval(left, right) +end diff --git a/src/extent.jl b/src/extent.jl new file mode 100644 index 0000000..0afe46c --- /dev/null +++ b/src/extent.jl @@ -0,0 +1,14 @@ +const Shape3D = Union{PolylineZ,PolygonZ,MultiPointZ,MultiPatch} + +# extent +GI.is3d(::GI.AbstractGeometryTrait, ::AbstractShape) = false +GI.is3d(::GI.AbstractGeometryTrait, ::Shape3D) = true + +function GI.extent(x::AbstractShape) + rect = x.MBR + return Extents.Extent(X=(rect.left, rect.right), Y=(rect.bottom, rect.top)) +end +function GI.extent(x::Union{Shape3D,Handle}) + rect = x.MBR + return Extents.Extent(X=(rect.left, rect.right), Y=(rect.bottom, rect.top), Z=(x.zrange.left, x.zrange.right)) +end diff --git a/src/geo_interface.jl b/src/geo_interface.jl deleted file mode 100644 index 4463b37..0000000 --- a/src/geo_interface.jl +++ /dev/null @@ -1,83 +0,0 @@ -GeoInterface.coordinates(::GeoInterface.PointTrait, obj::Union{Point,PointM}) = [obj.x, obj.y] -GeoInterface.coordinates(::GeoInterface.PointTrait, obj::PointZ) = [obj.x, obj.y, obj.z] - -GeoInterface.coordinates(obj::Union{MultiPoint,MultiPointM,MultiPointZ}) = - Vector{Float64}[GeoInterface.coordinates(p) for p in obj.points] - -function GeoInterface.coordinates(::GeoInterface.MultiLineStringTrait, obj::Union{Polyline,PolylineM,PolylineZ}) - npoints = length(obj.points) - nparts = length(obj.parts) - @assert obj.parts[nparts] <= npoints - push!(obj.parts, npoints) - coords = Vector{Vector{Float64}}[Vector{Float64}[GeoInterface.coordinates(obj.points[j+1]) for j = obj.parts[i]:(obj.parts[i+1]-1)] for i = 1:nparts] - pop!(obj.parts) - return coords -end - -# Only supports 2D geometries for now -# pt is [x,y] and ring is [[x,y], [x,y],..] -# ported from https://github.com/Turfjs/turf/blob/3d0efd96d59878f82c6d6baf8ed75695e0b6dbc0/packages/turf-inside/index.js#L76-L91 -function inring(pt::Vector{Float64}, ring::Vector{Vector{Float64}}) - intersect(i::Vector{Float64}, j::Vector{Float64}) = - (i[2] >= pt[2]) != (j[2] >= pt[2]) && (pt[1] <= (j[1] - i[1]) * - (pt[2] - i[2]) / - (j[2] - i[2]) + i[1]) - isinside = intersect(ring[1], ring[end]) - for k = 2:length(ring) - isinside = intersect(ring[k], ring[k-1]) ? !isinside : isinside - end - return isinside -end - -# ported from https://github.com/Esri/terraformer-arcgis-parser/blob/master/terraformer-arcgis-parser.js#L168-L253 -function GeoInterface.coordinates(::GeoInterface.PolygonTrait, obj::Union{Polygon,PolygonZ,PolygonM}) - npoints = length(obj.points) - nparts = length(obj.parts) - @assert obj.parts[end] <= npoints - coords = Vector{Vector{Vector{Float64}}}[] - holes = Vector{Vector{Vector{Float64}}}() - push!(obj.parts, npoints) - for i = 1:nparts - ringlength = obj.parts[i+1] - obj.parts[i] - ringlength < 4 && continue - test = 0.0 - ring = Vector{Vector{Float64}}() - for j = (obj.parts[i]+1):(obj.parts[i+1]-1) - prev = obj.points[j] - cur = obj.points[j+1] - test += (cur.x - prev.x) * (cur.y + prev.y) - push!(ring, GeoInterface.coordinates(prev)) - end - push!(ring, GeoInterface.coordinates(obj.points[obj.parts[i+1]])) - @assert length(ring) == ringlength - if test > 0 # clockwise - push!(coords, Vector{Vector{Float64}}[ring]) # create new polygon - else # anti-clockwise - push!(holes, ring) - end - end - pop!(obj.parts) - nrings = length(coords) - for hole in holes - found = false - for i = 1:nrings - if inring(hole[1], coords[i][1]) - push!(coords[i], hole) - found = true - break - end - end - if !found - push!(coords, [hole]) # hole is not inside any ring; make it a polygon - end - end - return coords -end - -# bbox -const HasMBR = Union{AbstractPolyline,AbstractPolygon,AbstractMultiPoint,MultiPatch,Handle} - -function GeoInterface.extent(x::HasMBR) - rect = x.MBR - return Extents.Extent(X=(rect.left, rect.right), Y=(rect.bottom, rect.top)) -end diff --git a/src/handle.jl b/src/handle.jl new file mode 100644 index 0000000..e23ae0c --- /dev/null +++ b/src/handle.jl @@ -0,0 +1,94 @@ +""" + Handle + + Handle(path::AbstractString, [indexpath::AbstractString]) + +Load a shapefile into GeoInterface compatible objects. This can be plotted +with Plots.jl `plot`. + +The Vector of shape object can be accessed with `shapes(handle)`. + +`Handle` may have a known bounding box, which can be retrieved with `GeoInterface.bbox`. +""" +mutable struct Handle{T<:Union{<:AbstractShape,Missing}} + code::Int32 + length::Int32 + version::Int32 + shapeType::Int32 + MBR::Rect + zrange::Interval + mrange::Interval + shapes::Vector{T} + crs::Union{Nothing,GeoFormatTypes.ESRIWellKnownText} +end +function Handle(path::AbstractString, index=nothing) + open(path) do io + read(io, Handle, index; path = path) + end +end +function Handle(path::AbstractString, indexpath::AbstractString) + index = open(indexpath) do io + read(io, IndexHandle) + end + Handle(path, index) +end + +shapes(h::Handle) = h.shapes + +GeoInterface.crs(h::Handle) = h.crs + +Base.length(shp::Handle) = length(shapes(shp)) + +function Base.read(io::IO, ::Type{Handle}, index = nothing; path = nothing) + code = bswap(read(io, Int32)) + read!(io, Vector{Int32}(undef, 5)) + fileSize = bswap(read(io, Int32)) + version = read(io, Int32) + shapeType = read(io, Int32) + MBR = read(io, Rect) + zmin = read(io, Float64) + zmax = read(io, Float64) + mmin = read(io, Float64) + mmax = read(io, Float64) + jltype = SHAPETYPE[shapeType] + shapes = Vector{Union{jltype,Missing}}(undef, 0) + crs = nothing + if path !== nothing + prjfile = string(splitext(path)[1], ".prj") + if isfile(prjfile) + try + crs = GeoFormatTypes.ESRIWellKnownText(read(open(prjfile), String)) + catch + @warn "Projection file $prjfile appears to be corrupted. `nothing` used for `crs`" + end + end + end + num = Int32(0) + while (!eof(io)) + seeknext(io, num, index) + num = bswap(read(io, Int32)) + rlength = bswap(read(io, Int32)) + shapeType = read(io, Int32) + if shapeType === Int32(0) + push!(shapes, missing) + else + push!(shapes, read(io, jltype)) + end + end + return Handle( + code, + fileSize, + version, + shapeType, + MBR, + Interval(zmin, zmax), + Interval(mmin, mmax), + shapes, + crs, + ) +end + +function seeknext(io, num, index::IndexHandle) + seek(io, index.indices[num+1].offset * 2) +end +seeknext(io, num, ::Nothing) = nothing diff --git a/src/multipatch.jl b/src/multipatch.jl new file mode 100644 index 0000000..8c7e7bb --- /dev/null +++ b/src/multipatch.jl @@ -0,0 +1,44 @@ +""" + MultiPatch + +Stores a collection of patch representing the boundary of a 3d object. + +# Fields +- `points`: a `Vector` of [`Point`](@ref) represents a one or multiple spatial objects. +- `parts`: a `Vector` of `Int32` indicating the object each point belongs to. +- `parttypes`: a `Vector` of `Int32` indicating the type of object each point belongs to. +- `zvalues`: a `Vector` of `Float64` indicating absolute or relative heights. +- `MBR`: `nothing` or the known bounding box. Can be retrieved with `GeoInterface.bbox`. +""" +struct MultiPatch <: AbstractShape + MBR::Rect + parts::Vector{Int32} + parttypes::Vector{Int32} + points::Vector{Point} + zrange::Interval + zvalues::Vector{Float64} + # measures::Vector{Float64} # (optional) +end + +GeoInterface.geomtype(::MultiPatch) = GeoInterface.MultiPolygonTrait() +GeoInterface.ncoord(::MultiPatch) = 3 + +function Base.read(io::IO, ::Type{MultiPatch}) + box = read(io, Rect) + numparts = read(io, Int32) + numpoints = read(io, Int32) + parts = Vector{Int32}(undef, numparts) + read!(io, parts) + parttypes = Vector{Int32}(undef, numparts) + read!(io, parttypes) + points = Vector{Point}(undef, numpoints) + read!(io, points) + zrange = read(io, Interval) + zvalues = Vector{Float64}(undef, numpoints) + read!(io, zvalues) + # mrange = Vector{Float64}(2) + # read!(io, mrange) + # measures = Vector{Float64}(numpoints) + # read!(io, measures) + MultiPatch(box, parts, parttypes, points, zrange, zvalues) #,measures) +end diff --git a/src/multipoints.jl b/src/multipoints.jl new file mode 100644 index 0000000..56d3dbd --- /dev/null +++ b/src/multipoints.jl @@ -0,0 +1,107 @@ + +abstract type AbstractMultiPoint{T} <: AbstractShape end + +GI.geomtype(::AbstractMultiPoint) = GI.MultiPointTrait() +GI.ngeom(::GI.MultiPointTrait, geom::AbstractMultiPoint) = length(geom.points) +GI.ncoord(::GI.MultiPointTrait, geom::AbstractMultiPoint{T}) where T = _ncoord(T) +GI.npoint(::GI.MultiPointTrait, geom::AbstractMultiPoint) = length(geom.points) + +""" + MultiPoint <: AbstractMultiPoint + +Collection of points, from a shape file. + +# Fields +- `points`: a `Vector` of [`Point`](@ref). +- `MBR`: `nothing` or the known bounding box. Can be retrieved with `GeoInterface.bbox`. +""" +struct MultiPoint <: AbstractMultiPoint{Point} + MBR::Rect + points::Vector{Point} +end + +function Base.read(io::IO, ::Type{MultiPoint}) + box = read(io, Rect) + numpoints = read(io, Int32) + points = Vector{Point}(undef, numpoints) + read!(io, points) + MultiPoint(box, points) +end + +GI.getgeom(::GI.MultiPointTrait, geom::MultiPoint, i::Integer) = geom.points[i] + +""" + MultiPointM <: AbstractMultiPoint + +Collection of points, from a shape file. + +Includes a `measures` field, holding values from each point. + +May have a known bounding box, which can be retrieved with `GeoInterface.bbox`. + +# Fields +- `points`: a `Vector` of [`Point`](@ref). +- `MBR`: `nothing` or the known bounding box. Can be retrieved with `GeoInterface.bbox`. +- `measures`: holds values from each point. +""" +struct MultiPointM <: AbstractMultiPoint{PointM} + MBR::Rect + points::Vector{Point} + mrange::Interval + measures::Vector{Float64} +end + +function Base.read(io::IO, ::Type{MultiPointM}) + box = read(io, Rect) + numpoints = read(io, Int32) + points = Vector{Point}(undef, numpoints) + read!(io, points) + mrange = read(io, Interval) + measures = Vector{Float64}(undef, numpoints) + read!(io, measures) + MultiPointM(box, points, mrange, measures) +end + +GI.getgeom(::GI.MultiPointTrait, geom::MultiPointM, i::Integer) = + PointM(geom.points[i], geom.measures[i]) + +""" + MultiPointZ <: AbstractMultiPoint + +Collection of 3d points, from a shape file. + +Includes a `measures` field, holding values from each point. + +May have a known bounding box, which can be retrieved with `GeoInterface.bbox`. + +# Fields +- `points`: a `Vector` of [`Point`](@ref). +- `zvalues`: a `Vector` of `Float64` representing the z dimension values. +- `MBR`: `nothing` or the known bounding box. Can be retrieved with `GeoInterface.bbox`. +- `measures`: holds values from each point. +""" +struct MultiPointZ <: AbstractMultiPoint{PointZ} + MBR::Rect + points::Vector{Point} + zrange::Interval + zvalues::Vector{Float64} + mrange::Interval + measures::Vector{Float64} +end + +GI.getgeom(::GI.MultiPointTrait, geom::MultiPointZ, i::Integer) = + PointZ(geom.points[i], geom.zvalues[i], geom.measures[i]) + +function Base.read(io::IO, ::Type{MultiPointZ}) + box = read(io, Rect) + numpoints = read(io, Int32) + points = Vector{Point}(undef, numpoints) + read!(io, points) + zrange = read(io, Interval) + zvalues = Vector{Float64}(undef, numpoints) + read!(io, zvalues) + mrange = read(io, Interval) + measures = Vector{Float64}(undef, numpoints) + read!(io, measures) + MultiPointZ(box, points, zrange, zvalues, mrange, measures) +end diff --git a/src/plotrecipes.jl b/src/plotrecipes.jl index 48b4aa2..7f04e77 100644 --- a/src/plotrecipes.jl +++ b/src/plotrecipes.jl @@ -5,3 +5,5 @@ end @recipe function f(shp::Handle) shapes(shp) end + +GeoInterfaceRecipes.@enable_geo_plots Shapefile.AbstractShape diff --git a/src/points.jl b/src/points.jl new file mode 100644 index 0000000..f050368 --- /dev/null +++ b/src/points.jl @@ -0,0 +1,84 @@ + +abstract type AbstractPoint <: AbstractShape end + +GI.geomtype(::AbstractPoint) = GI.PointTrait() +GI.x(::GI.PointTrait, point::AbstractPoint) = point.x +GI.y(::GI.PointTrait, point::AbstractPoint) = point.y +GI.getcoord(::GI.PointTrait, p::AbstractPoint, i::Integer) = getfield(p, i) +GI.ncoord(::GI.PointTrait, p::T) where T <: AbstractPoint = _ncoord(T) + +_ncoord(::Type{<:AbstractPoint}) = 2 + +""" + Point <: AbstractPoint + +Point from a shape file. + +Fields `x`, `y` hold the spatial location. +""" +struct Point <: AbstractPoint + x::Float64 + y::Float64 +end + +function Base.read(io::IO, ::Type{Point}) + x = read(io, Float64) + y = read(io, Float64) + Point(x, y) +end + +""" + PointM <: AbstractPoint + +Point from a shape file. + +Fields `x`, `y` hold the spatial location. + +Includes a measure field `m`, holding a value for the point. +""" +struct PointM <: AbstractPoint + x::Float64 + y::Float64 + m::Float64 # measure +end +PointM(p::Point, m::Float64) = PointM(p.x, p.y, m) + +function Base.read(io::IO, ::Type{PointM}) + x = read(io, Float64) + y = read(io, Float64) + m = read(io, Float64) + PointM(x, y, m) +end + +GI.m(::GI.PointTrait, point::PointM) = point.m + +""" + PointZ <: AbstractPoint + +Three dimensional point, from a shape file. + +Fields `x`, `y`, `z` hold the spatial location. + +Includes a measure field `m`, holding a value for the point. +""" +struct PointZ <: AbstractPoint + x::Float64 + y::Float64 + z::Float64 + m::Float64 # measure +end +PointZ(p::Point, z, m) = PointZ(p.x, p.y, z, m) + +function Base.read(io::IO, ::Type{PointZ}) + x = read(io, Float64) + y = read(io, Float64) + z = read(io, Float64) + m = read(io, Float64) + PointZ(x, y, z, m) +end + +GI.m(::GI.PointTrait, point::PointZ) = point.m +GI.z(::GI.PointTrait, point::PointZ) = point.z + +_ncoord(::Type{PointZ}) = 3 + diff --git a/src/polygons.jl b/src/polygons.jl new file mode 100644 index 0000000..57b567b --- /dev/null +++ b/src/polygons.jl @@ -0,0 +1,233 @@ + +# Shapefiles Polygons are actually MultiPolygons Made up of unsorted rings. +# To handle this we add additional types not in the Shapefile specification. +# LinearRing represents individual rings while `SubPolygon` is a single external +# ring and n internal rings (holes). +# Generating polygons is expensive - it must be calculated if they are +# exterior or interior, and which interiors are inside which exteriors + +struct LinearRing{P,Z,M} <: AbstractVector{P} + xy::SubArray{Point,1,Vector{Point},Tuple{UnitRange{Int64}},true} + z::Z + m::M +end +LinearRing{P}(xy, z::Z, m::M) where {P,Z,M} = LinearRing{P,Z,M}(xy, z, m) +LinearRing{P}(xy; z=nothing, m=nothing) where P = LinearRing{P}(xy, z, m) + +GI.isgeometry(::Type{<:LinearRing}) = true +GI.geomtype(::LinearRing) = GI.LinearRingTrait() +GI.ncoord(lr::LinearRing{P}) where P = _ncoord(P) +GI.ngeom(::GI.LinearRingTrait, lr::LinearRing) = length(lr) +# GI.getgeom(::GI.LinearRingTrait, lr::LinearRing) = (getindex(lr, i) for i in 1:ngeom(lr)) + +GI.getgeom(::GI.LinearRingTrait, lr::LinearRing{Point}, i::Integer) = lr.xy[i] +GI.getgeom(::GI.LinearRingTrait, lr::LinearRing{PointM}, i::Integer) = + PointM(lr.xy[i], lr.m[i]) +GI.getgeom(::GI.LinearRingTrait, lr::LinearRing{PointZ}, i::Integer) = + PointZ(lr.xy[i], lr.z[i], lr.m[i]) + +Base.parent(lr::LinearRing) = lr.xy +Base.getindex(lr::LinearRing, i) = getindex(parent(lr), i) +Base.size(p::LinearRing) = (length(p),) +Base.length(lr::LinearRing) = length(parent(lr)) + +struct SubPolygon{L<:LinearRing} <: AbstractVector{L} + rings::Vector{L} +end +GI.isgeometry(::Type{<:SubPolygon}) = true +GI.geomtype(::SubPolygon) = GI.PolygonTrait() +GI.ncoord(::GI.PolygonTrait, ::SubPolygon{LinearRing{P}}) where P = _ncoord(P) +GI.ngeom(::GI.PolygonTrait, sp::SubPolygon) = length(sp) +GI.getgeom(::GI.PolygonTrait, sp::SubPolygon, i::Integer) = getindex(sp, i) + +Base.parent(p::SubPolygon) = p.rings +Base.getindex(p::SubPolygon, i) = getindex(parent(p), i) +Base.length(p::SubPolygon) = length(parent(p)) +Base.size(p::SubPolygon) = (length(p),) +Base.push!(p::SubPolygon, x) = Base.push!(parent(p), x) + + +abstract type AbstractPolygon{T} <: AbstractShape end + +# Shapefile polygons are OGC multipolygons +GI.geomtype(geom::AbstractPolygon) = GI.MultiPolygonTrait() +GI.nring(::GI.MultiPolygonTrait, geom::AbstractPolygon) = length(geom.parts) +function GI.ngeom(::GI.MultiPolygonTrait, geom::AbstractPolygon) + n = 0 + for ring in GI.getring(geom) + if _isclockwise(ring) + n += 1 + end + end + return n +end + +function _isclockwise(ring) + clockwise_test = 0.0 + for i in 1:GI.npoint(ring)-1 + prev = ring[i] + cur = ring[i + 1] + clockwise_test += (cur.x - prev.x) * (cur.y + prev.y) + end + clockwise_test > 0 +end + +function GI.getring(t::GI.MultiPolygonTrait, geom::AbstractPolygon) + return (GI.getring(t, geom, i) for i in eachindex(geom.parts)) +end +function GI.getring(::GI.MultiPolygonTrait, geom::AbstractPolygon{P}, i::Integer) where P + pa = geom.parts + range = pa[i]+1:(i == lastindex(pa) ? lastindex(geom.points) : pa[i+1]) + xy = view(geom.points, range) + m = P <: Union{PointM,PointZ} ? view(geom.measures, range) : nothing + z = P <: Union{PointZ} ? view(geom.zvalues, range) : nothing + LinearRing{P}(xy, z, m) +end + +# function GI.gethole(::GI.MultiPolygonTrait, geom::AbstractPolygon) +# end + +# Warning: getgeom is very slow for a Shapefile. Use `getpolygons`, +# or if you don't need exteriors and holes to be separated, `getring`. +GI.getgeom(::GI.MultiPolygonTrait, geom::AbstractPolygon, i::Integer) = collect(GI.getgeom(geom))[i] +function GI.getgeom(::GI.MultiPolygonTrait, geom::AbstractPolygon{T}) where T + r1 = GI.getring(geom, 1) + polygons = SubPolygon{typeof(r1)}[] + holes = typeof(r1)[] + for ring in GI.getring(geom) + if _isclockwise(ring) + push!(polygons, SubPolygon([ring])) # create new polygon + else + push!(holes, ring) + end + end + for hole in holes + found = false + for i = 1:length(polygons) + if _inring(hole[1], polygons[i][1]) + push!(polygons[i], hole) + found = true + break + end + end + if !found + # TODO: does this follow the spec? this should not happen with a correct file. + push!(polygons, SubPolygon([hole])) # hole is not inside any ring; make it a polygon + end + end + return polygons +end + +function _inring(pt::Point, ring::LinearRing) + intersect(i, j) = + (i.y >= pt.y) != (j.y >= pt.y) && + (pt.x <= (j.x - i.x) * (pt.y - i.y) / (j.y - i.y) + i.x) + isinside = intersect(ring[1], ring[end]) + for k = 2:length(ring) + isinside = intersect(ring[k], ring[k-1]) ? !isinside : isinside + end + return isinside +end + + +""" + Polygon <: AbstractPolygon + +Represents a Polygon from a shape file. + +# Fields +- `points`: a `Vector` of [`Point`](@ref) represents a one or multiple closed areas. +- `parts`: a `Vector` of `Int32` indicating the polygon each point belongs to. +- `MBR`: `nothing` or the known bounding box. Can be retrieved with `GI.bbox`. +""" +struct Polygon <: AbstractPolygon{Point} + MBR::Rect + parts::Vector{Int32} + points::Vector{Point} +end + +Base.show(io::IO, p::Polygon) = print(io, "Polygon(", length(p.points), " Points)") + +function Base.read(io::IO, ::Type{Polygon}) + box = read(io, Rect) + numparts = read(io, Int32) + numpoints = read(io, Int32) + parts = Vector{Int32}(undef, numparts) + read!(io, parts) + points = Vector{Point}(undef, numpoints) + read!(io, points) + Polygon(box, parts, points) +end + +""" + PolygonM <: AbstractPolygon + +Represents a polygon from a shape file + +# Fields +- `points`: a `Vector` of [`Point`](@ref) represents a one or multiple closed areas. +- `parts`: a `Vector` of `Int32` indicating the polygon each point belongs to. +- `MBR`: `nothing` or the known bounding box. Can be retrieved with `GI.bbox`. +- `measures`: holds values from each point. +""" +struct PolygonM <: AbstractPolygon{PointM} + MBR::Rect + parts::Vector{Int32} + points::Vector{Point} + mrange::Interval + measures::Vector{Float64} +end + +function Base.read(io::IO, ::Type{PolygonM}) + box = read(io, Rect) + numparts = read(io, Int32) + numpoints = read(io, Int32) + parts = Vector{Int32}(undef, numparts) + read!(io, parts) + points = Vector{Point}(undef, numpoints) + read!(io, points) + mrange = read(io, Interval) + measures = Vector{Float64}(undef, numpoints) + read!(io, measures) + PolygonM(box, parts, points, mrange, measures) +end + +""" + PolygonZ <: AbstractPolygon + +A three dimensional polygon from a shape file. + +# Fields +- `points`: a `Vector` of [`Point`](@ref) represents a one or multiple closed areas. +- `parts`: a `Vector` of `Int32` indicating the polygon each point belongs to. +- `zvalues`: a `Vector` of `Float64` representing the z dimension values. +- `MBR`: `nothing` or the known bounding box. Can be retrieved with `GI.bbox`. +- `measures`: holds values from each point. +""" +struct PolygonZ <: AbstractPolygon{PointZ} + MBR::Rect + parts::Vector{Int32} + points::Vector{Point} + zrange::Interval + zvalues::Vector{Float64} + mrange::Interval + measures::Vector{Float64} +end + +function Base.read(io::IO, ::Type{PolygonZ}) + box = read(io, Rect) + numparts = read(io, Int32) + numpoints = read(io, Int32) + parts = Vector{Int32}(undef, numparts) + read!(io, parts) + points = Vector{Point}(undef, numpoints) + read!(io, points) + zrange = read(io, Interval) + zvalues = Vector{Float64}(undef, numpoints) + read!(io, zvalues) + mrange = read(io, Interval) + measures = Vector{Float64}(undef, numpoints) + read!(io, measures) + PolygonZ(box, parts, points, zrange, zvalues, mrange, measures) +end + diff --git a/src/polylines.jl b/src/polylines.jl new file mode 100644 index 0000000..49fa9ee --- /dev/null +++ b/src/polylines.jl @@ -0,0 +1,147 @@ + +struct LineString{P,Z,M} <: AbstractVector{P} + xy::SubArray{Point,1,Vector{Point},Tuple{UnitRange{Int64}},true} + z::Z + m::M +end +LineString{P}(xy, z::Z, m::M) where {P,Z,M} = LineString{P,Z,M}(xy, z, m) +LineString{P}(xy; z=nothing, m=nothing) where P = LineString{P}(xy, z, m) + +GI.isgeometry(::Type{<:LineString}) = true +GI.geomtype(::LineString) = GI.LineStringTrait() +GI.ncoord(lr::LineString{T}) where T = _ncoord(T) +GI.ngeom(::GI.LineStringTrait, lr::LineString) = length(lr) +# GI.getgeom(::GI.LineStringTrait, lr::LineString) = (getindex(lr, i) for i in 1:ngeom(lr)) +GI.getgeom(::GI.LineStringTrait, lr::LineString, i::Integer) = getindex(lr, i) + +Base.parent(lr::LineString) = lr.xy +Base.size(p::LineString) = (length(p),) +Base.length(lr::LineString) = length(parent(lr)) +Base.getindex(lr::LineString{Point}, i) = lr.xy[i] +Base.getindex(lr::LineString{PointM}, i) = PointM(lr.xy[i], lr.m[i]) +Base.getindex(lr::LineString{PointZ}, i) = PointZ(lr.xy[i], lr.m[i], lr.z[i]) + + +abstract type AbstractPolyline{T} <: AbstractShape end + +GI.geomtype(::AbstractPolyline) = GI.MultiLineStringTrait() +GI.ngeom(::GI.MultiLineStringTrait, geom::AbstractPolyline) = length(geom.parts) +GI.ncoord(::GI.MultiLineStringTrait, ::AbstractPolyline{T}) where T = _ncoord(T) +function GI.getgeom(::GI.MultiLineStringTrait, geom::AbstractPolyline{P}, i::Integer) where P + range = if i < GI.ngeom(geom) + # The linestring is the points between two parts + # C integers start at zero so we add 1 to the range start + (geom.parts[i] + 1):(geom.parts[i+1]) + else + # For the last inestring use the vector lastindex as the last point + (geom.parts[i] + 1):lastindex(geom.points) + end + return LineString(geom, range) +end + +""" + Polyline <: AbstractPolyline + +Represents a single or multiple polylines from a shape file. + +# Fields +- `points`: a `Vector` of [`Point`]() represents a one or multiple lines. +- `parts`: a `Vector` of `Int32` indicating the line each point belongs to. +- `MBR`: `nothing` or the known bounding box. Can be retrieved with `GeoInterface.bbox`. +""" +struct Polyline <: AbstractPolyline{Point} + MBR::Rect + parts::Vector{Int32} + points::Vector{Point} +end + +function Base.read(io::IO, ::Type{Polyline}) + box = read(io, Rect) + numparts = read(io, Int32) + numpoints = read(io, Int32) + parts = Vector{Int32}(undef, numparts) + read!(io, parts) + points = Vector{Point}(undef, numpoints) + read!(io, points) + Polyline(box, parts, points) +end + +LineString(geom::Polyline, range) = @views LineString{Point}(geom.points[range]) + +""" + PolylineM <: AbstractPolyline + +Polyline from a shape file, with measures. + +# Fields +- `points`: a `Vector` of [`Point`](@ref) represents a one or multiple lines. +- `parts`: a `Vector` of `Int32` indicating the line each point belongs to. +- `MBR`: `nothing` or the known bounding box. Can be retrieved with `GeoInterface.bbox`. +- `measures`: holds values from each point. +""" +struct PolylineM <: AbstractPolyline{PointM} + MBR::Rect + parts::Vector{Int32} + points::Vector{Point} + mrange::Interval + measures::Vector{Float64} +end + +LineString(geom::PolylineM, range) = + @views LineString{PointM}(geom.points[range]; m=geom.measures[range]) + +function Base.read(io::IO, ::Type{PolylineM}) + box = read(io, Rect) + numparts = read(io, Int32) + numpoints = read(io, Int32) + parts = Vector{Int32}(undef, numparts) + read!(io, parts) + points = Vector{Point}(undef, numpoints) + read!(io, points) + mrange = read(io, Interval) + measures = Vector{Float64}(undef, numpoints) + read!(io, measures) + PolylineM(box, parts, points, mrange, measures) +end + +""" + PolylineZ <: AbstractPolyline + +Three dimensional polyline of from a shape file. + +# Fields +- `points`: a `Vector` of [`Point`](@ref) represents a one or multiple lines. +- `parts`: a `Vector` of `Int32` indicating the line each point belongs to. +- `zvalues`: a `Vector` of `Float64` representing the z dimension values. +- `MBR`: `nothing` or the known bounding box. Can be retrieved with `GeoInterface.bbox`. +- `measures`: holds values from each point. +""" +struct PolylineZ <: AbstractPolyline{PointZ} + MBR::Rect + parts::Vector{Int32} + points::Vector{Point} + zrange::Interval + zvalues::Vector{Float64} + mrange::Interval + measures::Vector{Float64} +end + +LineString(geom::PolylineZ, range) = + @views LineString{PointZ}(geom.points[range]; z=geom.zvalues[range], m=geom.measures[range]) + +function Base.read(io::IO, ::Type{PolylineZ}) + box = read(io, Rect) + numparts = read(io, Int32) + numpoints = read(io, Int32) + parts = Vector{Int32}(undef, numparts) + read!(io, parts) + points = Vector{Point}(undef, numpoints) + read!(io, points) + zrange = read(io, Interval) + zvalues = Vector{Float64}(undef, numpoints) + read!(io, zvalues) + mrange = read(io, Interval) + measures = Vector{Float64}(undef, numpoints) + read!(io, measures) + PolylineZ(box, parts, points, zrange, zvalues, mrange, measures) +end diff --git a/src/table.jl b/src/table.jl index e2d089f..cfcad17 100644 --- a/src/table.jl +++ b/src/table.jl @@ -145,4 +145,5 @@ Get a vector of the geometries in a shapefile `Table`, without any metadata. """ shapes(t::Table) = shapes(getshp(t)) +GeoInterface.extent(t::Table) = GeoInterface.extent(getshp(t)) GeoInterface.crs(t::Table) = GeoInterface.crs(getshp(t)) diff --git a/test/runtests.jl b/test/runtests.jl index 9dcb85a..a0c7c5a 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -1,89 +1,99 @@ using Shapefile, GeoInterface, Plots, Extents +using Shapefile using Test -plot(Shapefile.Point(1, 2)) +shp = Shapefile.Handle(joinpath("shapelib_testcases", "test.shp")) test_tuples = [ ( path=joinpath("shapelib_testcases", "test.shp"), geomtype=Shapefile.Polygon, - coordinates=Array{Array{Array{Array{Float64,1},1},1},1}[Array{Array{Array{Float64,1},1},1}[Array{Array{Float64,1},1}[Array{Float64,1}[[20.0,20.0],[20.0,30.0],[30.0,30.0],[20.0,20.0]]],Array{Array{Float64,1},1}[Array{Float64,1}[[0.0,0.0],[100.0,0.0],[100.0,100.0],[0.0,100.0],[0.0,0.0]]]],Array{Array{Array{Float64,1},1},1}[Array{Array{Float64,1},1}[Array{Float64,1}[[150.0,150.0],[160.0,150.0],[180.0,170.0],[150.0,150.0]]]],Array{Array{Array{Float64,1},1},1}[Array{Array{Float64,1},1}[Array{Float64,1}[[150.0,150.0],[160.0,150.0],[180.0,170.0],[150.0,150.0]]]]], - extent=Extent(X=(0.0, 180.0), Y=(0.0, 170.0)), + coordinates=[[[[[20.0,20.0],[20.0,30.0],[30.0,30.0],[20.0,20.0]]],[[[0.0,0.0],[100.0,0.0],[100.0,100.0],[0.0,100.0],[0.0,0.0]]]],[[[[150.0,150.0],[160.0,150.0],[180.0,170.0],[150.0,150.0]]]],[[[[150.0,150.0],[160.0,150.0],[180.0,170.0],[150.0,150.0]]]]], + extent=Extent(X=(0.0, 180.0), Y=(0.0, 170.0), Z=(0.0, 0.0)), ),( path=joinpath("shapelib_testcases", "test0.shp"), geomtype=Missing, coordinates=nothing, - extent=Extent(X=(0.0, 10.0), Y=(0.0, 20.0)), + extent=Extent(X=(0.0, 10.0), Y=(0.0, 20.0), Z=(0.0, 0.0)), ),( path=joinpath("shapelib_testcases", "test1.shp"), geomtype=Shapefile.Point, - coordinates=Array{Float64,1}[[1.0,2.0],[10.0,20.0]], - extent=Extent(X=(1.0, 10.0), Y=(2.0, 20.0)), + coordinates=[[1.0,2.0],[10.0,20.0]], + extent=Extent(X=(1.0, 10.0), Y=(2.0, 20.0), Z=(0.0, 0.0)), ),( path=joinpath("shapelib_testcases", "test2.shp"), geomtype=Shapefile.PointZ, - coordinates=Array{Float64,1}[[1.0,2.0,3.0],[10.0,20.0,30.0]], - extent=Extent(X=(1.0, 10.0), Y=(2.0, 20.0)), + coordinates=[[1.0,2.0,3.0],[10.0,20.0,30.0]], + extent=Extent(X=(1.0, 10.0), Y=(2.0, 20.0), Z=(3.0, 30.0)), ),( path=joinpath("shapelib_testcases", "test3.shp"), geomtype=Shapefile.PointM, - coordinates=Array{Float64,1}[[1.0,2.0],[10.0,20.0]], - extent=Extent(X=(1.0, 10.0), Y=(2.0, 20.0)), + coordinates=[[1.0,2.0],[10.0,20.0]], + extent=Extent(X=(1.0, 10.0), Y=(2.0, 20.0), Z=(0.0, 0.0)), ),( path=joinpath("shapelib_testcases", "test4.shp"), geomtype=Shapefile.MultiPoint, - coordinates=Array{Array{Float64,1},1}[Array{Float64,1}[[1.15,2.25],[2.15,3.25],[3.15,4.25],[4.15,5.25]],Array{Float64,1}[[11.15,12.25],[12.15,13.25],[13.15,14.25],[14.15,15.25]],Array{Float64,1}[[21.15,22.25],[22.15,23.25],[23.15,24.25],[24.15,25.25]]], - extent=Extent(X=(1.15, 24.15), Y=(2.25, 25.25)), + coordinates=[[[1.15,2.25],[2.15,3.25],[3.15,4.25],[4.15,5.25]],[[11.15,12.25],[12.15,13.25],[13.15,14.25],[14.15,15.25]],[[21.15,22.25],[22.15,23.25],[23.15,24.25],[24.15,25.25]]], + extent=Extent(X=(1.15, 24.15), Y=(2.25, 25.25), Z=(0.0, 0.0)), ),( path=joinpath("shapelib_testcases", "test5.shp"), geomtype=Shapefile.MultiPointZ, - coordinates=Array{Array{Float64,1},1}[Array{Float64,1}[[1.15,2.25],[2.15,3.25],[3.15,4.25],[4.15,5.25]],Array{Float64,1}[[11.15,12.25],[12.15,13.25],[13.15,14.25],[14.15,15.25]],Array{Float64,1}[[21.15,22.25],[22.15,23.25],[23.15,24.25],[24.15,25.25]]], - extent=Extent(X=(1.15, 24.15), Y=(2.25, 25.25)), + coordinates=[[[1.15,2.25,3.35],[2.15,3.25,4.35],[3.15,4.25,5.35],[4.15,5.25,6.35]],[[11.15,12.25,13.35],[12.15,13.25,14.35],[13.15,14.25,15.35],[14.15,15.25,16.35]],[[21.15,22.25,23.35],[22.15,23.25,24.35],[23.15,24.25,25.35],[24.15,25.25,26.35]]], + extent=Extent(X=(1.15, 24.15), Y=(2.25, 25.25), Z=(3.35, 26.35)), ),( path=joinpath("shapelib_testcases", "test6.shp"), geomtype=Shapefile.MultiPointM, - coordinates=Array{Array{Float64,1},1}[Array{Float64,1}[[1.15,2.25],[2.15,3.25],[3.15,4.25],[4.15,5.25]],Array{Float64,1}[[11.15,12.25],[12.15,13.25],[13.15,14.25],[14.15,15.25]],Array{Float64,1}[[21.15,22.25],[22.15,23.25],[23.15,24.25],[24.15,25.25]]], - extent=Extent(X=(1.15, 24.15), Y=(2.25, 25.25)), + coordinates=[[[1.15,2.25],[2.15,3.25],[3.15,4.25],[4.15,5.25]],[[11.15,12.25],[12.15,13.25],[13.15,14.25],[14.15,15.25]],[[21.15,22.25],[22.15,23.25],[23.15,24.25],[24.15,25.25]]], + extent=Extent(X=(1.15, 24.15), Y=(2.25, 25.25), Z=(0.0, 0.0)), ),( path=joinpath("shapelib_testcases", "test7.shp"), geomtype=Shapefile.Polyline, - coordinates=Array{Array{Array{Float64,1},1},1}[Array{Array{Float64,1},1}[Array{Float64,1}[[1.0,1.0],[2.0,1.0],[2.0,2.0],[1.0,2.0],[1.0,1.0]]],Array{Array{Float64,1},1}[Array{Float64,1}[[1.0,4.0],[2.0,4.0],[2.0,5.0],[1.0,5.0],[1.0,4.0]]],Array{Array{Float64,1},1}[Array{Float64,1}[[1.0,7.0],[2.0,7.0],[2.0,8.0],[1.0,8.0],[1.0,7.0]]],Array{Array{Float64,1},1}[Array{Float64,1}[[0.0,0.0],[0.0,100.0],[100.0,100.0],[100.0,0.0],[0.0,0.0]],Array{Float64,1}[[10.0,20.0],[30.0,20.0],[30.0,40.0],[10.0,40.0],[10.0,20.0]],Array{Float64,1}[[60.0,20.0],[90.0,20.0],[90.0,40.0],[60.0,40.0],[60.0,20.0]]]], - extent=Extent(X=(0.0, 100.0), Y=(0.0, 100.0)), + coordinates=[[[[1.0,1.0],[2.0,1.0],[2.0,2.0],[1.0,2.0],[1.0,1.0]]],[[[1.0,4.0],[2.0,4.0],[2.0,5.0],[1.0,5.0],[1.0,4.0]]],[[[1.0,7.0],[2.0,7.0],[2.0,8.0],[1.0,8.0],[1.0,7.0]]],[[[0.0,0.0],[0.0,100.0],[100.0,100.0],[100.0,0.0],[0.0,0.0]],Array{Float64,1}[[10.0,20.0],[30.0,20.0],[30.0,40.0],[10.0,40.0],[10.0,20.0]],Array{Float64,1}[[60.0,20.0],[90.0,20.0],[90.0,40.0],[60.0,40.0],[60.0,20.0]]]], + extent=Extent(X=(0.0, 100.0), Y=(0.0, 100.0), Z=(0.0, 0.0)), ),( path=joinpath("shapelib_testcases", "test8.shp"), geomtype=Shapefile.PolylineZ, - coordinates=Array{Array{Array{Float64,1},1},1}[Array{Array{Float64,1},1}[Array{Float64,1}[[1.0,1.0],[2.0,1.0],[2.0,2.0],[1.0,2.0],[1.0,1.0]]],Array{Array{Float64,1},1}[Array{Float64,1}[[1.0,4.0],[2.0,4.0],[2.0,5.0],[1.0,5.0],[1.0,4.0]]],Array{Array{Float64,1},1}[Array{Float64,1}[[1.0,7.0],[2.0,7.0],[2.0,8.0],[1.0,8.0],[1.0,7.0]]],Array{Array{Float64,1},1}[Array{Float64,1}[[0.0,0.0],[0.0,100.0],[100.0,100.0],[100.0,0.0],[0.0,0.0]],Array{Float64,1}[[10.0,20.0],[30.0,20.0],[30.0,40.0],[10.0,40.0],[10.0,20.0]],Array{Float64,1}[[60.0,20.0],[90.0,20.0],[90.0,40.0],[60.0,40.0],[60.0,20.0]]]], - extent=Extent(X=(0.0, 100.0), Y=(0.0, 100.0)), + coordinates=[[[[1.0,1.0,4.45],[2.0,1.0,5.45],[2.0,2.0,6.45],[1.0,2.0,7.45],[1.0,1.0,8.45]]],[[[1.0,4.0,14.45],[2.0,4.0,15.45],[2.0,5.0,16.45],[1.0,5.0,17.45],[1.0,4.0,18.45]]],[[[1.0,7.0,24.45],[2.0,7.0,25.45],[2.0,8.0,26.45],[1.0,8.0,27.45],[1.0,7.0,28.45]]],[[[0.0,0.0,0.0],[0.0,100.0,2.0],[100.0,100.0,4.0],[100.0,0.0,6.0],[0.0,0.0,8.0]],[[10.0,20.0,10.0],[30.0,20.0,12.0],[30.0,40.0,14.0],[10.0,40.0,16.0],[10.0,20.0,18.0]],[[60.0,20.0,20.0],[90.0,20.0,22.0],[90.0,40.0,24.0],[60.0,40.0,26.0],[60.0,20.0,28.0]]]], + extent=Extent(X=(0.0, 100.0), Y=(0.0, 100.0), Z=(0.0, 27.35)), ),( path=joinpath("shapelib_testcases", "test9.shp"), geomtype=Shapefile.PolylineM, - coordinates=Array{Array{Array{Float64,1},1},1}[Array{Array{Float64,1},1}[Array{Float64,1}[[1.0,1.0],[2.0,1.0],[2.0,2.0],[1.0,2.0],[1.0,1.0]]],Array{Array{Float64,1},1}[Array{Float64,1}[[1.0,4.0],[2.0,4.0],[2.0,5.0],[1.0,5.0],[1.0,4.0]]],Array{Array{Float64,1},1}[Array{Float64,1}[[1.0,7.0],[2.0,7.0],[2.0,8.0],[1.0,8.0],[1.0,7.0]]],Array{Array{Float64,1},1}[Array{Float64,1}[[0.0,0.0],[0.0,100.0],[100.0,100.0],[100.0,0.0],[0.0,0.0]],Array{Float64,1}[[10.0,20.0],[30.0,20.0],[30.0,40.0],[10.0,40.0],[10.0,20.0]],Array{Float64,1}[[60.0,20.0],[90.0,20.0],[90.0,40.0],[60.0,40.0],[60.0,20.0]]]], - extent=Extent(X=(0.0, 100.0), Y=(0.0, 100.0)), + coordinates=[[[[1.0,1.0],[2.0,1.0],[2.0,2.0],[1.0,2.0],[1.0,1.0]]],[[[1.0,4.0],[2.0,4.0],[2.0,5.0],[1.0,5.0],[1.0,4.0]]],[[[1.0,7.0],[2.0,7.0],[2.0,8.0],[1.0,8.0],[1.0,7.0]]],[[[0.0,0.0],[0.0,100.0],[100.0,100.0],[100.0,0.0],[0.0,0.0]],[[10.0,20.0],[30.0,20.0],[30.0,40.0],[10.0,40.0],[10.0,20.0]],Array{Float64,1}[[60.0,20.0],[90.0,20.0],[90.0,40.0],[60.0,40.0],[60.0,20.0]]]], + extent=Extent(X=(0.0, 100.0), Y=(0.0, 100.0), Z=(0.0, 0.0)), ),( path=joinpath("shapelib_testcases", "test10.shp"), geomtype=Shapefile.Polygon, - coordinates=[[[[[1.0, 1.0], [2.0, 1.0], [2.0, 2.0], [1.0, 2.0], [1.0, 1.0]]]], [[[[1.0, 4.0], [2.0, 4.0], [2.0, 5.0], [1.0, 5.0], [1.0, 4.0]]]], [[[[1.0, 7.0], [2.0, 7.0], [2.0, 8.0], [1.0, 8.0], [1.0, 7.0]]]], [[[[0.0, 0.0], [0.0, 100.0], [100.0, 100.0], [100.0, 0.0], [0.0, 0.0]], [[10.0, 20.0], [30.0, 20.0], [30.0, 40.0], [10.0, 40.0], [10.0, 20.0]], [[60.0, 20.0], [90.0, 20.0], [90.0, 40.0], [60.0, 40.0], [60.0, 20.0]]]]], - extent=Extent(X=(0.0, 100.0), Y=(0.0, 100.0)), + coordinates=[[[[[1.0,1.0],[2.0,1.0],[2.0,2.0],[1.0,2.0],[1.0,1.0]]]],[[[[1.0,4.0],[2.0,4.0],[2.0,5.0],[1.0,5.0],[1.0,4.0]]]],[[[[1.0,7.0],[2.0,7.0],[2.0,8.0],[1.0,8.0],[1.0,7.0]]]],[[[[0.0,0.0],[0.0,100.0],[100.0,100.0],[100.0,0.0],[0.0,0.0]],[[10.0,20.0],[30.0,20.0],[30.0,40.0],[10.0,40.0],[10.0,20.0]],[[60.0,20.0],[90.0,20.0],[90.0,40.0],[60.0,40.0],[60.0,20.0]]]]], + extent=Extent(X=(0.0, 100.0), Y=(0.0, 100.0), Z=(0.0, 0.0)), ),( path=joinpath("shapelib_testcases", "test11.shp"), geomtype=Shapefile.PolygonZ, - coordinates=[[[[[1.0, 1.0], [2.0, 1.0], [2.0, 2.0], [1.0, 2.0], [1.0, 1.0]]]], [[[[1.0, 4.0], [2.0, 4.0], [2.0, 5.0], [1.0, 5.0], [1.0, 4.0]]]], [[[[1.0, 7.0], [2.0, 7.0], [2.0, 8.0], [1.0, 8.0], [1.0, 7.0]]]], [[[[0.0, 0.0], [0.0, 100.0], [100.0, 100.0], [100.0, 0.0], [0.0, 0.0]], [[10.0, 20.0], [30.0, 20.0], [30.0, 40.0], [10.0, 40.0], [10.0, 20.0]], [[60.0, 20.0], [90.0, 20.0], [90.0, 40.0], [60.0, 40.0], [60.0, 20.0]]]]], - extent=Extent(X=(0.0, 100.0), Y=(0.0, 100.0)), + coordinates=[[[[[1.0,1.0,3.35],[2.0,1.0,4.35],[2.0,2.0,5.35],[1.0,2.0,6.35],[1.0,1.0,7.35]]]],[[[[1.0,4.0,13.35],[2.0,4.0,14.35],[2.0,5.0,15.35],[1.0,5.0,16.35],[1.0,4.0,17.35]]]],[[[[1.0,7.0,23.35],[2.0,7.0,24.35],[2.0,8.0,25.35],[1.0,8.0,26.35],[1.0,7.0,27.35]]]],[[[[0.0,0.0,0.0],[0.0,100.0,1.0],[100.0,100.0,2.0],[100.0,0.0,3.0],[0.0,0.0,4.0]],[[10.0,20.0,5.0],[30.0,20.0,6.0],[30.0,40.0,7.0],[10.0,40.0,8.0],[10.0,20.0,9.0]],[[60.0,20.0,10.0],[90.0,20.0,11.0],[90.0,40.0,12.0],[60.0,40.0,13.0],[60.0,20.0,14.0]]]]], + extent=Extent(X=(0.0, 100.0), Y=(0.0, 100.0), Z=(0.0, 27.35)), ),( path=joinpath("shapelib_testcases", "test12.shp"), geomtype=Shapefile.PolygonM, coordinates=[[[[[1.0, 1.0], [2.0, 1.0], [2.0, 2.0], [1.0, 2.0], [1.0, 1.0]]]], [[[[1.0, 4.0], [2.0, 4.0], [2.0, 5.0], [1.0, 5.0], [1.0, 4.0]]]], [[[[1.0, 7.0], [2.0, 7.0], [2.0, 8.0], [1.0, 8.0], [1.0, 7.0]]]], [[[[0.0, 0.0], [0.0, 100.0], [100.0, 100.0], [100.0, 0.0], [0.0, 0.0]], [[10.0, 20.0], [30.0, 20.0], [30.0, 40.0], [10.0, 40.0], [10.0, 20.0]], [[60.0, 20.0], [90.0, 20.0], [90.0, 40.0], [60.0, 40.0], [60.0, 20.0]]]]], - extent=Extent(X=(0.0, 100.0), Y=(0.0, 100.0)), + extent=Extent(X=(0.0, 100.0), Y=(0.0, 100.0), Z=(0.0, 0.0)), ),( path=joinpath("shapelib_testcases/test13.shp"), geomtype=Shapefile.MultiPatch, coordinates=nothing, - extent=Extent(X=(0.0, 100.0), Y=(0.0, 100.0)), + extent=Extent(X=(0.0, 100.0), Y=(0.0, 100.0), Z=(0.0, 27.35)), ) ] @testset "Shapefile" begin +for t in test_tuples + if !(t.geomtype <: Union{Missing,Shapefile.MultiPatch}) + sh = Shapefile.Handle(t.path) + p = sh.shapes + display(plot(p)) + sleep(1) + end +end + for test in test_tuples for use_shx in (false, true) shp = if use_shx From d8a73fe83d2189684a8fb6945a22cb5c3fe8bf9d Mon Sep 17 00:00:00 2001 From: Maarten Pronk Date: Wed, 18 May 2022 08:49:45 +0200 Subject: [PATCH 03/16] Updated to latest GeoInterface. --- src/multipoints.jl | 14 +++++++------- src/points.jl | 5 ++--- src/polygons.jl | 19 +++++++++---------- src/polylines.jl | 28 ++++++++++++++-------------- 4 files changed, 32 insertions(+), 34 deletions(-) diff --git a/src/multipoints.jl b/src/multipoints.jl index 56d3dbd..32e89df 100644 --- a/src/multipoints.jl +++ b/src/multipoints.jl @@ -1,9 +1,9 @@ abstract type AbstractMultiPoint{T} <: AbstractShape end -GI.geomtype(::AbstractMultiPoint) = GI.MultiPointTrait() +GI.geomtrait(::AbstractMultiPoint) = GI.MultiPointTrait() GI.ngeom(::GI.MultiPointTrait, geom::AbstractMultiPoint) = length(geom.points) -GI.ncoord(::GI.MultiPointTrait, geom::AbstractMultiPoint{T}) where T = _ncoord(T) +GI.ncoord(::GI.MultiPointTrait, geom::AbstractMultiPoint{T}) where {T} = _ncoord(T) GI.npoint(::GI.MultiPointTrait, geom::AbstractMultiPoint) = length(geom.points) """ @@ -12,7 +12,7 @@ GI.npoint(::GI.MultiPointTrait, geom::AbstractMultiPoint) = length(geom.points) Collection of points, from a shape file. # Fields -- `points`: a `Vector` of [`Point`](@ref). +- `points`: a `Vector` of [`Point`](@ref). - `MBR`: `nothing` or the known bounding box. Can be retrieved with `GeoInterface.bbox`. """ struct MultiPoint <: AbstractMultiPoint{Point} @@ -33,14 +33,14 @@ GI.getgeom(::GI.MultiPointTrait, geom::MultiPoint, i::Integer) = geom.points[i] """ MultiPointM <: AbstractMultiPoint -Collection of points, from a shape file. +Collection of points, from a shape file. Includes a `measures` field, holding values from each point. May have a known bounding box, which can be retrieved with `GeoInterface.bbox`. # Fields -- `points`: a `Vector` of [`Point`](@ref). +- `points`: a `Vector` of [`Point`](@ref). - `MBR`: `nothing` or the known bounding box. Can be retrieved with `GeoInterface.bbox`. - `measures`: holds values from each point. """ @@ -68,14 +68,14 @@ GI.getgeom(::GI.MultiPointTrait, geom::MultiPointM, i::Integer) = """ MultiPointZ <: AbstractMultiPoint -Collection of 3d points, from a shape file. +Collection of 3d points, from a shape file. Includes a `measures` field, holding values from each point. May have a known bounding box, which can be retrieved with `GeoInterface.bbox`. # Fields -- `points`: a `Vector` of [`Point`](@ref). +- `points`: a `Vector` of [`Point`](@ref). - `zvalues`: a `Vector` of `Float64` representing the z dimension values. - `MBR`: `nothing` or the known bounding box. Can be retrieved with `GeoInterface.bbox`. - `measures`: holds values from each point. diff --git a/src/points.jl b/src/points.jl index f050368..fb324ab 100644 --- a/src/points.jl +++ b/src/points.jl @@ -1,11 +1,11 @@ abstract type AbstractPoint <: AbstractShape end -GI.geomtype(::AbstractPoint) = GI.PointTrait() +GI.geomtrait(::AbstractPoint) = GI.PointTrait() GI.x(::GI.PointTrait, point::AbstractPoint) = point.x GI.y(::GI.PointTrait, point::AbstractPoint) = point.y GI.getcoord(::GI.PointTrait, p::AbstractPoint, i::Integer) = getfield(p, i) -GI.ncoord(::GI.PointTrait, p::T) where T <: AbstractPoint = _ncoord(T) +GI.ncoord(::GI.PointTrait, p::T) where {T<:AbstractPoint} = _ncoord(T) _ncoord(::Type{<:AbstractPoint}) = 2 @@ -81,4 +81,3 @@ GI.m(::GI.PointTrait, point::PointZ) = point.m GI.z(::GI.PointTrait, point::PointZ) = point.z _ncoord(::Type{PointZ}) = 3 - diff --git a/src/polygons.jl b/src/polygons.jl index 57b567b..f5b6696 100644 --- a/src/polygons.jl +++ b/src/polygons.jl @@ -12,11 +12,11 @@ struct LinearRing{P,Z,M} <: AbstractVector{P} m::M end LinearRing{P}(xy, z::Z, m::M) where {P,Z,M} = LinearRing{P,Z,M}(xy, z, m) -LinearRing{P}(xy; z=nothing, m=nothing) where P = LinearRing{P}(xy, z, m) +LinearRing{P}(xy; z=nothing, m=nothing) where {P} = LinearRing{P}(xy, z, m) GI.isgeometry(::Type{<:LinearRing}) = true -GI.geomtype(::LinearRing) = GI.LinearRingTrait() -GI.ncoord(lr::LinearRing{P}) where P = _ncoord(P) +GI.geomtrait(::LinearRing) = GI.LinearRingTrait() +GI.ncoord(lr::LinearRing{P}) where {P} = _ncoord(P) GI.ngeom(::GI.LinearRingTrait, lr::LinearRing) = length(lr) # GI.getgeom(::GI.LinearRingTrait, lr::LinearRing) = (getindex(lr, i) for i in 1:ngeom(lr)) @@ -35,8 +35,8 @@ struct SubPolygon{L<:LinearRing} <: AbstractVector{L} rings::Vector{L} end GI.isgeometry(::Type{<:SubPolygon}) = true -GI.geomtype(::SubPolygon) = GI.PolygonTrait() -GI.ncoord(::GI.PolygonTrait, ::SubPolygon{LinearRing{P}}) where P = _ncoord(P) +GI.geomtrait(::SubPolygon) = GI.PolygonTrait() +GI.ncoord(::GI.PolygonTrait, ::SubPolygon{LinearRing{P}}) where {P} = _ncoord(P) GI.ngeom(::GI.PolygonTrait, sp::SubPolygon) = length(sp) GI.getgeom(::GI.PolygonTrait, sp::SubPolygon, i::Integer) = getindex(sp, i) @@ -50,7 +50,7 @@ Base.push!(p::SubPolygon, x) = Base.push!(parent(p), x) abstract type AbstractPolygon{T} <: AbstractShape end # Shapefile polygons are OGC multipolygons -GI.geomtype(geom::AbstractPolygon) = GI.MultiPolygonTrait() +GI.geomtrait(geom::AbstractPolygon) = GI.MultiPolygonTrait() GI.nring(::GI.MultiPolygonTrait, geom::AbstractPolygon) = length(geom.parts) function GI.ngeom(::GI.MultiPolygonTrait, geom::AbstractPolygon) n = 0 @@ -66,7 +66,7 @@ function _isclockwise(ring) clockwise_test = 0.0 for i in 1:GI.npoint(ring)-1 prev = ring[i] - cur = ring[i + 1] + cur = ring[i+1] clockwise_test += (cur.x - prev.x) * (cur.y + prev.y) end clockwise_test > 0 @@ -75,7 +75,7 @@ end function GI.getring(t::GI.MultiPolygonTrait, geom::AbstractPolygon) return (GI.getring(t, geom, i) for i in eachindex(geom.parts)) end -function GI.getring(::GI.MultiPolygonTrait, geom::AbstractPolygon{P}, i::Integer) where P +function GI.getring(::GI.MultiPolygonTrait, geom::AbstractPolygon{P}, i::Integer) where {P} pa = geom.parts range = pa[i]+1:(i == lastindex(pa) ? lastindex(geom.points) : pa[i+1]) xy = view(geom.points, range) @@ -90,7 +90,7 @@ end # Warning: getgeom is very slow for a Shapefile. Use `getpolygons`, # or if you don't need exteriors and holes to be separated, `getring`. GI.getgeom(::GI.MultiPolygonTrait, geom::AbstractPolygon, i::Integer) = collect(GI.getgeom(geom))[i] -function GI.getgeom(::GI.MultiPolygonTrait, geom::AbstractPolygon{T}) where T +function GI.getgeom(::GI.MultiPolygonTrait, geom::AbstractPolygon{T}) where {T} r1 = GI.getring(geom, 1) polygons = SubPolygon{typeof(r1)}[] holes = typeof(r1)[] @@ -230,4 +230,3 @@ function Base.read(io::IO, ::Type{PolygonZ}) read!(io, measures) PolygonZ(box, parts, points, zrange, zvalues, mrange, measures) end - diff --git a/src/polylines.jl b/src/polylines.jl index 49fa9ee..493cd97 100644 --- a/src/polylines.jl +++ b/src/polylines.jl @@ -5,11 +5,11 @@ struct LineString{P,Z,M} <: AbstractVector{P} m::M end LineString{P}(xy, z::Z, m::M) where {P,Z,M} = LineString{P,Z,M}(xy, z, m) -LineString{P}(xy; z=nothing, m=nothing) where P = LineString{P}(xy, z, m) +LineString{P}(xy; z=nothing, m=nothing) where {P} = LineString{P}(xy, z, m) GI.isgeometry(::Type{<:LineString}) = true -GI.geomtype(::LineString) = GI.LineStringTrait() -GI.ncoord(lr::LineString{T}) where T = _ncoord(T) +GI.geomtrait(::LineString) = GI.LineStringTrait() +GI.ncoord(lr::LineString{T}) where {T} = _ncoord(T) GI.ngeom(::GI.LineStringTrait, lr::LineString) = length(lr) # GI.getgeom(::GI.LineStringTrait, lr::LineString) = (getindex(lr, i) for i in 1:ngeom(lr)) GI.getgeom(::GI.LineStringTrait, lr::LineString, i::Integer) = getindex(lr, i) @@ -24,17 +24,17 @@ Base.getindex(lr::LineString{PointZ}, i) = PointZ(lr.xy[i], lr.m[i], lr.z[i]) abstract type AbstractPolyline{T} <: AbstractShape end -GI.geomtype(::AbstractPolyline) = GI.MultiLineStringTrait() +GI.geomtrait(::AbstractPolyline) = GI.MultiLineStringTrait() GI.ngeom(::GI.MultiLineStringTrait, geom::AbstractPolyline) = length(geom.parts) -GI.ncoord(::GI.MultiLineStringTrait, ::AbstractPolyline{T}) where T = _ncoord(T) -function GI.getgeom(::GI.MultiLineStringTrait, geom::AbstractPolyline{P}, i::Integer) where P +GI.ncoord(::GI.MultiLineStringTrait, ::AbstractPolyline{T}) where {T} = _ncoord(T) +function GI.getgeom(::GI.MultiLineStringTrait, geom::AbstractPolyline{P}, i::Integer) where {P} range = if i < GI.ngeom(geom) # The linestring is the points between two parts # C integers start at zero so we add 1 to the range start - (geom.parts[i] + 1):(geom.parts[i+1]) + (geom.parts[i]+1):(geom.parts[i+1]) else # For the last inestring use the vector lastindex as the last point - (geom.parts[i] + 1):lastindex(geom.points) + (geom.parts[i]+1):lastindex(geom.points) end return LineString(geom, range) end @@ -45,7 +45,7 @@ end Represents a single or multiple polylines from a shape file. # Fields -- `points`: a `Vector` of [`Point`]() represents a one or multiple lines. +- `points`: a `Vector` of [`Point`]() represents a one or multiple lines. - `parts`: a `Vector` of `Int32` indicating the line each point belongs to. - `MBR`: `nothing` or the known bounding box. Can be retrieved with `GeoInterface.bbox`. """ @@ -74,7 +74,7 @@ LineString(geom::Polyline, range) = @views LineString{Point}(geom.points[range]) Polyline from a shape file, with measures. # Fields -- `points`: a `Vector` of [`Point`](@ref) represents a one or multiple lines. +- `points`: a `Vector` of [`Point`](@ref) represents a one or multiple lines. - `parts`: a `Vector` of `Int32` indicating the line each point belongs to. - `MBR`: `nothing` or the known bounding box. Can be retrieved with `GeoInterface.bbox`. - `measures`: holds values from each point. @@ -87,7 +87,7 @@ struct PolylineM <: AbstractPolyline{PointM} measures::Vector{Float64} end -LineString(geom::PolylineM, range) = +LineString(geom::PolylineM, range) = @views LineString{PointM}(geom.points[range]; m=geom.measures[range]) function Base.read(io::IO, ::Type{PolylineM}) @@ -107,10 +107,10 @@ end """ PolylineZ <: AbstractPolyline -Three dimensional polyline of from a shape file. +Three dimensional polyline of from a shape file. # Fields -- `points`: a `Vector` of [`Point`](@ref) represents a one or multiple lines. +- `points`: a `Vector` of [`Point`](@ref) represents a one or multiple lines. - `parts`: a `Vector` of `Int32` indicating the line each point belongs to. - `zvalues`: a `Vector` of `Float64` representing the z dimension values. - `MBR`: `nothing` or the known bounding box. Can be retrieved with `GeoInterface.bbox`. @@ -126,7 +126,7 @@ struct PolylineZ <: AbstractPolyline{PointZ} measures::Vector{Float64} end -LineString(geom::PolylineZ, range) = +LineString(geom::PolylineZ, range) = @views LineString{PointZ}(geom.points[range]; z=geom.zvalues[range], m=geom.measures[range]) function Base.read(io::IO, ::Type{PolylineZ}) From 329cfcf3d1267251d43c2c2614dd31632b120f48 Mon Sep 17 00:00:00 2001 From: rafaqz Date: Tue, 10 May 2022 00:40:58 +0400 Subject: [PATCH 04/16] minor cleanup --- src/polygons.jl | 21 +++++++++++---------- test/runtests.jl | 20 +++++++++++--------- 2 files changed, 22 insertions(+), 19 deletions(-) diff --git a/src/polygons.jl b/src/polygons.jl index f5b6696..247041e 100644 --- a/src/polygons.jl +++ b/src/polygons.jl @@ -52,6 +52,7 @@ abstract type AbstractPolygon{T} <: AbstractShape end # Shapefile polygons are OGC multipolygons GI.geomtrait(geom::AbstractPolygon) = GI.MultiPolygonTrait() GI.nring(::GI.MultiPolygonTrait, geom::AbstractPolygon) = length(geom.parts) +GI.npoint(::GI.MultiPolygonTrait, geom::AbstractPolygon) = length(geom.points) function GI.ngeom(::GI.MultiPolygonTrait, geom::AbstractPolygon) n = 0 for ring in GI.getring(geom) @@ -62,16 +63,6 @@ function GI.ngeom(::GI.MultiPolygonTrait, geom::AbstractPolygon) return n end -function _isclockwise(ring) - clockwise_test = 0.0 - for i in 1:GI.npoint(ring)-1 - prev = ring[i] - cur = ring[i+1] - clockwise_test += (cur.x - prev.x) * (cur.y + prev.y) - end - clockwise_test > 0 -end - function GI.getring(t::GI.MultiPolygonTrait, geom::AbstractPolygon) return (GI.getring(t, geom, i) for i in eachindex(geom.parts)) end @@ -129,6 +120,16 @@ function _inring(pt::Point, ring::LinearRing) return isinside end +function _isclockwise(ring) + clockwise_test = 0.0 + for i in 1:GI.npoint(ring)-1 + prev = ring[i] + cur = ring[i + 1] + clockwise_test += (cur.x - prev.x) * (cur.y + prev.y) + end + clockwise_test > 0 +end + """ Polygon <: AbstractPolygon diff --git a/test/runtests.jl b/test/runtests.jl index a0c7c5a..8588a62 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -83,16 +83,18 @@ test_tuples = [ ) ] -@testset "Shapefile" begin +# Visual plot check +# for t in test_tuples +# if !(t.geomtype <: Union{Missing,Shapefile.MultiPatch}) +# @show t.path t.geomtype +# sh = Shapefile.Handle(t.path) +# p = sh.shapes +# display(plot(p; opacity=.5)) +# sleep(1) +# end +# end -for t in test_tuples - if !(t.geomtype <: Union{Missing,Shapefile.MultiPatch}) - sh = Shapefile.Handle(t.path) - p = sh.shapes - display(plot(p)) - sleep(1) - end -end +@testset "Shapefile" begin for test in test_tuples for use_shx in (false, true) From b0e566642b9cb463aa6880bc0e46217c250a0ab9 Mon Sep 17 00:00:00 2001 From: rafaqz Date: Wed, 18 May 2022 18:36:30 +0400 Subject: [PATCH 05/16] add interface tests --- src/common.jl | 2 +- src/multipatch.jl | 19 +++++++++++--- test/runtests.jl | 65 ++++++++++++++++++++++++++++++++++------------- test/table.jl | 2 +- 4 files changed, 65 insertions(+), 23 deletions(-) diff --git a/src/common.jl b/src/common.jl index e52c75a..667cb23 100644 --- a/src/common.jl +++ b/src/common.jl @@ -1,7 +1,7 @@ abstract type AbstractShape end -isgeometry(::AbstractShape) = true +GeoInterface.isgeometry(::Type{<:AbstractShape}) = true GeoInterface.ncoord(::GI.AbstractGeometryTrait, ::AbstractShape) = 2 # With specific methods when 3 """ diff --git a/src/multipatch.jl b/src/multipatch.jl index 8c7e7bb..2795dcb 100644 --- a/src/multipatch.jl +++ b/src/multipatch.jl @@ -4,11 +4,12 @@ Stores a collection of patch representing the boundary of a 3d object. # Fields -- `points`: a `Vector` of [`Point`](@ref) represents a one or multiple spatial objects. +- `MBR`: `nothing` or the known bounding box. Can be retrieved with `GeoInterface.bbox`. - `parts`: a `Vector` of `Int32` indicating the object each point belongs to. - `parttypes`: a `Vector` of `Int32` indicating the type of object each point belongs to. +- `points`: a `Vector` of [`Point`](@ref) represents a one or multiple spatial objects. +- `zrange`: and `Interval` of bounds for the `zvalues`. - `zvalues`: a `Vector` of `Float64` indicating absolute or relative heights. -- `MBR`: `nothing` or the known bounding box. Can be retrieved with `GeoInterface.bbox`. """ struct MultiPatch <: AbstractShape MBR::Rect @@ -20,9 +21,21 @@ struct MultiPatch <: AbstractShape # measures::Vector{Float64} # (optional) end -GeoInterface.geomtype(::MultiPatch) = GeoInterface.MultiPolygonTrait() +GI.geomtrait(::MultiPatch) = GI.GeometryCollectionTrait() +GI.ngeom(geom::MultiPatch) = length(geom.parts) GeoInterface.ncoord(::MultiPatch) = 3 +# TODO implement the rest of this +function GI.getgeom(::GI.GeometryCollectionTrait, geom::MultiPatch, i::Integer) + error("`getgeom` not implemented for `MultiPatch`: open a github issue at JuliaGeo/Shapefile.jl if you need this.") + # T = SHAPETYPE[geom.parttypes[i]] +end +function GI.subtrait(geom::MultiPatch) + error("`subtype` not implemented for `MultiPatch`: open a github issue at JuliaGeo/Shapefile.jl if you need this.") + # types = map(i -> SHAPETYPE[i], union(geom.parttypes)) +end + + function Base.read(io::IO, ::Type{MultiPatch}) box = read(io, Rect) numparts = read(io, Int32) diff --git a/test/runtests.jl b/test/runtests.jl index 8588a62..b364355 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -2,12 +2,16 @@ using Shapefile, GeoInterface, Plots, Extents using Shapefile using Test +using Shapefile: Point, PointM, PointZ, Polygon, PolygonM, PolygonZ, Polyline, + PolylineM, PolylineZ, MultiPoint, MultiPointM, MultiPointZ, + MultiPatch, LineString, LinearRing, SubPolygon, Rect, Interval + shp = Shapefile.Handle(joinpath("shapelib_testcases", "test.shp")) test_tuples = [ ( path=joinpath("shapelib_testcases", "test.shp"), - geomtype=Shapefile.Polygon, + geomtype=Polygon, coordinates=[[[[[20.0,20.0],[20.0,30.0],[30.0,30.0],[20.0,20.0]]],[[[0.0,0.0],[100.0,0.0],[100.0,100.0],[0.0,100.0],[0.0,0.0]]]],[[[[150.0,150.0],[160.0,150.0],[180.0,170.0],[150.0,150.0]]]],[[[[150.0,150.0],[160.0,150.0],[180.0,170.0],[150.0,150.0]]]]], extent=Extent(X=(0.0, 180.0), Y=(0.0, 170.0), Z=(0.0, 0.0)), ),( @@ -17,67 +21,67 @@ test_tuples = [ extent=Extent(X=(0.0, 10.0), Y=(0.0, 20.0), Z=(0.0, 0.0)), ),( path=joinpath("shapelib_testcases", "test1.shp"), - geomtype=Shapefile.Point, + geomtype=Point, coordinates=[[1.0,2.0],[10.0,20.0]], extent=Extent(X=(1.0, 10.0), Y=(2.0, 20.0), Z=(0.0, 0.0)), ),( path=joinpath("shapelib_testcases", "test2.shp"), - geomtype=Shapefile.PointZ, + geomtype=PointZ, coordinates=[[1.0,2.0,3.0],[10.0,20.0,30.0]], extent=Extent(X=(1.0, 10.0), Y=(2.0, 20.0), Z=(3.0, 30.0)), ),( path=joinpath("shapelib_testcases", "test3.shp"), - geomtype=Shapefile.PointM, + geomtype=PointM, coordinates=[[1.0,2.0],[10.0,20.0]], extent=Extent(X=(1.0, 10.0), Y=(2.0, 20.0), Z=(0.0, 0.0)), ),( path=joinpath("shapelib_testcases", "test4.shp"), - geomtype=Shapefile.MultiPoint, + geomtype=MultiPoint, coordinates=[[[1.15,2.25],[2.15,3.25],[3.15,4.25],[4.15,5.25]],[[11.15,12.25],[12.15,13.25],[13.15,14.25],[14.15,15.25]],[[21.15,22.25],[22.15,23.25],[23.15,24.25],[24.15,25.25]]], extent=Extent(X=(1.15, 24.15), Y=(2.25, 25.25), Z=(0.0, 0.0)), ),( path=joinpath("shapelib_testcases", "test5.shp"), - geomtype=Shapefile.MultiPointZ, + geomtype=MultiPointZ, coordinates=[[[1.15,2.25,3.35],[2.15,3.25,4.35],[3.15,4.25,5.35],[4.15,5.25,6.35]],[[11.15,12.25,13.35],[12.15,13.25,14.35],[13.15,14.25,15.35],[14.15,15.25,16.35]],[[21.15,22.25,23.35],[22.15,23.25,24.35],[23.15,24.25,25.35],[24.15,25.25,26.35]]], extent=Extent(X=(1.15, 24.15), Y=(2.25, 25.25), Z=(3.35, 26.35)), ),( path=joinpath("shapelib_testcases", "test6.shp"), - geomtype=Shapefile.MultiPointM, + geomtype=MultiPointM, coordinates=[[[1.15,2.25],[2.15,3.25],[3.15,4.25],[4.15,5.25]],[[11.15,12.25],[12.15,13.25],[13.15,14.25],[14.15,15.25]],[[21.15,22.25],[22.15,23.25],[23.15,24.25],[24.15,25.25]]], extent=Extent(X=(1.15, 24.15), Y=(2.25, 25.25), Z=(0.0, 0.0)), ),( path=joinpath("shapelib_testcases", "test7.shp"), - geomtype=Shapefile.Polyline, + geomtype=Polyline, coordinates=[[[[1.0,1.0],[2.0,1.0],[2.0,2.0],[1.0,2.0],[1.0,1.0]]],[[[1.0,4.0],[2.0,4.0],[2.0,5.0],[1.0,5.0],[1.0,4.0]]],[[[1.0,7.0],[2.0,7.0],[2.0,8.0],[1.0,8.0],[1.0,7.0]]],[[[0.0,0.0],[0.0,100.0],[100.0,100.0],[100.0,0.0],[0.0,0.0]],Array{Float64,1}[[10.0,20.0],[30.0,20.0],[30.0,40.0],[10.0,40.0],[10.0,20.0]],Array{Float64,1}[[60.0,20.0],[90.0,20.0],[90.0,40.0],[60.0,40.0],[60.0,20.0]]]], extent=Extent(X=(0.0, 100.0), Y=(0.0, 100.0), Z=(0.0, 0.0)), ),( path=joinpath("shapelib_testcases", "test8.shp"), - geomtype=Shapefile.PolylineZ, + geomtype=PolylineZ, coordinates=[[[[1.0,1.0,4.45],[2.0,1.0,5.45],[2.0,2.0,6.45],[1.0,2.0,7.45],[1.0,1.0,8.45]]],[[[1.0,4.0,14.45],[2.0,4.0,15.45],[2.0,5.0,16.45],[1.0,5.0,17.45],[1.0,4.0,18.45]]],[[[1.0,7.0,24.45],[2.0,7.0,25.45],[2.0,8.0,26.45],[1.0,8.0,27.45],[1.0,7.0,28.45]]],[[[0.0,0.0,0.0],[0.0,100.0,2.0],[100.0,100.0,4.0],[100.0,0.0,6.0],[0.0,0.0,8.0]],[[10.0,20.0,10.0],[30.0,20.0,12.0],[30.0,40.0,14.0],[10.0,40.0,16.0],[10.0,20.0,18.0]],[[60.0,20.0,20.0],[90.0,20.0,22.0],[90.0,40.0,24.0],[60.0,40.0,26.0],[60.0,20.0,28.0]]]], extent=Extent(X=(0.0, 100.0), Y=(0.0, 100.0), Z=(0.0, 27.35)), ),( path=joinpath("shapelib_testcases", "test9.shp"), - geomtype=Shapefile.PolylineM, + geomtype=PolylineM, coordinates=[[[[1.0,1.0],[2.0,1.0],[2.0,2.0],[1.0,2.0],[1.0,1.0]]],[[[1.0,4.0],[2.0,4.0],[2.0,5.0],[1.0,5.0],[1.0,4.0]]],[[[1.0,7.0],[2.0,7.0],[2.0,8.0],[1.0,8.0],[1.0,7.0]]],[[[0.0,0.0],[0.0,100.0],[100.0,100.0],[100.0,0.0],[0.0,0.0]],[[10.0,20.0],[30.0,20.0],[30.0,40.0],[10.0,40.0],[10.0,20.0]],Array{Float64,1}[[60.0,20.0],[90.0,20.0],[90.0,40.0],[60.0,40.0],[60.0,20.0]]]], extent=Extent(X=(0.0, 100.0), Y=(0.0, 100.0), Z=(0.0, 0.0)), ),( path=joinpath("shapelib_testcases", "test10.shp"), - geomtype=Shapefile.Polygon, + geomtype=Polygon, coordinates=[[[[[1.0,1.0],[2.0,1.0],[2.0,2.0],[1.0,2.0],[1.0,1.0]]]],[[[[1.0,4.0],[2.0,4.0],[2.0,5.0],[1.0,5.0],[1.0,4.0]]]],[[[[1.0,7.0],[2.0,7.0],[2.0,8.0],[1.0,8.0],[1.0,7.0]]]],[[[[0.0,0.0],[0.0,100.0],[100.0,100.0],[100.0,0.0],[0.0,0.0]],[[10.0,20.0],[30.0,20.0],[30.0,40.0],[10.0,40.0],[10.0,20.0]],[[60.0,20.0],[90.0,20.0],[90.0,40.0],[60.0,40.0],[60.0,20.0]]]]], extent=Extent(X=(0.0, 100.0), Y=(0.0, 100.0), Z=(0.0, 0.0)), ),( path=joinpath("shapelib_testcases", "test11.shp"), - geomtype=Shapefile.PolygonZ, + geomtype=PolygonZ, coordinates=[[[[[1.0,1.0,3.35],[2.0,1.0,4.35],[2.0,2.0,5.35],[1.0,2.0,6.35],[1.0,1.0,7.35]]]],[[[[1.0,4.0,13.35],[2.0,4.0,14.35],[2.0,5.0,15.35],[1.0,5.0,16.35],[1.0,4.0,17.35]]]],[[[[1.0,7.0,23.35],[2.0,7.0,24.35],[2.0,8.0,25.35],[1.0,8.0,26.35],[1.0,7.0,27.35]]]],[[[[0.0,0.0,0.0],[0.0,100.0,1.0],[100.0,100.0,2.0],[100.0,0.0,3.0],[0.0,0.0,4.0]],[[10.0,20.0,5.0],[30.0,20.0,6.0],[30.0,40.0,7.0],[10.0,40.0,8.0],[10.0,20.0,9.0]],[[60.0,20.0,10.0],[90.0,20.0,11.0],[90.0,40.0,12.0],[60.0,40.0,13.0],[60.0,20.0,14.0]]]]], extent=Extent(X=(0.0, 100.0), Y=(0.0, 100.0), Z=(0.0, 27.35)), ),( path=joinpath("shapelib_testcases", "test12.shp"), - geomtype=Shapefile.PolygonM, + geomtype=PolygonM, coordinates=[[[[[1.0, 1.0], [2.0, 1.0], [2.0, 2.0], [1.0, 2.0], [1.0, 1.0]]]], [[[[1.0, 4.0], [2.0, 4.0], [2.0, 5.0], [1.0, 5.0], [1.0, 4.0]]]], [[[[1.0, 7.0], [2.0, 7.0], [2.0, 8.0], [1.0, 8.0], [1.0, 7.0]]]], [[[[0.0, 0.0], [0.0, 100.0], [100.0, 100.0], [100.0, 0.0], [0.0, 0.0]], [[10.0, 20.0], [30.0, 20.0], [30.0, 40.0], [10.0, 40.0], [10.0, 20.0]], [[60.0, 20.0], [90.0, 20.0], [90.0, 40.0], [60.0, 40.0], [60.0, 20.0]]]]], extent=Extent(X=(0.0, 100.0), Y=(0.0, 100.0), Z=(0.0, 0.0)), ),( path=joinpath("shapelib_testcases/test13.shp"), - geomtype=Shapefile.MultiPatch, + geomtype=MultiPatch, coordinates=nothing, extent=Extent(X=(0.0, 100.0), Y=(0.0, 100.0), Z=(0.0, 27.35)), ) @@ -85,7 +89,7 @@ test_tuples = [ # Visual plot check # for t in test_tuples -# if !(t.geomtype <: Union{Missing,Shapefile.MultiPatch}) +# if !(t.geomtype <: Union{Missing,MultiPatch}) # @show t.path t.geomtype # sh = Shapefile.Handle(t.path) # p = sh.shapes @@ -94,7 +98,32 @@ test_tuples = [ # end # end -@testset "Shapefile" begin +@testset "GeoInterface compatability" begin + points = [Point(1, 2), Point(2, 2), Point(2, 1), Point(1, 1)] + + shapes = ( + Point(1, 2), + PointM(1, 2, 3), + PointZ(1, 2, 3, 4), + MultiPoint(Rect(1, 2, 2, 1), points), + MultiPointM(Rect(1, 2, 2, 1), points, Interval(10.0, 40.0), [10.0, 20.0, 30.0, 40.0]), + MultiPointZ(Rect(1, 2, 2, 1), points, Interval(1, 4), [1, 2, 3, 4], Interval(10.0, 40.0), [10.0, 20.0, 30.0, 40.0]), + Polygon(Rect(1, 2, 2, 1), [1], points), + PolygonM(Rect(1, 2, 2, 1), [1], points, Interval(10.0, 40.0), [10.0, 20.0, 30.0, 40.0]), + PolygonZ(Rect(1, 2, 2, 1), [1], points, Interval(1, 4), [1, 2, 3, 4], Interval(10.0, 40.0), [10.0, 20.0, 30.0, 40.0]), + Polyline(Rect(1, 2, 2, 1), [1], points), + PolylineM(Rect(1, 2, 2, 1), [1], points, Interval(10.0, 40.0), [10.0, 20.0, 30.0, 40.0]), + PolylineZ(Rect(1, 2, 2, 1), [1], points, Interval(1, 4), [1, 2, 3, 4], Interval(10.0, 40.0), [10.0, 20.0, 30.0, 40.0]), + LineString{Point}(view(points, 1:4)), + LinearRing{Point}(view(points, 1:4)), + SubPolygon([LinearRing{Point}(view(points, 1:4))]), + ) + + @test all(s -> GeoInterface.testgeometry(s), shapes) + @test_broken GeoInterface.testgeometry(MultiPatch(Rect(1, 2, 2, 1), [1], [1], points, Interval(1, 4), [1, 2, 3, 4])) +end + +@testset "Loading Shapefiles" begin for test in test_tuples for use_shx in (false, true) @@ -173,6 +202,6 @@ for test in test_tuples end -include("table.jl") +end # @testset "Loading Shapefiles" -end # @testset "Shapefile" +include("table.jl") diff --git a/test/table.jl b/test/table.jl index c9d7792..31a928b 100644 --- a/test/table.jl +++ b/test/table.jl @@ -150,7 +150,7 @@ end @test Shapefile.shape(r) isa Shapefile.Point @test r.featurecla in classes end - show_result = "Shapefile.Table{Union{Missing, Shapefile.Point}} with 243 rows and the following 39 columns:\n\t\ngeometry, scalerank, natscale, labelrank, featurecla, name, namepar, namealt, diffascii, nameascii, adm0cap, capalt, capin, worldcity, megacity, sov0name, sov_a3, adm0name, adm0_a3, adm1name, iso_a2, note, latitude, longitude, changed, namediff, diffnote, pop_max, pop_min, pop_other, rank_max, rank_min, geonameid, meganame, ls_name, ls_match, checkme, min_zoom, ne_id\n" + show_result = "Shapefile.Table{Union{Missing, Point}} with 243 rows and the following 39 columns:\n\t\ngeometry, scalerank, natscale, labelrank, featurecla, name, namepar, namealt, diffascii, nameascii, adm0cap, capalt, capin, worldcity, megacity, sov0name, sov_a3, adm0name, adm0_a3, adm1name, iso_a2, note, latitude, longitude, changed, namediff, diffnote, pop_max, pop_min, pop_other, rank_max, rank_min, geonameid, meganame, ls_name, ls_match, checkme, min_zoom, ne_id\n" if VERSION < v"1.1" show_result = replace(show_result, "Shapefile.Point" => "Point") end From ef569fc803f314246fa64e728749495ec6f2330e Mon Sep 17 00:00:00 2001 From: rafaqz Date: Thu, 19 May 2022 09:10:36 +0400 Subject: [PATCH 06/16] cleanup --- src/polygons.jl | 8 ++------ 1 file changed, 2 insertions(+), 6 deletions(-) diff --git a/src/polygons.jl b/src/polygons.jl index 247041e..5fa4d1c 100644 --- a/src/polygons.jl +++ b/src/polygons.jl @@ -18,7 +18,6 @@ GI.isgeometry(::Type{<:LinearRing}) = true GI.geomtrait(::LinearRing) = GI.LinearRingTrait() GI.ncoord(lr::LinearRing{P}) where {P} = _ncoord(P) GI.ngeom(::GI.LinearRingTrait, lr::LinearRing) = length(lr) -# GI.getgeom(::GI.LinearRingTrait, lr::LinearRing) = (getindex(lr, i) for i in 1:ngeom(lr)) GI.getgeom(::GI.LinearRingTrait, lr::LinearRing{Point}, i::Integer) = lr.xy[i] GI.getgeom(::GI.LinearRingTrait, lr::LinearRing{PointM}, i::Integer) = @@ -75,11 +74,8 @@ function GI.getring(::GI.MultiPolygonTrait, geom::AbstractPolygon{P}, i::Integer LinearRing{P}(xy, z, m) end -# function GI.gethole(::GI.MultiPolygonTrait, geom::AbstractPolygon) -# end - -# Warning: getgeom is very slow for a Shapefile. Use `getpolygons`, -# or if you don't need exteriors and holes to be separated, `getring`. +# Warning: getgeom is very slow for a Shapefile. +# If you don't need exteriors and holes to be separated, use `getring`. GI.getgeom(::GI.MultiPolygonTrait, geom::AbstractPolygon, i::Integer) = collect(GI.getgeom(geom))[i] function GI.getgeom(::GI.MultiPolygonTrait, geom::AbstractPolygon{T}) where {T} r1 = GI.getring(geom, 1) From 364275a2ae309c344a07ab249e2ad00d7eecbd95 Mon Sep 17 00:00:00 2001 From: rafaqz Date: Thu, 19 May 2022 10:42:08 +0400 Subject: [PATCH 07/16] bugfixes --- src/polygons.jl | 17 +++++++---------- 1 file changed, 7 insertions(+), 10 deletions(-) diff --git a/src/polygons.jl b/src/polygons.jl index 5fa4d1c..8e6d5bd 100644 --- a/src/polygons.jl +++ b/src/polygons.jl @@ -19,16 +19,13 @@ GI.geomtrait(::LinearRing) = GI.LinearRingTrait() GI.ncoord(lr::LinearRing{P}) where {P} = _ncoord(P) GI.ngeom(::GI.LinearRingTrait, lr::LinearRing) = length(lr) -GI.getgeom(::GI.LinearRingTrait, lr::LinearRing{Point}, i::Integer) = lr.xy[i] -GI.getgeom(::GI.LinearRingTrait, lr::LinearRing{PointM}, i::Integer) = - PointM(lr.xy[i], lr.m[i]) -GI.getgeom(::GI.LinearRingTrait, lr::LinearRing{PointZ}, i::Integer) = - PointZ(lr.xy[i], lr.z[i], lr.m[i]) +GI.getgeom(::GI.LinearRingTrait, lr::LinearRing, i::Integer) = lr[i] -Base.parent(lr::LinearRing) = lr.xy -Base.getindex(lr::LinearRing, i) = getindex(parent(lr), i) -Base.size(p::LinearRing) = (length(p),) -Base.length(lr::LinearRing) = length(parent(lr)) +Base.getindex(lr::LinearRing{Point}, i) = lr.xy[i] +Base.getindex(lr::LinearRing{PointM}, i) = PointM(lr.xy[i], lr.m[i]) +Base.getindex(lr::LinearRing{PointZ}, i) = PointZ(lr.xy[i], lr.m[i], lr.z[i]) +Base.size(lr::LinearRing) = (length(lr),) +Base.length(lr::LinearRing) = length(lr.xy) struct SubPolygon{L<:LinearRing} <: AbstractVector{L} rings::Vector{L} @@ -105,7 +102,7 @@ function GI.getgeom(::GI.MultiPolygonTrait, geom::AbstractPolygon{T}) where {T} return polygons end -function _inring(pt::Point, ring::LinearRing) +function _inring(pt::AbstractPoint, ring::LinearRing) intersect(i, j) = (i.y >= pt.y) != (j.y >= pt.y) && (pt.x <= (j.x - i.x) * (pt.y - i.y) / (j.y - i.y) + i.x) From 596fc46258aa3eb7e4248862c33c38eb581d37bd Mon Sep 17 00:00:00 2001 From: rafaqz Date: Thu, 19 May 2022 10:42:20 +0400 Subject: [PATCH 08/16] fix 3d polygon coords --- test/runtests.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/test/runtests.jl b/test/runtests.jl index b364355..5a4b830 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -72,7 +72,7 @@ test_tuples = [ ),( path=joinpath("shapelib_testcases", "test11.shp"), geomtype=PolygonZ, - coordinates=[[[[[1.0,1.0,3.35],[2.0,1.0,4.35],[2.0,2.0,5.35],[1.0,2.0,6.35],[1.0,1.0,7.35]]]],[[[[1.0,4.0,13.35],[2.0,4.0,14.35],[2.0,5.0,15.35],[1.0,5.0,16.35],[1.0,4.0,17.35]]]],[[[[1.0,7.0,23.35],[2.0,7.0,24.35],[2.0,8.0,25.35],[1.0,8.0,26.35],[1.0,7.0,27.35]]]],[[[[0.0,0.0,0.0],[0.0,100.0,1.0],[100.0,100.0,2.0],[100.0,0.0,3.0],[0.0,0.0,4.0]],[[10.0,20.0,5.0],[30.0,20.0,6.0],[30.0,40.0,7.0],[10.0,40.0,8.0],[10.0,20.0,9.0]],[[60.0,20.0,10.0],[90.0,20.0,11.0],[90.0,40.0,12.0],[60.0,40.0,13.0],[60.0,20.0,14.0]]]]], + coordinates=[[[[[1.0,1.0,4.45],[2.0,1.0,5.45],[2.0,2.0,6.45],[1.0,2.0,7.45],[1.0,1.0,8.45]]]],[[[[1.0,4.0,14.45],[2.0,4.0,15.45],[2.0,5.0,16.45],[1.0,5.0,17.45],[1.0,4.0,18.45]]]],[[[[1.0,7.0,24.45],[2.0,7.0,25.45],[2.0,8.0,26.45],[1.0,8.0,27.45],[1.0,7.0,28.45]]]],[[[[0.0,0.0,0.0],[0.0,100.0,2.0],[100.0,100.0,4.0],[100.0,0.0,6.0],[0.0,0.0,8.0]],[[10.0,20.0,10.0],[30.0,20.0,12.0],[30.0,40.0,14.0],[10.0,40.0,16.0],[10.0,20.0,18.0]],[[60.0,20.0,20.0],[90.0,20.0,22.0],[90.0,40.0,24.0],[60.0,40.0,26.0],[60.0,20.0,28.0]]]]], extent=Extent(X=(0.0, 100.0), Y=(0.0, 100.0), Z=(0.0, 27.35)), ),( path=joinpath("shapelib_testcases", "test12.shp"), From 806a3763a7e32ca374cb3b4186ce9fa90583044b Mon Sep 17 00:00:00 2001 From: rafaqz Date: Thu, 19 May 2022 15:58:37 +0400 Subject: [PATCH 09/16] add convert methods --- README.md | 3 +- src/Shapefile.jl | 1 + src/common.jl | 5 +++ src/extent.jl | 18 ++++++-- src/multipatch.jl | 10 ++--- src/multipoints.jl | 29 ++++++------- src/points.jl | 13 ++++++ src/polygons.jl | 38 ++++++++--------- src/polylines.jl | 37 +++++++---------- src/utils.jl | 100 +++++++++++++++++++++++++++++++++++++++++++++ test/runtests.jl | 74 +++++++++++++++++++++++---------- 11 files changed, 235 insertions(+), 93 deletions(-) create mode 100644 src/utils.jl diff --git a/README.md b/README.md index 792ac86..fc224db 100644 --- a/README.md +++ b/README.md @@ -20,7 +20,7 @@ table = Shapefile.Table(path) geoms = Shapefile.shapes(table) # whole columns can be retrieved by their name -table.Descriptio # => Union{String, Missing}["Square with triangle missing", "Smaller triangle", missing] +table.Description # => Union{String, Missing}["Square with triangle missing", "Smaller triangle", missing] # example function that iterates over the rows and gathers shapes that meet specific criteria function selectshapes(table) @@ -50,6 +50,7 @@ julia> GeoInterface.coordinates(Shapefile.shape(first(table))) ``` ## Alternative packages + If you want another lightweight pure Julia package for reading feature files, consider also [GeoJSON.jl](https://github.com/JuliaGeo/GeoJSON.jl). diff --git a/src/Shapefile.jl b/src/Shapefile.jl index 06fea61..10981d6 100644 --- a/src/Shapefile.jl +++ b/src/Shapefile.jl @@ -12,6 +12,7 @@ include("polygons.jl") include("polylines.jl") include("multipoints.jl") include("multipatch.jl") +include("utils.jl") # Handle/Table diff --git a/src/common.jl b/src/common.jl index 667cb23..dc8d98e 100644 --- a/src/common.jl +++ b/src/common.jl @@ -4,6 +4,11 @@ abstract type AbstractShape end GeoInterface.isgeometry(::Type{<:AbstractShape}) = true GeoInterface.ncoord(::GI.AbstractGeometryTrait, ::AbstractShape) = 2 # With specific methods when 3 +Base.convert(T::Type{<:AbstractShape}, geom) = Base.convert(T, GI.geomtrait(geom), geom) +Base.convert(::Type{T}, geom::T) where {T<:AbstractShape} = geom +Base.convert(::Type{<:AbstractShape}, ::Nothing, geom) = + throw(ArgumentError("object to convert is not a GeoInterface compatible geometry")) + """ Rect diff --git a/src/extent.jl b/src/extent.jl index 0afe46c..0552a02 100644 --- a/src/extent.jl +++ b/src/extent.jl @@ -1,14 +1,24 @@ -const Shape3D = Union{PolylineZ,PolygonZ,MultiPointZ,MultiPatch} +const ShapeZ = Union{PolylineZ,PolygonZ,MultiPointZ,MultiPatch} +const ShapeM = Union{ShapeZ,PolylineM,PolygonM,MultiPointM} -# extent +# 3d GI.is3d(::GI.AbstractGeometryTrait, ::AbstractShape) = false -GI.is3d(::GI.AbstractGeometryTrait, ::Shape3D) = true +GI.is3d(::GI.AbstractPointTrait, ::AbstractPoint) = false +GI.is3d(::GI.AbstractGeometryTrait, ::ShapeZ) = true +GI.is3d(::GI.AbstractPointTrait, ::PointZ) = true + +# measured +GI.ismeasured(::GI.AbstractGeometryTrait, ::AbstractShape) = false +GI.ismeasured(::GI.AbstractGeometryTrait, ::ShapeM) = true +GI.ismeasured(::GI.AbstractPointTrait, ::Point) = false +GI.ismeasured(::GI.AbstractPointTrait, ::Union{PointM,PointZ}) = true +# extent function GI.extent(x::AbstractShape) rect = x.MBR return Extents.Extent(X=(rect.left, rect.right), Y=(rect.bottom, rect.top)) end -function GI.extent(x::Union{Shape3D,Handle}) +function GI.extent(x::Union{ShapeZ,Handle}) rect = x.MBR return Extents.Extent(X=(rect.left, rect.right), Y=(rect.bottom, rect.top), Z=(x.zrange.left, x.zrange.right)) end diff --git a/src/multipatch.jl b/src/multipatch.jl index 2795dcb..d125b48 100644 --- a/src/multipatch.jl +++ b/src/multipatch.jl @@ -40,15 +40,11 @@ function Base.read(io::IO, ::Type{MultiPatch}) box = read(io, Rect) numparts = read(io, Int32) numpoints = read(io, Int32) - parts = Vector{Int32}(undef, numparts) - read!(io, parts) + parts = _readparts(io, numparts) parttypes = Vector{Int32}(undef, numparts) read!(io, parttypes) - points = Vector{Point}(undef, numpoints) - read!(io, points) - zrange = read(io, Interval) - zvalues = Vector{Float64}(undef, numpoints) - read!(io, zvalues) + points = _readpoints(io, numpoints) + zrange, zvalues = _readfloats(io, numpoints) # mrange = Vector{Float64}(2) # read!(io, mrange) # measures = Vector{Float64}(numpoints) diff --git a/src/multipoints.jl b/src/multipoints.jl index 32e89df..80d5542 100644 --- a/src/multipoints.jl +++ b/src/multipoints.jl @@ -1,6 +1,11 @@ abstract type AbstractMultiPoint{T} <: AbstractShape end +_pointtype(::Type{<:AbstractMultiPoint{T}}) where T = T + +Base.convert(::Type{T}, ::GI.MultiPointTrait, geom) where T<:AbstractMultiPoint = + T(_convert(_pointtype(T), geom)...) + GI.geomtrait(::AbstractMultiPoint) = GI.MultiPointTrait() GI.ngeom(::GI.MultiPointTrait, geom::AbstractMultiPoint) = length(geom.points) GI.ncoord(::GI.MultiPointTrait, geom::AbstractMultiPoint{T}) where {T} = _ncoord(T) @@ -23,13 +28,13 @@ end function Base.read(io::IO, ::Type{MultiPoint}) box = read(io, Rect) numpoints = read(io, Int32) - points = Vector{Point}(undef, numpoints) - read!(io, points) - MultiPoint(box, points) + points = _readpoints(io, numpoints) + return MultiPoint(box, points) end GI.getgeom(::GI.MultiPointTrait, geom::MultiPoint, i::Integer) = geom.points[i] + """ MultiPointM <: AbstractMultiPoint @@ -54,11 +59,8 @@ end function Base.read(io::IO, ::Type{MultiPointM}) box = read(io, Rect) numpoints = read(io, Int32) - points = Vector{Point}(undef, numpoints) - read!(io, points) - mrange = read(io, Interval) - measures = Vector{Float64}(undef, numpoints) - read!(io, measures) + points = _readpoints(io, numpoints) + mrange, measures = _readm(io, numpoints) MultiPointM(box, points, mrange, measures) end @@ -95,13 +97,8 @@ GI.getgeom(::GI.MultiPointTrait, geom::MultiPointZ, i::Integer) = function Base.read(io::IO, ::Type{MultiPointZ}) box = read(io, Rect) numpoints = read(io, Int32) - points = Vector{Point}(undef, numpoints) - read!(io, points) - zrange = read(io, Interval) - zvalues = Vector{Float64}(undef, numpoints) - read!(io, zvalues) - mrange = read(io, Interval) - measures = Vector{Float64}(undef, numpoints) - read!(io, measures) + points = _readpoints(io, numpoints) + zrange, zvalues = _readz(io, numpoints) + mrange, measures = _readm(io, numpoints) MultiPointZ(box, points, zrange, zvalues, mrange, measures) end diff --git a/src/points.jl b/src/points.jl index fb324ab..a7a6786 100644 --- a/src/points.jl +++ b/src/points.jl @@ -27,6 +27,8 @@ function Base.read(io::IO, ::Type{Point}) Point(x, y) end +Base.convert(::Type{<:Point}, ::GI.PointTrait, geom) = Point(GI.x(geom), GI.y(geom)) + """ PointM <: AbstractPoint @@ -50,6 +52,11 @@ function Base.read(io::IO, ::Type{PointM}) PointM(x, y, m) end +function Base.convert(::Type{<:PointM}, ::GI.PointTrait, geom) + m = GI.ismeasured(geom) ? GI.m(geom) : 0.0 + PointM(GI.x(geom), GI.y(geom), m) +end + GI.m(::GI.PointTrait, point::PointM) = point.m """ @@ -77,6 +84,12 @@ function Base.read(io::IO, ::Type{PointZ}) PointZ(x, y, z, m) end +function Base.convert(::Type{<:PointZ}, ::GI.PointTrait, geom) + z = GI.is3d(geom) ? GI.z(geom) : 0.0 + m = GI.ismeasured(geom) ? GI.m(geom) : 0.0 + PointZ(GI.x(geom), GI.y(geom), z, m) +end + GI.m(::GI.PointTrait, point::PointZ) = point.m GI.z(::GI.PointTrait, point::PointZ) = point.z diff --git a/src/polygons.jl b/src/polygons.jl index 8e6d5bd..ca41cc1 100644 --- a/src/polygons.jl +++ b/src/polygons.jl @@ -23,7 +23,7 @@ GI.getgeom(::GI.LinearRingTrait, lr::LinearRing, i::Integer) = lr[i] Base.getindex(lr::LinearRing{Point}, i) = lr.xy[i] Base.getindex(lr::LinearRing{PointM}, i) = PointM(lr.xy[i], lr.m[i]) -Base.getindex(lr::LinearRing{PointZ}, i) = PointZ(lr.xy[i], lr.m[i], lr.z[i]) +Base.getindex(lr::LinearRing{PointZ}, i) = PointZ(lr.xy[i], lr.z[i], lr.m[i]) Base.size(lr::LinearRing) = (length(lr),) Base.length(lr::LinearRing) = length(lr.xy) @@ -45,6 +45,11 @@ Base.push!(p::SubPolygon, x) = Base.push!(parent(p), x) abstract type AbstractPolygon{T} <: AbstractShape end +_pointtype(::Type{<:AbstractPolygon{T}}) where T = T + +Base.convert(::Type{T}, ::GI.MultiPolygonTrait, geom) where T<:AbstractPolygon = + T(_convertparts(_pointtype(T), geom)...) + # Shapefile polygons are OGC multipolygons GI.geomtrait(geom::AbstractPolygon) = GI.MultiPolygonTrait() GI.nring(::GI.MultiPolygonTrait, geom::AbstractPolygon) = length(geom.parts) @@ -59,6 +64,7 @@ function GI.ngeom(::GI.MultiPolygonTrait, geom::AbstractPolygon) return n end + function GI.getring(t::GI.MultiPolygonTrait, geom::AbstractPolygon) return (GI.getring(t, geom, i) for i in eachindex(geom.parts)) end @@ -146,10 +152,8 @@ function Base.read(io::IO, ::Type{Polygon}) box = read(io, Rect) numparts = read(io, Int32) numpoints = read(io, Int32) - parts = Vector{Int32}(undef, numparts) - read!(io, parts) - points = Vector{Point}(undef, numpoints) - read!(io, points) + parts = _readparts(io, numparts) + points = _readpoints(io, numpoints) Polygon(box, parts, points) end @@ -176,13 +180,9 @@ function Base.read(io::IO, ::Type{PolygonM}) box = read(io, Rect) numparts = read(io, Int32) numpoints = read(io, Int32) - parts = Vector{Int32}(undef, numparts) - read!(io, parts) - points = Vector{Point}(undef, numpoints) - read!(io, points) - mrange = read(io, Interval) - measures = Vector{Float64}(undef, numpoints) - read!(io, measures) + parts = _readparts(io, numparts) + points = _readpoints(io, numpoints) + mrange, measures = _readm(io, numpoints) PolygonM(box, parts, points, mrange, measures) end @@ -212,15 +212,9 @@ function Base.read(io::IO, ::Type{PolygonZ}) box = read(io, Rect) numparts = read(io, Int32) numpoints = read(io, Int32) - parts = Vector{Int32}(undef, numparts) - read!(io, parts) - points = Vector{Point}(undef, numpoints) - read!(io, points) - zrange = read(io, Interval) - zvalues = Vector{Float64}(undef, numpoints) - read!(io, zvalues) - mrange = read(io, Interval) - measures = Vector{Float64}(undef, numpoints) - read!(io, measures) + parts = _readparts(io, numparts) + points = _readpoints(io, numpoints) + zrange, zvalues = _readz(io, numpoints) + mrange, measures = _readm(io, numpoints) PolygonZ(box, parts, points, zrange, zvalues, mrange, measures) end diff --git a/src/polylines.jl b/src/polylines.jl index 493cd97..ced4945 100644 --- a/src/polylines.jl +++ b/src/polylines.jl @@ -19,11 +19,16 @@ Base.size(p::LineString) = (length(p),) Base.length(lr::LineString) = length(parent(lr)) Base.getindex(lr::LineString{Point}, i) = lr.xy[i] Base.getindex(lr::LineString{PointM}, i) = PointM(lr.xy[i], lr.m[i]) -Base.getindex(lr::LineString{PointZ}, i) = PointZ(lr.xy[i], lr.m[i], lr.z[i]) +Base.getindex(lr::LineString{PointZ}, i) = PointZ(lr.xy[i], lr.z[i], lr.m[i]) abstract type AbstractPolyline{T} <: AbstractShape end +_pointtype(::Type{<:AbstractPolyline{T}}) where T = T + +Base.convert(::Type{T}, ::GI.MultiLineStringTrait, geom) where T<:AbstractPolyline = + T(_convertparts(_pointtype(T), geom)...) + GI.geomtrait(::AbstractPolyline) = GI.MultiLineStringTrait() GI.ngeom(::GI.MultiLineStringTrait, geom::AbstractPolyline) = length(geom.parts) GI.ncoord(::GI.MultiLineStringTrait, ::AbstractPolyline{T}) where {T} = _ncoord(T) @@ -59,10 +64,8 @@ function Base.read(io::IO, ::Type{Polyline}) box = read(io, Rect) numparts = read(io, Int32) numpoints = read(io, Int32) - parts = Vector{Int32}(undef, numparts) - read!(io, parts) - points = Vector{Point}(undef, numpoints) - read!(io, points) + parts = _readparts(io, numparts) + points = _readpoints(io, numpoints) Polyline(box, parts, points) end @@ -94,13 +97,9 @@ function Base.read(io::IO, ::Type{PolylineM}) box = read(io, Rect) numparts = read(io, Int32) numpoints = read(io, Int32) - parts = Vector{Int32}(undef, numparts) - read!(io, parts) - points = Vector{Point}(undef, numpoints) - read!(io, points) - mrange = read(io, Interval) - measures = Vector{Float64}(undef, numpoints) - read!(io, measures) + parts = _readparts(io, numparts) + points = _readpoints(io, numpoints) + mrange, measures = _readm(io, numpoints) PolylineM(box, parts, points, mrange, measures) end @@ -133,15 +132,9 @@ function Base.read(io::IO, ::Type{PolylineZ}) box = read(io, Rect) numparts = read(io, Int32) numpoints = read(io, Int32) - parts = Vector{Int32}(undef, numparts) - read!(io, parts) - points = Vector{Point}(undef, numpoints) - read!(io, points) - zrange = read(io, Interval) - zvalues = Vector{Float64}(undef, numpoints) - read!(io, zvalues) - mrange = read(io, Interval) - measures = Vector{Float64}(undef, numpoints) - read!(io, measures) + parts = _readparts(io, numparts) + points = _readpoints(io, numpoints) + zrange, zvalues = _readz(io, numpoints) + mrange, measures = _readm(io, numpoints) PolylineZ(box, parts, points, zrange, zvalues, mrange, measures) end diff --git a/src/utils.jl b/src/utils.jl new file mode 100644 index 0000000..8948acc --- /dev/null +++ b/src/utils.jl @@ -0,0 +1,100 @@ +_partvec(n) = Vector{Int32}(undef, n) +_pointvec(n) = Vector{Point}(undef, n) +_floatvec(n) = Vector{Float64}(undef, n) + +function _readparts(io, n) + points = _partvec(n) + read!(io, points) + return points +end + +function _readpoints(io, n) + points = _pointvec(n) + read!(io, points) + return points +end + +_readm(io, n) = _readfloats(io, n) +_readz(io, n) = _readfloats(io, n) +function _readfloats(io, n) + interval = read(io, Interval) + values = _floatvec(n) + read!(io, values) + return interval, values +end + +# Most objects have similar structure, so we share conversion internals +function _convert(pointmode::Type{<:Point}, geom) + n = GI.npoint(geom) + points = _pointvec(n) + for (i, point) in enumerate(GI.getpoint(geom)) + points[i] = Point(GI.x(point), GI.y(point)) + end + MBR = _getbounds(points) + return MBR, points +end +function _convert(::Type{<:PointM}, geom) + hasm = GI.ismeasured(first(GI.getpoint(geom))) + n = GI.npoint(geom) + points = _pointvec(n) + measures = _floatvec(n) + for (i, point) in enumerate(GI.getpoint(geom)) + points[i] = Point(GI.x(point), GI.y(point)) + measures[i] = hasm ? GI.m(point) : 0.0 + end + @show measures geom.measures + mrange = Interval(extrema(measures)...) + MBR = _getbounds(points) + return MBR, points, mrange, measures +end +function _convert(::Type{<:PointZ}, geom) + hasz = GI.is3d(first(GI.getpoint(geom))) + hasm = GI.ismeasured(first(GI.getpoint(geom))) + n = GI.npoint(geom) + points = _pointvec(n) + zvalues = _floatvec(n) + measures = _floatvec(n) + for (i, point) in enumerate(GI.getpoint(geom)) + points[i] = Point(GI.x(point), GI.y(point)) + zvalues[i] = hasz ? GI.z(point) : 0.0 + measures[i] = hasm ? GI.m(point) : 0.0 + end + mrange = Interval(extrema(measures)...) + zrange = Interval(extrema(zvalues)...) + MBR = _getbounds(points) + return MBR, points, zrange, zvalues, mrange, measures +end + +# _convert with additional parts field +function _convertparts(T, geom) + MBR, args... = _convert(T, geom) + parts = _getparts(geom) + return (MBR, parts, args...) +end + +_getparts(geom) = _getparts(GI.geomtrait(geom), geom) +# Special-case multi polygons because we only store the rings +function _getparts(::GI.MultiPolygonTrait, geom) + parts = _partvec(GI.nring(geom)) + n = 0 + for (i, ring) in enumerate(GI.getring(geom)) + parts[i] = n + n += GI.npoint(ring) + end + return parts +end +function _getparts(::GI.AbstractGeometryTrait, geom) + parts = _partvec(GI.ngeom(geom)) + n = 0 + for (i, g) in enumerate(GI.getgeom(geom)) + parts[i] = n + n += GI.npoint(g) + end + return parts +end + +function _getbounds(points::Vector{Point}) + xrange = extrema(GI.x(p) for p in points) + yrange = extrema(GI.y(p) for p in points) + return Rect(xrange[1], yrange[1], xrange[2], yrange[2]) +end diff --git a/test/runtests.jl b/test/runtests.jl index 5a4b830..37d11f6 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -98,29 +98,61 @@ test_tuples = [ # end # end +points = [Point(1, 3), Point(2, 3), Point(2, 4), Point(1, 4)] +test_shapes = Dict( + Point => Point(1, 2), + PointM => PointM(1, 2, 3), + PointZ => PointZ(1, 2, 3, 4), + MultiPoint => MultiPoint(Rect(1, 3, 2, 4), points), + MultiPointM => MultiPointM(Rect(1, 3, 2, 4), points, Interval(10.0, 40.0), [10.0, 20.0, 30.0, 40.0]), + MultiPointZ => MultiPointZ(Rect(1, 3, 2, 4), points, Interval(1, 4), [1, 2, 3, 4], Interval(10.0, 40.0), [10.0, 20.0, 30.0, 40.0]), + Polygon => Polygon(Rect(1, 3, 2, 4), [0], points), + PolygonM => PolygonM(Rect(1, 3, 2, 4), [0], points, Interval(10.0, 40.0), [10.0, 20.0, 30.0, 40.0]), + PolygonZ => PolygonZ(Rect(1, 3, 2, 4), [0], points, Interval(1, 4), [1, 2, 3, 4], Interval(10.0, 40.0), [10.0, 20.0, 30.0, 40.0]), + Polyline => Polyline(Rect(1, 3, 2, 4), [0], points), + PolylineM => PolylineM(Rect(1, 3, 2, 4), [0], points, Interval(10.0, 40.0), [10.0, 20.0, 30.0, 40.0]), + PolylineZ => PolylineZ(Rect(1, 3, 2, 4), [0], points, Interval(1, 4), [1, 2, 3, 4], Interval(10.0, 40.0), [10.0, 20.0, 30.0, 40.0]), + LineString => LineString{Point}(view(points, 1:4)), + LinearRing => LinearRing{Point}(view(points, 1:4)), + SubPolygon => SubPolygon([LinearRing{Point}(view(points, 1:4))]), +) + @testset "GeoInterface compatability" begin - points = [Point(1, 2), Point(2, 2), Point(2, 1), Point(1, 1)] - - shapes = ( - Point(1, 2), - PointM(1, 2, 3), - PointZ(1, 2, 3, 4), - MultiPoint(Rect(1, 2, 2, 1), points), - MultiPointM(Rect(1, 2, 2, 1), points, Interval(10.0, 40.0), [10.0, 20.0, 30.0, 40.0]), - MultiPointZ(Rect(1, 2, 2, 1), points, Interval(1, 4), [1, 2, 3, 4], Interval(10.0, 40.0), [10.0, 20.0, 30.0, 40.0]), - Polygon(Rect(1, 2, 2, 1), [1], points), - PolygonM(Rect(1, 2, 2, 1), [1], points, Interval(10.0, 40.0), [10.0, 20.0, 30.0, 40.0]), - PolygonZ(Rect(1, 2, 2, 1), [1], points, Interval(1, 4), [1, 2, 3, 4], Interval(10.0, 40.0), [10.0, 20.0, 30.0, 40.0]), - Polyline(Rect(1, 2, 2, 1), [1], points), - PolylineM(Rect(1, 2, 2, 1), [1], points, Interval(10.0, 40.0), [10.0, 20.0, 30.0, 40.0]), - PolylineZ(Rect(1, 2, 2, 1), [1], points, Interval(1, 4), [1, 2, 3, 4], Interval(10.0, 40.0), [10.0, 20.0, 30.0, 40.0]), - LineString{Point}(view(points, 1:4)), - LinearRing{Point}(view(points, 1:4)), - SubPolygon([LinearRing{Point}(view(points, 1:4))]), - ) + @test all(s -> GeoInterface.testgeometry(s), values(test_shapes)) + @test_broken GeoInterface.testgeometry(MultiPatch(Rect(1, 3, 2, 4), [0], [1], points, Interval(1, 4), [1, 2, 3, 4])) +end - @test all(s -> GeoInterface.testgeometry(s), shapes) - @test_broken GeoInterface.testgeometry(MultiPatch(Rect(1, 2, 2, 1), [1], [1], points, Interval(1, 4), [1, 2, 3, 4])) +@testset "convert" begin + @test convert(Point, PointZ(1, 2, 3, 4)) == Point(1, 2) + @test convert(PointM, PointZ(1, 2, 3, 4)) == PointM(1, 2, 4) + @test convert(PointZ, Point(1, 2)) == PointZ(1, 2, 0, 0) + pl = convert(Polyline, test_shapes[PolylineZ]) + @test pl.MBR == test_shapes[Polyline].MBR + @test pl.points == test_shapes[Polyline].points + pl = convert(PolylineM, test_shapes[PolylineZ]) + @test pl.MBR == test_shapes[PolylineM].MBR + @test pl.points == test_shapes[PolylineM].points + @test pl.measures == test_shapes[PolylineM].measures + pg = convert(Polygon, test_shapes[PolygonZ]) + @test pg.MBR == test_shapes[Polygon].MBR + @test pg.points == test_shapes[Polygon].points + pg = convert(PolygonM, test_shapes[PolygonZ]) + @test pg.MBR == test_shapes[PolygonM].MBR + @test pg.points == test_shapes[PolygonM].points + @test pg.measures == test_shapes[PolygonM].measures + mp = convert(MultiPoint, test_shapes[MultiPointZ]) + @test mp.MBR == test_shapes[MultiPoint].MBR + @test mp.points == test_shapes[MultiPoint].points + mp = convert(MultiPointM, test_shapes[MultiPointZ]) + @test mp.MBR == test_shapes[PolylineM].MBR + @test mp.points == test_shapes[PolylineM].points + @test mp.measures == test_shapes[PolylineM].measures + + mp = convert(MultiPointZ, test_shapes[MultiPointM]) + @test mp.MBR == test_shapes[PolylineZ].MBR + @test mp.points == test_shapes[PolylineZ].points + @test mp.measures == test_shapes[PolylineM].measures + @test mp.zvalues == [0.0, 0.0, 0.0, 0.0] # Missing Z filled with zeros: is this the best default? end @testset "Loading Shapefiles" begin From b196a1dd10ae8e33db5d8fde069df514b8a5756b Mon Sep 17 00:00:00 2001 From: rafaqz Date: Thu, 19 May 2022 17:19:04 +0400 Subject: [PATCH 10/16] drop support below 1.6 --- .github/workflows/CI.yml | 2 +- Project.toml | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/.github/workflows/CI.yml b/.github/workflows/CI.yml index a814b18..0c77eba 100644 --- a/.github/workflows/CI.yml +++ b/.github/workflows/CI.yml @@ -13,7 +13,7 @@ jobs: fail-fast: true matrix: version: - - '1.0' + - '1.6' - '1' os: - ubuntu-latest diff --git a/Project.toml b/Project.toml index b5382f7..57e7532 100644 --- a/Project.toml +++ b/Project.toml @@ -20,7 +20,7 @@ GeoInterfaceRecipes = "1.0" Extents = "0.1" RecipesBase = "1" Tables = "0.2, 1" -julia = "1" +julia = "1.6" [extras] DataFrames = "a93c6f00-e57d-5684-b7b6-d8193f3e46c0" From 36e6d8a88ed209ed499f9ed011123b9312cf5275 Mon Sep 17 00:00:00 2001 From: rafaqz Date: Thu, 19 May 2022 17:19:40 +0400 Subject: [PATCH 11/16] revert typo its not a typo --- README.md | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/README.md b/README.md index fc224db..7f15a80 100644 --- a/README.md +++ b/README.md @@ -20,7 +20,7 @@ table = Shapefile.Table(path) geoms = Shapefile.shapes(table) # whole columns can be retrieved by their name -table.Description # => Union{String, Missing}["Square with triangle missing", "Smaller triangle", missing] +table.Descriptio # => Union{String, Missing}["Square with triangle missing", "Smaller triangle", missing] # example function that iterates over the rows and gathers shapes that meet specific criteria function selectshapes(table) From eafd274bc3a1d391c0e865c1a5604b9faf548e30 Mon Sep 17 00:00:00 2001 From: rafaqz Date: Mon, 23 May 2022 17:33:08 +0400 Subject: [PATCH 12/16] bugfix and cleanup --- src/multipatch.jl | 11 ----------- src/utils.jl | 1 - test/runtests.jl | 4 ++-- 3 files changed, 2 insertions(+), 14 deletions(-) diff --git a/src/multipatch.jl b/src/multipatch.jl index d125b48..c64e894 100644 --- a/src/multipatch.jl +++ b/src/multipatch.jl @@ -25,17 +25,6 @@ GI.geomtrait(::MultiPatch) = GI.GeometryCollectionTrait() GI.ngeom(geom::MultiPatch) = length(geom.parts) GeoInterface.ncoord(::MultiPatch) = 3 -# TODO implement the rest of this -function GI.getgeom(::GI.GeometryCollectionTrait, geom::MultiPatch, i::Integer) - error("`getgeom` not implemented for `MultiPatch`: open a github issue at JuliaGeo/Shapefile.jl if you need this.") - # T = SHAPETYPE[geom.parttypes[i]] -end -function GI.subtrait(geom::MultiPatch) - error("`subtype` not implemented for `MultiPatch`: open a github issue at JuliaGeo/Shapefile.jl if you need this.") - # types = map(i -> SHAPETYPE[i], union(geom.parttypes)) -end - - function Base.read(io::IO, ::Type{MultiPatch}) box = read(io, Rect) numparts = read(io, Int32) diff --git a/src/utils.jl b/src/utils.jl index 8948acc..44e1f77 100644 --- a/src/utils.jl +++ b/src/utils.jl @@ -42,7 +42,6 @@ function _convert(::Type{<:PointM}, geom) points[i] = Point(GI.x(point), GI.y(point)) measures[i] = hasm ? GI.m(point) : 0.0 end - @show measures geom.measures mrange = Interval(extrema(measures)...) MBR = _getbounds(points) return MBR, points, mrange, measures diff --git a/test/runtests.jl b/test/runtests.jl index 37d11f6..40ebdc5 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -57,7 +57,7 @@ test_tuples = [ ),( path=joinpath("shapelib_testcases", "test8.shp"), geomtype=PolylineZ, - coordinates=[[[[1.0,1.0,4.45],[2.0,1.0,5.45],[2.0,2.0,6.45],[1.0,2.0,7.45],[1.0,1.0,8.45]]],[[[1.0,4.0,14.45],[2.0,4.0,15.45],[2.0,5.0,16.45],[1.0,5.0,17.45],[1.0,4.0,18.45]]],[[[1.0,7.0,24.45],[2.0,7.0,25.45],[2.0,8.0,26.45],[1.0,8.0,27.45],[1.0,7.0,28.45]]],[[[0.0,0.0,0.0],[0.0,100.0,2.0],[100.0,100.0,4.0],[100.0,0.0,6.0],[0.0,0.0,8.0]],[[10.0,20.0,10.0],[30.0,20.0,12.0],[30.0,40.0,14.0],[10.0,40.0,16.0],[10.0,20.0,18.0]],[[60.0,20.0,20.0],[90.0,20.0,22.0],[90.0,40.0,24.0],[60.0,40.0,26.0],[60.0,20.0,28.0]]]], + coordinates=[[[[1.0,1.0,3.35],[2.0,1.0,4.35],[2.0,2.0,5.35],[1.0,2.0,6.35],[1.0,1.0,7.35]]],[[[1.0,4.0,13.35],[2.0,4.0,14.35],[2.0,5.0,15.35],[1.0,5.0,16.35],[1.0,4.0,17.35]]],[[[1.0,7.0,23.35],[2.0,7.0,24.35],[2.0,8.0,25.35],[1.0,8.0,26.35],[1.0,7.0,27.35]]],[[[0.0,0.0,0.0],[0.0,100.0,1.0],[100.0,100.0,2.0],[100.0,0.0,3.0],[0.0,0.0,4.0]],[[10.0,20.0,5.0],[30.0,20.0,6.0],[30.0,40.0,7.0],[10.0,40.0,8.0],[10.0,20.0,9.0]],[[60.0,20.0,10.0],[90.0,20.0,11.0],[90.0,40.0,12.0],[60.0,40.0,13.0],[60.0,20.0,14.0]]]], extent=Extent(X=(0.0, 100.0), Y=(0.0, 100.0), Z=(0.0, 27.35)), ),( path=joinpath("shapelib_testcases", "test9.shp"), @@ -72,7 +72,7 @@ test_tuples = [ ),( path=joinpath("shapelib_testcases", "test11.shp"), geomtype=PolygonZ, - coordinates=[[[[[1.0,1.0,4.45],[2.0,1.0,5.45],[2.0,2.0,6.45],[1.0,2.0,7.45],[1.0,1.0,8.45]]]],[[[[1.0,4.0,14.45],[2.0,4.0,15.45],[2.0,5.0,16.45],[1.0,5.0,17.45],[1.0,4.0,18.45]]]],[[[[1.0,7.0,24.45],[2.0,7.0,25.45],[2.0,8.0,26.45],[1.0,8.0,27.45],[1.0,7.0,28.45]]]],[[[[0.0,0.0,0.0],[0.0,100.0,2.0],[100.0,100.0,4.0],[100.0,0.0,6.0],[0.0,0.0,8.0]],[[10.0,20.0,10.0],[30.0,20.0,12.0],[30.0,40.0,14.0],[10.0,40.0,16.0],[10.0,20.0,18.0]],[[60.0,20.0,20.0],[90.0,20.0,22.0],[90.0,40.0,24.0],[60.0,40.0,26.0],[60.0,20.0,28.0]]]]], + coordinates=[[[[[1.0,1.0,3.35],[2.0,1.0,4.35],[2.0,2.0,5.35],[1.0,2.0,6.35],[1.0,1.0,7.35]]]],[[[[1.0,4.0,13.35],[2.0,4.0,14.35],[2.0,5.0,15.35],[1.0,5.0,16.35],[1.0,4.0,17.35]]]],[[[[1.0,7.0,23.35],[2.0,7.0,24.35],[2.0,8.0,25.35],[1.0,8.0,26.35],[1.0,7.0,27.35]]]],[[[[0.0,0.0,0.0],[0.0,100.0,1.0],[100.0,100.0,2.0],[100.0,0.0,3.0],[0.0,0.0,4.0]],[[10.0,20.0,5.0],[30.0,20.0,6.0],[30.0,40.0,7.0],[10.0,40.0,8.0],[10.0,20.0,9.0]],[[60.0,20.0,10.0],[90.0,20.0,11.0],[90.0,40.0,12.0],[60.0,40.0,13.0],[60.0,20.0,14.0]]]]], extent=Extent(X=(0.0, 100.0), Y=(0.0, 100.0), Z=(0.0, 27.35)), ),( path=joinpath("shapelib_testcases", "test12.shp"), From 2fff2ab7b3826cbde5ed32599b6aa77209b24b0f Mon Sep 17 00:00:00 2001 From: rafaqz Date: Mon, 23 May 2022 17:55:43 +0400 Subject: [PATCH 13/16] add geointerface methods for Row --- src/table.jl | 51 ++++++++++++++++++++++++--------------------------- 1 file changed, 24 insertions(+), 27 deletions(-) diff --git a/src/table.jl b/src/table.jl index cfcad17..708fffb 100644 --- a/src/table.jl +++ b/src/table.jl @@ -1,5 +1,28 @@ import GeoInterface, DBFTables, Tables +""" + Row + + Row(geometry, record::DBFTables.Row) + +A struct representing a single record in a shapefile. + +Property names accessable by `row.x` are `geometry` for the +geometry object, and the names of the columns in `record`. +""" +struct Row{T} + geometry::T + record::DBFTables.Row +end + + +Base.getproperty(row::Row, name::Symbol) = getproperty(getfield(row, :record), name) +Base.propertynames(row::Row) = propertynames(getfield(row, :record)) + +GeoInterface.isfeature(t::Row) = true +GeoInterface.geometry(t::Row) = getfield(t, :geometry) +GeoInterface.properties(t::Row) = getfield(t, :record) + """ Table @@ -55,21 +78,6 @@ function Table(path::AbstractString) return Shapefile.Table(shp, dbf) end -""" - Row - - Row(geometry, record::DBFTables.Row) - -A struct representing a single record in a shapefile. - -Property names accessable by `row.x` are `geometry` for the -geometry object, and the names of the columns in `record`. -""" -struct Row{T} - geometry::T - record::DBFTables.Row -end - getshp(t::Table) = getfield(t, :shp) getdbf(t::Table) = getfield(t, :dbf) @@ -88,13 +96,11 @@ Iterate over the rows of a Shapefile.Table, yielding a Shapefile.Row for each ro """ function Base.iterate(t::Table, st = 1) st > length(t) && return nothing - geom = @inbounds shapes(t)[st] + geom = shapes(t)[st] record = DBFTables.Row(getdbf(t), st) return Row(geom, record), st + 1 end -Base.getproperty(row::Row, name::Symbol) = getproperty(getfield(row, :record), name) - function Base.getproperty(t::Table, name::Symbol) if name === :geometry shapes(t) @@ -103,8 +109,6 @@ function Base.getproperty(t::Table, name::Symbol) end end -Base.propertynames(row::Row) = propertynames(getfield(row, :record)) - function Base.propertynames(t::Table) names = propertynames(getdbf(t)) pushfirst!(names, :geometry) @@ -131,13 +135,6 @@ end # TODO generalize these with a future GeoInterface/GeoTables # should probably be geometry/geometries but don't want to claim these names yet -""" - shapes(row::Row) - -Get the geometry associated with a `Row` from a shapefile `Table`. -""" -shape(row::Row) = getfield(row, :geometry) - """ shapes(t::Table) From 0ef37b9e9830f9af346bce4dd1e67acea5c088e2 Mon Sep 17 00:00:00 2001 From: rafaqz Date: Mon, 23 May 2022 18:01:50 +0400 Subject: [PATCH 14/16] better info for partially supported Multipatch --- src/multipatch.jl | 9 ++++++++- 1 file changed, 8 insertions(+), 1 deletion(-) diff --git a/src/multipatch.jl b/src/multipatch.jl index c64e894..11f92f0 100644 --- a/src/multipatch.jl +++ b/src/multipatch.jl @@ -18,14 +18,20 @@ struct MultiPatch <: AbstractShape points::Vector{Point} zrange::Interval zvalues::Vector{Float64} - # measures::Vector{Float64} # (optional) + # Optional, not implemented + # measures::Vector{Float64} end GI.geomtrait(::MultiPatch) = GI.GeometryCollectionTrait() GI.ngeom(geom::MultiPatch) = length(geom.parts) GeoInterface.ncoord(::MultiPatch) = 3 +function GI.getgeom(::GI.GeometryCollectionTrait, geom::MultiPatch, i::Integer) + error("`getgeom` not implemented for `MultiPatch`: open a github issue at JuliaGeo/Shapefile.jl if you need this.") +end + function Base.read(io::IO, ::Type{MultiPatch}) + @info "MultiPatch objects are not fully supported in Shapefile.jl" box = read(io, Rect) numparts = read(io, Int32) numpoints = read(io, Int32) @@ -34,6 +40,7 @@ function Base.read(io::IO, ::Type{MultiPatch}) read!(io, parttypes) points = _readpoints(io, numpoints) zrange, zvalues = _readfloats(io, numpoints) + # Optional, not implemented. This may cause `read` bugs if measures are present. # mrange = Vector{Float64}(2) # read!(io, mrange) # measures = Vector{Float64}(numpoints) From 1d39740cbbff255d5aa1bfda6b9d55f2b95321d1 Mon Sep 17 00:00:00 2001 From: rafaqz Date: Mon, 23 May 2022 18:11:43 +0400 Subject: [PATCH 15/16] fix missing shape method --- src/table.jl | 7 +++++++ 1 file changed, 7 insertions(+) diff --git a/src/table.jl b/src/table.jl index 708fffb..a702393 100644 --- a/src/table.jl +++ b/src/table.jl @@ -23,6 +23,13 @@ GeoInterface.isfeature(t::Row) = true GeoInterface.geometry(t::Row) = getfield(t, :geometry) GeoInterface.properties(t::Row) = getfield(t, :record) +""" + shape(row::Row) + +Get the geometry associated with a `Row` from a shapefile `Table`. +""" +shape(row::Row) = getfield(row, :geometry) + """ Table From 40d1c215b5c1d6a7ed2c1121021889368f965604 Mon Sep 17 00:00:00 2001 From: rafaqz Date: Mon, 23 May 2022 20:16:09 +0400 Subject: [PATCH 16/16] update readme example output --- README.md | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/README.md b/README.md index 7f15a80..432a710 100644 --- a/README.md +++ b/README.md @@ -44,9 +44,9 @@ Use `GeoInterface.coordinates` to fully decompose the shape data into parts. ```julia # Example of converting the 1st shape of the file into parts (array of coordinates) julia> GeoInterface.coordinates(Shapefile.shape(first(table))) -2-element Array{Array{Array{Array{Float64,1},1},1},1}: - Array{Array{Float64,1},1}[Array{Float64,1}[[20.0, 20.0], ...]] - Array{Array{Float64,1},1}[Array{Float64,1}[[0.0, 0.0], ...]] +2-element Vector{Vector{Vector{Vector{Float64}}}}: + [[[20.0, 20.0], [20.0, 30.0], [30.0, 30.0], [20.0, 20.0]]] + [[[0.0, 0.0], [100.0, 0.0], [100.0, 100.0], [0.0, 100.0], [0.0, 0.0]]] ``` ## Alternative packages