Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

fix performance and correctness of base points, and test #89

Merged
merged 3 commits into from
Jan 18, 2023
Merged

fix performance and correctness of base points, and test #89

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
3 commits
Select commit Hold shift + click to select a range
20b90ce
fix base performance and correctness, and test
rafaqz Jan 17, 2023
6abf141
more clearly document x/y/z/m for points
rafaqz Jan 17, 2023
a40dccd
test ismeasured and is3d for Tuple
rafaqz Jan 17, 2023
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
79 changes: 62 additions & 17 deletions src/base.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -1,47 +1,89 @@
# Implementation of GeoInterface for Base Types


const PointTuple2 = Tuple{<:Real,<:Real}
const PointTuple3 = Tuple{<:Real,<:Real,<:Real}
const PointTuple4 = Tuple{<:Real,<:Real,<:Real,<:Real}
const PointTuple = Union{PointTuple2,PointTuple3,PointTuple4}
# AbstractVector{<:Real} length 2 - 4 (where length checking is possible and efficient)

GeoInterface.isgeometry(::Type{<:AbstractVector{<:Real}}) = true
GeoInterface.geomtrait(::AbstractVector{<:Real}) = PointTrait()
GeoInterface.ncoord(::PointTrait, geom::AbstractVector{<:Real}) = Base.length(geom)
GeoInterface.ncoord(::PointTrait, geom::AbstractVector{<:Real}) = Base.length(geom) # TODO should this error for length > 4 ?
GeoInterface.getcoord(::PointTrait, geom::AbstractVector{<:Real}, i) = getindex(geom, i)
GeoInterface.is3d(::PointTrait, geom::AbstractVector{<:Real}) = Base.length(geom) in (3, 4)
GeoInterface.ismeasured(::PointTrait, geom::AbstractVector{<:Real}) = Base.length(geom) == 4 # Measured must have Z
rafaqz marked this conversation as resolved.
Show resolved Hide resolved

# x/y/z/m methods were 100x slower without these custom definitions
# Also allow @inbounds to make them slightly faster when you know the sizes
Base.@propagate_inbounds function GeoInterface.x(::PointTrait, geom::AbstractVector{<:Real})
@boundscheck Base.length(geom) in (2, 3, 4) || _xy_error(Base.length(geom))
@inbounds geom[begin]
end
Base.@propagate_inbounds function GeoInterface.y(::PointTrait, geom::AbstractVector{<:Real})
@boundscheck Base.length(geom) in (2, 3, 4) || _xy_error(Base.length(geom))
@inbounds geom[begin + 1]
end
Base.@propagate_inbounds function GeoInterface.z(::PointTrait, geom::AbstractVector{<:Real})
@boundscheck Base.length(geom) in (3, 4) || _z_error(Base.length(geom))
@inbounds geom[begin + 2]
end
Base.@propagate_inbounds function GeoInterface.m(::PointTrait, geom::AbstractVector{<:Real})
@boundscheck Base.length(geom) == 4 || _m_error(Base.length(geom))
@inbounds geom[end]
end

@noinline _xy_error(l) = throw(ArgumentError("Length of point must be 2, 3 or 4 to use `GeoInterface.x(point)` or `GeoInterface.y(point)`, got $l"))
@noinline _z_error(l) = throw(ArgumentError("Length of point must be 3 or 4 to use `GeoInterface.z(point)`, got $l"))
@noinline _m_error(l) = throw(ArgumentError("Length of point must be 4 to use `GeoInterface.m(point)`, got $l"))


# Tuple length 2 - 4

const PointTuple2 = Tuple{<:Real,<:Real}
const PointTuple3 = Tuple{<:Real,<:Real,<:Real}
const PointTuple4 = Tuple{<:Real,<:Real,<:Real,<:Real}
const PointTuple = Union{PointTuple2,PointTuple3,PointTuple4}

GeoInterface.isgeometry(::Type{<:PointTuple}) = true
GeoInterface.geomtrait(::PointTuple) = PointTrait()
GeoInterface.ncoord(::PointTrait, geom::PointTuple) = Base.length(geom)
GeoInterface.getcoord(::PointTrait, geom::PointTuple, i) = getindex(geom, i)
GeoInterface.is3d(::PointTrait, geom::PointTuple2) = false
GeoInterface.is3d(::PointTrait, geom::Union{PointTuple3,PointTuple4}) = true
GeoInterface.ismeasured(::PointTrait, geom::Union{PointTuple2,PointTuple3}) = false
GeoInterface.ismeasured(::PointTrait, geom::PointTuple4) = true

GeoInterface.x(::PointTrait, geom::PointTuple) = geom[1]
GeoInterface.y(::PointTrait, geom::PointTuple) = geom[2]
GeoInterface.z(::PointTrait, geom::PointTuple2) = _z_error(Base.length(geom))
GeoInterface.z(::PointTrait, geom::Union{PointTuple3,PointTuple4}) = geom[3]
GeoInterface.m(::PointTrait, geom::PointTuple4) = geom[4]
GeoInterface.m(::PointTrait, geom::Union{PointTuple2,PointTuple3}) = _m_error(Base.length(geom))

for (i, pointtype) in enumerate((PointTuple2, PointTuple3, PointTuple4))
keys = default_coord_names[1:i+1]
sig = NamedTuple{keys,<:pointtype}
@eval GeoInterface.isgeometry(::Type{<:$sig}) = true
@eval GeoInterface.geomtrait(::$sig) = PointTrait()
@eval GeoInterface.ncoord(::PointTrait, geom::$sig) = $i + 1
@eval GeoInterface.getcoord(::PointTrait, geom::$sig, i) = getindex(geom, i)
end

# NamedTuple
# with X/Y/Z/M names in any order

# Define all possible NamedTuple points
const NamedTuplePoint = Union{
const NamedTuplePointXY = Union{
NamedTuple{(:X, :Y),<:PointTuple2},
NamedTuple{(:Y, :X),<:PointTuple2},
}

const NamedTuplePointZ = Union{
NamedTuple{(:X, :Y, :Z),<:PointTuple3},
NamedTuple{(:X, :Z, :Y),<:PointTuple3},
NamedTuple{(:Z, :Y, :X),<:PointTuple3},
NamedTuple{(:Z, :X, :Y),<:PointTuple3},
NamedTuple{(:Y, :X, :Z),<:PointTuple3},
NamedTuple{(:Y, :Z, :X),<:PointTuple3},
}
const NamedTuplePointM = Union{
NamedTuple{(:X, :Y, :M),<:PointTuple3},
NamedTuple{(:X, :M, :Y),<:PointTuple3},
NamedTuple{(:M, :Y, :X),<:PointTuple3},
NamedTuple{(:M, :X, :Y),<:PointTuple3},
NamedTuple{(:Y, :X, :Z),<:PointTuple3},
NamedTuple{(:Y, :Z, :X),<:PointTuple3},
}
const NamedTuplePointZM = Union{
NamedTuple{(:X, :Y, :Z, :M),<:PointTuple4},
NamedTuple{(:X, :Y, :M, :Z),<:PointTuple4},
NamedTuple{(:X, :Z, :Y, :M),<:PointTuple4},
Expand All @@ -68,16 +110,19 @@ const NamedTuplePoint = Union{
NamedTuple{(:M, :X, :Y, :Z),<:PointTuple4},
}

const NamedTuplePoint = Union{NamedTuplePointXY,NamedTuplePointZ,NamedTuplePointM,NamedTuplePointZM}

_keys(::Type{<:NamedTuple{K}}) where K = K
GeoInterface.isgeometry(::Type{T}) where {T<:NamedTuplePoint} = all(in(default_coord_names), _keys(T))
GeoInterface.isgeometry(::Type{T}) where {T<:NamedTuplePoint} = true
GeoInterface.geomtrait(::NamedTuplePoint) = PointTrait()
GeoInterface.ncoord(::PointTrait, geom::NamedTuplePoint) = Base.length(geom)
GeoInterface.getcoord(::PointTrait, geom::NamedTuplePoint, i) = getindex(geom, i)
GeoInterface.coordnames(::PointTrait, geom::NamedTuplePoint) = _keys(typeof(geom))
GeoInterface.x(::PointTrait, geom::NamedTuplePoint) = geom.X
GeoInterface.y(::PointTrait, geom::NamedTuplePoint) = geom.Y
GeoInterface.z(::PointTrait, geom::NamedTuplePoint) = geom.Z
GeoInterface.m(::PointTrait, geom::NamedTuplePoint) = geom.M
GeoInterface.z(::PointTrait, geom::Union{NamedTuplePointZ,NamedTuplePointZM}) = geom.Z
GeoInterface.z(::PointTrait, geom::NamedTuplePointXY) = throw(ArgumentError("NamedTuple point has no Z field"))
GeoInterface.m(::PointTrait, geom::Union{NamedTuplePointXY,NamedTuplePointZ}) = throw(ArgumentError("NamedTuple point has no M field"))


# Default features using NamedTuple and AbstractArray
Expand Down
14 changes: 13 additions & 1 deletion src/interface.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -561,6 +561,8 @@ convexhull(geom) = convexhull(geomtrait(geom), geom)

Return the :X coordinate of the given `geom`.
Note that this is only valid for [`AbstractPointTrait`](@ref)s.

For `Tuple` and `Vector` points, the first value is returned.
"""
x(geom) = x(geomtrait(geom), geom)

Expand All @@ -569,6 +571,8 @@ x(geom) = x(geomtrait(geom), geom)

Return the :Y coordinate of the given `geom`.
Note that this is only valid for [`AbstractPointTrait`](@ref)s.

For `Tuple` and `Vector` points, the second value is returned.
"""
y(geom) = y(geomtrait(geom), geom)

Expand All @@ -577,14 +581,22 @@ y(geom) = y(geomtrait(geom), geom)

Return the :Z coordinate of the given `geom`.
Note that this is only valid for [`AbstractPointTrait`](@ref)s.

For length 3 `Tuple` and `Vector` points, the third value is returned.
"""
z(geom) = z(geomtrait(geom), geom)

"""
m(geom) -> Number

Return the :M coordinate of the given `geom`.
Return the :M (measured) coordinate of the given `geom`.
Note that this is only valid for [`AbstractPointTrait`](@ref)s.

For length 4 `Tuple` and `Vector` points, the fouth value
is returned.

Length 3 `Tuple` and `Vector` points can *not* represent measured points,
and will throw an `ArgumentError`.
"""
m(geom) = m(geomtrait(geom), geom)

Expand Down
102 changes: 94 additions & 8 deletions test/test_primitives.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -308,21 +308,90 @@ end

@testset "Vector" begin
geom = [1, 2]
@test !GeoInterface.is3d(geom)
@test !GeoInterface.ismeasured(geom)
@test testgeometry(geom)
@test GeoInterface.x(geom) == 1
@test collect(GeoInterface.getcoord(geom)) == geom
@test GeoInterface.ncoord(geom) == 2
@test GeoInterface.x(geom) == 1
@test GeoInterface.y(geom) == 2
@test_throws ArgumentError GeoInterface.z(geom)
@test_throws ArgumentError GeoInterface.m(geom)
geom = [1, 2, 3]
@test testgeometry(geom)
@test collect(GeoInterface.getcoord(geom)) == geom
@test GeoInterface.is3d(geom)
@test !GeoInterface.ismeasured(geom)
@test GeoInterface.x(geom) == 1
@test GeoInterface.y(geom) == 2
@test GeoInterface.z(geom) == 3
@test_throws ArgumentError GeoInterface.m(geom)
geom = [1, 2, 3, 4]
@test testgeometry(geom)
@test collect(GeoInterface.getcoord(geom)) == geom
@test GeoInterface.is3d(geom)
@test GeoInterface.ismeasured(geom)
@test GeoInterface.ncoord(geom) == 4
@test GeoInterface.x(geom) == 1
@test GeoInterface.y(geom) == 2
@test GeoInterface.z(geom) == 3
@test GeoInterface.m(geom) == 4
geom = [1, 2, 3, 4, 5]
@test !GeoInterface.is3d(geom)
@test !GeoInterface.ismeasured(geom)
@test_throws ArgumentError GeoInterface.x(geom)
@test_throws ArgumentError GeoInterface.y(geom)
@test_throws ArgumentError GeoInterface.z(geom)
@test_throws ArgumentError GeoInterface.m(geom)
end

@testset "Tuple" begin
geom = (1, 2.0f0)
@test GeoInterface.trait(geom) isa PointTrait
@test GeoInterface.geomtrait(geom) isa PointTrait
@test testgeometry(geom)
@test GeoInterface.x(geom) == 1
@test GeoInterface.y(geom) == 2.0f0
@test !GeoInterface.is3d(geom)
@test !GeoInterface.ismeasured(geom)
@test GeoInterface.x(geom) === 1
@test GeoInterface.y(geom) === 2.0f0
@test_throws ArgumentError GeoInterface.z(geom)
@test_throws ArgumentError GeoInterface.m(geom)
@test GeoInterface.ncoord(geom) == 2
@test collect(GeoInterface.getcoord(geom)) == [1, 2]
geom = (1, 2, 3.0)
@test GeoInterface.trait(geom) isa PointTrait
@test GeoInterface.geomtrait(geom) isa PointTrait
@test testgeometry(geom)
@test GeoInterface.is3d(geom)
@test !GeoInterface.ismeasured(geom)
@test GeoInterface.x(geom) === 1
@test GeoInterface.y(geom) === 2
@test GeoInterface.z(geom) === 3.0
@test_throws ArgumentError GeoInterface.m(geom)
@test GeoInterface.ncoord(geom) == 3
@test collect(GeoInterface.getcoord(geom)) == [1, 2, 3]
geom = (1, 2, 3, 4.0)
@test GeoInterface.trait(geom) isa PointTrait
@test GeoInterface.geomtrait(geom) isa PointTrait
@test testgeometry(geom)
@test GeoInterface.is3d(geom)
@test GeoInterface.ismeasured(geom)
@test GeoInterface.x(geom) === 1
@test GeoInterface.y(geom) === 2
@test GeoInterface.z(geom) === 3
@test GeoInterface.m(geom) === 4.0
@test GeoInterface.ncoord(geom) == 4
@test collect(GeoInterface.getcoord(geom)) == [1, 2, 3, 4]
geom = (1, 2, 3, 4.0, 5)
@test GeoInterface.isgeometry(geom) == false
@test GeoInterface.trait(geom) isa Nothing
@test GeoInterface.geomtrait(geom) isa Nothing
@test_throws MethodError GeoInterface.x(geom)
@test_throws MethodError GeoInterface.y(geom)
@test_throws MethodError GeoInterface.z(geom)
@test_throws MethodError GeoInterface.m(geom)
@test_throws MethodError GeoInterface.ncoord(geom) == 4
@test_throws MethodError GeoInterface.getcoord(geom)
end

@testset "NamedTuple" begin
Expand All @@ -334,8 +403,20 @@ end
@test testgeometry(geom)
@test GeoInterface.x(geom) == 1
@test GeoInterface.y(geom) == 2.0
@test_throws ArgumentError GeoInterface.z(geom)
@test_throws ArgumentError GeoInterface.m(geom)
@test collect(GeoInterface.getcoord(geom)) == [1, 2]

geom = (; X=1.0, Y=2.0, Z=0x03)
@test testgeometry(geom)
@test GeoInterface.is3d(geom)
@test GeoInterface.ismeasured(geom) == false
@test GeoInterface.coordnames(geom) == (:X, :Y, :Z)
@test GeoInterface.x(geom) == 1
@test GeoInterface.y(geom) == 2.0
@test GeoInterface.z(geom) == 0x03
@test_throws ArgumentError GeoInterface.m(geom)

geom = (; X=1, Y=2, Z=3.0f0, M=4.0)
@test GeoInterface.trait(geom) isa PointTrait
@test GeoInterface.geomtrait(geom) isa PointTrait
Expand All @@ -355,15 +436,20 @@ end
@test GeoInterface.getcoord(geom, 4) === 4.0
@test collect(GeoInterface.getcoord(geom)) == [1, 2, 3, 4]

geom = (; X=1.0, Y=2.0, Z=0x03)
@test testgeometry(geom)
@test GeoInterface.is3d(geom)
@test GeoInterface.ismeasured(geom) == false
@test GeoInterface.coordnames(geom) == (:X, :Y, :Z)
geom = (; Z=3, X=1, Y=2, M=4)
@test testgeometry(geom)
@test collect(GeoInterface.getcoord(geom)) == [3, 1, 2, 4]
@test GeoInterface.coordnames(geom) == (:Z, :X, :Y, :M)

geom = (; A=0, X=1, Y=2, Z=3.0f0, M=4.0)
@test GeoInterface.trait(geom) isa Nothing
@test GeoInterface.geomtrait(geom) isa Nothing
@test_throws MethodError GeoInterface.x(geom)
@test_throws MethodError GeoInterface.y(geom)
@test_throws MethodError GeoInterface.z(geom)
@test_throws MethodError GeoInterface.m(geom)
@test_throws MethodError GeoInterface.ncoord(geom) == 4
@test_throws MethodError GeoInterface.getcoord(geom)
end

@testset "NamedTupleFeature" begin
Expand Down