From 4544514db4e00fef94894562119a8969bea388b3 Mon Sep 17 00:00:00 2001 From: Skylar A Gering Date: Thu, 28 Dec 2023 14:25:52 -0800 Subject: [PATCH] Sg/area dist debug (#33) * Fix up intersection point base calculation * Update intersects and add line tests * Add more tests and debug intersects * Add comments to point_in_poly * Remove CairoMakie * Update equals and overlaps * Remove use of findfirst for 1.6 compat * Updated geom, multi-geom equality * Update area for empty geoms * Polygon area base case * Update GI call * testing updates * Update empty check * Added area calculations for empty geoms and collections * Remove try file * Update area and dist with type param --- .gitignore | 1 + src/methods/area.jl | 66 +++++++------ src/methods/distance.jl | 196 +++++++++++++++++++-------------------- test/methods/area.jl | 41 +++++++- test/methods/distance.jl | 21 +++++ 5 files changed, 195 insertions(+), 130 deletions(-) diff --git a/.gitignore b/.gitignore index d6b4b1c05..6977e6f71 100644 --- a/.gitignore +++ b/.gitignore @@ -1,2 +1,3 @@ /Manifest.toml /docs/build/ +.vscode/ \ No newline at end of file diff --git a/src/methods/area.jl b/src/methods/area.jl index 23446e726..5b256dbbd 100644 --- a/src/methods/area.jl +++ b/src/methods/area.jl @@ -44,19 +44,25 @@ for polygons. =# """ - area(geom)::Real + area(geom, ::Type{T} = Float64)::T Returns the area of the geometry. This is computed slighly differently for different geometries: - - The area of a point is always zero. - - The area of a curve is always zero. + - The area of a point/multipoint is always zero. + - The area of a curve/multicurve is always zero. - The area of a polygon is the absolute value of the signed area. - The area multi-polygon is the sum of the areas of all of the sub-polygons. + - The area of a geometry collection is the sum of the areas of all of the + sub-geometries. + +Result will be of type T, where T is an optional argument with a default value +of Float64. """ -area(geom) = area(GI.trait(geom), geom) +area(geom, ::Type{T} = Float64) where T <: AbstractFloat = + _area(T, GI.trait(geom), geom) """ - signed_area(geom)::Real + signed_area(geom, ::Type{T} = Float64)::T Returns the signed area of the geometry, based on winding order. This is computed slighly differently for different geometries: @@ -67,38 +73,42 @@ computed slighly differently for different geometries: counterclockwise. - You cannot compute the signed area of a multipolygon as it doesn't have a meaning as each sub-polygon could have a different winding order. -""" -signed_area(geom) = signed_area(GI.trait(geom), geom) - -# Points -area(::GI.PointTrait, point) = zero(typeof(GI.x(point))) -signed_area(trait::GI.PointTrait, point) = area(trait, point) +Result will be of type T, where T is an optional argument with a default value +of Float64. +""" +signed_area(geom, ::Type{T} = Float64) where T <: AbstractFloat = + _signed_area(T, GI.trait(geom), geom) -# Curves -area(::CT, curve) where CT <: GI.AbstractCurveTrait = - zero(typeof(GI.x(GI.getpoint(curve, 1)))) +# Points, MultiPoints, Curves, MultiCurves +_area(::Type{T}, ::GI.AbstractGeometryTrait, geom) where T = zero(T) -signed_area(trait::CT, curve) where CT <: GI.AbstractCurveTrait = - area(trait, curve) +_signed_area(::Type{T}, ::GI.AbstractGeometryTrait, geom) where T = zero(T) # Polygons -area(trait::GI.PolygonTrait, geom) = abs(signed_area(trait, geom)) +_area(::Type{T}, trait::GI.PolygonTrait, poly) where T = + abs(_signed_area(T, trait, poly)) -function signed_area(::GI.PolygonTrait, poly) - s_area = _signed_area(GI.getexterior(poly)) +function _signed_area(::Type{T}, ::GI.PolygonTrait, poly) where T + GI.isempty(poly) && return zero(T) + s_area = _signed_area(T, GI.getexterior(poly)) area = abs(s_area) + area == 0 && return area # Remove hole areas from total for hole in GI.gethole(poly) - area -= abs(_signed_area(hole)) + area -= abs(_signed_area(T, hole)) end # Winding of exterior ring determines sign return area * sign(s_area) end -# MultiPolygons -area(::GI.MultiPolygonTrait, geom) = - sum((area(poly) for poly in GI.getpolygon(geom))) +# # MultiPolygons and GeometryCollections +_area( + ::Type{T}, + ::Union{GI.MultiPolygonTrait, GI.GeometryCollectionTrait}, + geoms, +) where T = + sum((area(geom, T) for geom in GI.getgeom(geoms)), init = zero(T)) #= Helper function: @@ -108,20 +118,20 @@ to find the area under the curve. Even if curve isn't explicitly closed by repeating the first point at the end of the coordinates, curve is still assumed to be closed. =# -function _signed_area(geom) - # Close curve, even if last point isn't explicitly repeated +function _signed_area(::Type{T}, geom) where T + area = zero(T) + # Close curve, even if last point isn't explicitly repeated np = GI.npoint(geom) + np == 0 && return area first_last_equal = equals(GI.getpoint(geom, 1), GI.getpoint(geom, np)) np -= first_last_equal ? 1 : 0 # Integrate the area under the curve p1 = GI.getpoint(geom, np) - T = typeof(GI.x(p1)) - area = zero(T) for i in 1:np p2 = GI.getpoint(geom, i) # Accumulate the area into `area` area += GI.x(p1) * GI.y(p2) - GI.y(p1) * GI.x(p2) p1 = p2 end - return area / 2 + return T(area / 2) end \ No newline at end of file diff --git a/src/methods/distance.jl b/src/methods/distance.jl index c14c295f0..2c2782fa5 100644 --- a/src/methods/distance.jl +++ b/src/methods/distance.jl @@ -54,7 +54,7 @@ polygons, so it isn't implemented for curves. =# """ - distance(point, geom)::Real + distance(point, geom, ::Type{T} = Float64)::T Calculates the ditance from the geometry `g1` to the `point`. The distance will always be positive or zero. @@ -62,8 +62,6 @@ will always be positive or zero. The method will differ based on the type of the geometry provided: - The distance from a point to a point is just the Euclidean distance between the points. - - The distance from a point to a multipolygon is the shortest distance from - a the given point to any point within the multipoint object. - The distance from a point to a line is the minimum distance from the point to the closest point on the given line. - The distance from a point to a linestring is the minimum distance from the @@ -73,18 +71,17 @@ The method will differ based on the type of the geometry provided: - The distance from a point to a polygon is zero if the point is within the polygon and otherwise is the minimum distance from the point to an edge of the polygon. This includes edges created by holes. - - The distance from a point to a multipolygon is zero if the point is within - the multipolygon and otherwise is the minimum distance from the point to the - closest edge of any of the polygons within the multipolygon. This includes - edges created by holes of the polygons as well. + - The distance from a point to a multigeometry or a geometry collection is + the minimum distance between the point and any of the sub-geometries. + +Result will be of type T, where T is an optional argument with a default value +of Float64. """ -distance(point, geom) = distance( - GI.trait(point), point, - GI.trait(geom), geom, -) +distance(point, geom, ::Type{T} = Float64) where T <: AbstractFloat = + _distance(T, GI.trait(point), point, GI.trait(geom), geom) """ - signed_distance(point, geom)::Real + signed_distance(point, geom, ::Type{T} = Float64)::T Calculates the signed distance from the geometry `geom` to the given point. Points within `geom` have a negative signed distance, and points outside of @@ -95,116 +92,114 @@ Points within `geom` have a negative signed distance, and points outside of within the polygon and is positive otherwise. The value of the distance is the minimum distance from the point to an edge of the polygon. This includes edges created by holes. - - The signed distance from a point to a mulitpolygon is negative if the - point is within one of the polygons that make up the multipolygon and is - positive otherwise. The value of the distance is the minimum distance from - the point to an edge of the multipolygon. This includes edges created by - holes of the polygons as well. + - The signed distance from a point to a multigeometry or a geometry + collection is the minimum signed distance between the point and any of the + sub-geometries. + +Result will be of type T, where T is an optional argument with a default value +of Float64. """ -signed_distance(point, geom) = signed_distance( - GI.trait(point), point, - GI.trait(geom), geom, -) +signed_distance(point, geom, ::Type{T} = Float64) where T<:AbstractFloat = + _signed_distance(T, GI.trait(point), point, GI.trait(geom), geom) -# # Distance # Swap argument order to point as first argument -distance(gtrait::GI.AbstractTrait, geom, ptrait::GI.PointTrait, point) = - distance(ptrait, point, gtrait, geom) +_distance( + ::Type{T}, + gtrait::GI.AbstractTrait, geom, + ptrait::GI.PointTrait, point, +) where T = _distance(T, ptrait, point, gtrait, geom) + +_signed_distance( + ::Type{T}, + gtrait::GI.AbstractTrait, geom, + ptrait::GI.PointTrait, point, +) where T = _signed_distance(T, ptrait, point, gtrait, geom) -# Point-Point -distance(::GI.PointTrait, point, ::GI.PointTrait, geom) = - euclid_distance(point, geom) +# Point-Point, Point-Line, Point-LineString, Point-LinearRing +_distance(::Type{T}, ::GI.PointTrait, point, ::GI.PointTrait, geom) where T = + _euclid_distance(T, point, geom) -# Point-MultiPoint -function distance(::GI.PointTrait, point, ::GI.MultiPointTrait, geom) - T = typeof(GI.x(point)) - min_dist = typemax(T) - for p in GI.getpoint(geom) - dist = euclid_distance(point, p) - min_dist = dist < min_dist ? dist : min_dist - end - return min_dist -end +_distance(::Type{T}, ::GI.PointTrait, point, ::GI.LineTrait, geom) where T = + _distance_line(T, point, GI.getpoint(geom, 1), GI.getpoint(geom, 2)) -# Point-Line -distance(::GI.PointTrait, point, ::GI.LineTrait, geom) = - _distance_line(point, GI.getpoint(geom, 1), GI.getpoint(geom, 2)) +_distance(::Type{T}, ::GI.PointTrait, point, ::GI.LineStringTrait, geom) where T = + _distance_curve(T, point, geom, close_curve = false) -# Point-LineString -distance(::GI.PointTrait, point, ::GI.LineStringTrait, geom) = - _distance_curve(point, geom, close_curve = false) +_distance(::Type{T}, ::GI.PointTrait, point, ::GI.LinearRingTrait, geom) where T = + _distance_curve(T, point, geom, close_curve = true) -# Point-LinearRing -distance(::GI.PointTrait, point, ::GI.LinearRingTrait, geom) = - _distance_curve(point, geom, close_curve = true) +_signed_distance(::Type{T}, ptrait::GI.PointTrait, point, gtrait::GI.AbstractTrait, geom) where T = + _distance(T, ptrait, point, gtrait, geom) # Point-Polygon -function distance(::GI.PointTrait, point, ::GI.PolygonTrait, geom) - T = typeof(GI.x(point)) +function _distance(::Type{T}, ::GI.PointTrait, point, ::GI.PolygonTrait, geom) where T GI.within(point, geom) && return zero(T) - return _distance_polygon(point, geom) + return _distance_polygon(T, point, geom) end -# Point-MultiPolygon -function distance(::GI.PointTrait, point, ::GI.MultiPolygonTrait, geom) - min_dist = distance(point, GI.getpolygon(geom, 1)) - for i in 2:GI.npolygon(geom) - min_dist == 0 && return min_dist # point inside of last polygon checked - dist = distance(point, GI.getpolygon(geom, i)) - min_dist = dist < min_dist ? dist : min_dist - end - return min_dist +function _signed_distance(::Type{T}, ::GI.PointTrait, point, ::GI.PolygonTrait, geom) where T + min_dist = _distance_polygon(T, point, geom) + # negative if point is inside polygon + return GI.within(point, geom) ? -min_dist : min_dist end -# # Signed Distance - -# Swap argument order to point as first argument -signed_distance(gtrait::GI.AbstractTrait, geom, ptrait::GI.PointTrait, point) = - signed_distance(ptrait, point, gtrait, geom) -# Point-Point, Point-Line, Point-LineString, Point-LinearRing -signed_distance(ptrait::GI.PointTrait, point, gtrait::GI.AbstractTrait, geom) = - distance(ptrait, point, gtrait, geom) - -# Point-Polygon -function signed_distance(::GI.PointTrait, point, ::GI.PolygonTrait, geom) - min_dist = _distance_polygon(point, geom) - # should be negative if point is inside polygon - return GI.within(point, geom) ? -min_dist : min_dist +# Point-MultiGeometries / Point-GeometryCollections +function _distance( + ::Type{T}, + ::GI.PointTrait, + point, + ::Union{ + GI.MultiPointTrait, GI.MultiCurveTrait, + GI.MultiPolygonTrait, GI.GeometryCollectionTrait, + }, + geoms, +) where T + min_dist = typemax(T) + for g in GI.getgeom(geoms) + dist = distance(point, g, T) + min_dist = dist < min_dist ? dist : min_dist + end + return min_dist end -# Point-Multipolygon -function signed_distance(::GI.PointTrait, point, ::GI.MultiPolygonTrait, geom) - min_dist = signed_distance(point, GI.getpolygon(geom, 1)) - for i in 2:GI.npolygon(geom) - dist = signed_distance(point, GI.getpolygon(geom, i)) +function _signed_distance( + ::Type{T}, + ::GI.PointTrait, + point, + ::Union{ + GI.MultiPointTrait, GI.MultiCurveTrait, + GI.MultiPolygonTrait, GI.GeometryCollectionTrait, + }, + geoms, +) where T + min_dist = typemax(T) + for g in GI.getgeom(geoms) + dist = signed_distance(point, g, T) min_dist = dist < min_dist ? dist : min_dist end return min_dist end - -""" - euclid_distance(p1::Point, p2::Point)::Real - -Returns the Euclidean distance between two points. -""" -Base.@propagate_inbounds euclid_distance(p1, p2) = _euclid_distance( - GeoInterface.x(p1), GeoInterface.y(p1), - GeoInterface.x(p2), GeoInterface.y(p2), -) +# Returns the Euclidean distance between two points. +Base.@propagate_inbounds _euclid_distance(::Type{T}, p1, p2) where T = + _euclid_distance( + T, + GeoInterface.x(p1), GeoInterface.y(p1), + GeoInterface.x(p2), GeoInterface.y(p2), + ) # Returns the Euclidean distance between two points given their x and y values. -Base.@propagate_inbounds _euclid_distance(x1, y1, x2, y2) = - sqrt((x2 - x1)^2 + (y2 - y1)^2) +Base.@propagate_inbounds _euclid_distance(::Type{T}, x1, y1, x2, y2) where T = + T(sqrt((x2 - x1)^2 + (y2 - y1)^2)) #= Returns the minimum distance from point p0 to the line defined by endpoints p1 and p2. =# -function _distance_line(p0, p1, p2) +function _distance_line(::Type{T}, p0, p1, p2) where T x0, y0 = GeoInterface.x(p0), GeoInterface.y(p0) x1, y1 = GeoInterface.x(p1), GeoInterface.y(p1) x2, y2 = GeoInterface.x(p2), GeoInterface.y(p2) @@ -221,16 +216,16 @@ function _distance_line(p0, p1, p2) c1 = sum(w .* v) if c1 <= 0 # p0 is closest to first endpoint - return _euclid_distance(x0, y0, xfirst, yfirst) + return _euclid_distance(T, x0, y0, xfirst, yfirst) end c2 = sum(v .* v) if c2 <= c1 # p0 is closest to last endpoint - return _euclid_distance(x0, y0, xlast, ylast) + return _euclid_distance(T, x0, y0, xlast, ylast) end b2 = c1 / c2 # projection fraction - return _euclid_distance(x0, y0, xfirst + (b2 * v[1]), yfirst + (b2 * v[2])) + return _euclid_distance(T, x0, y0, xfirst + (b2 * v[1]), yfirst + (b2 * v[2])) end @@ -239,19 +234,18 @@ Returns the minimum distance from the given point to the given curve. If close_curve is true, make sure to include the edge from the first to last point of the curve, even if it isn't explicitly repeated. =# -function _distance_curve(point, curve; close_curve = false) - # See if linear ring has explicitly repeated last point in coordinates +function _distance_curve(::Type{T}, point, curve; close_curve = false) where T + # see if linear ring has explicitly repeated last point in coordinates np = GI.npoint(curve) first_last_equal = equals(GI.getpoint(curve, 1), GI.getpoint(curve, np)) close_curve &= first_last_equal np -= first_last_equal ? 1 : 0 - # Find minimum distance - T = typeof(GI.x(point)) + # find minimum distance min_dist = typemax(T) p1 = GI.getpoint(curve, close_curve ? np : 1) for i in (close_curve ? 1 : 2):np p2 = GI.getpoint(curve, i) - dist = _distance_line(point, p1, p2) + dist = _distance_line(T, point, p1, p2) min_dist = dist < min_dist ? dist : min_dist p1 = p2 end @@ -263,10 +257,10 @@ Returns the minimum distance from the given point to an edge of the given polygon, including from edges created by holes. Assumes polygon isn't filled and treats the exterior and each hole as a linear ring. =# -function _distance_polygon(point, poly) - min_dist = _distance_curve(point, GI.getexterior(poly); close_curve = true) +function _distance_polygon(::Type{T}, point, poly) where T + min_dist = _distance_curve(T, point, GI.getexterior(poly); close_curve = true) @inbounds for hole in GI.gethole(poly) - dist = _distance_curve(point, hole; close_curve = true) + dist = _distance_curve(T, point, hole; close_curve = true) min_dist = dist < min_dist ? dist : min_dist end return min_dist diff --git a/test/methods/area.jl b/test/methods/area.jl index b69e2797e..99227c7cd 100644 --- a/test/methods/area.jl +++ b/test/methods/area.jl @@ -1,6 +1,14 @@ pt = LG.Point([0.0, 0.0]) +empty_pt = LG.readgeom("POINT EMPTY") +mpt = LG.MultiPoint([[0.0, 0.0], [1.0, 0.0]]) +empty_mpt = LG.readgeom("MULTIPOINT EMPTY") l1 = LG.LineString([[0.0, 0.0], [0.5, 0.5], [1.0, 0.5]]) +empty_l = LG.readgeom("LINESTRING EMPTY") +ml1 = LG.MultiLineString([[[0.0, 0.0], [0.5, 0.5], [1.0, 0.5]], [[0.0, 0.0], [0.0, 0.1]]]) +empty_ml = LG.readgeom("MULTILINESTRING EMPTY") +empty_l = LG.readgeom("LINESTRING EMPTY") r1 = LG.LinearRing([[0.0, 0.0], [1.0, 0.0], [1.0, 2.0], [0.0, 0.0]]) +empty_r = LG.readgeom("LINEARRING EMPTY") p1 = LG.Polygon([ [[10.0, 0.0], [30.0, 0.0], [30.0, 20.0], [10.0, 20.0], [10.0, 0.0]], ]) @@ -19,18 +27,35 @@ p4 = LG.Polygon([ [0.0, 5.0], ], ]) +empty_p = LG.readgeom("POLYGON EMPTY") mp1 = LG.MultiPolygon([p2, p4]) - +empty_mp = LG.readgeom("MULTIPOLYGON EMPTY") +c = LG.GeometryCollection([p1, p2, r1, empty_l]) +empty_c = LG.readgeom("GEOMETRYCOLLECTION EMPTY") # Points, lines, and rings have zero area @test GO.area(pt) == GO.signed_area(pt) == LG.area(pt) == 0 +@test GO.area(empty_pt) == LG.area(empty_pt) == 0 +@test GO.area(pt) isa Float64 +@test GO.signed_area(pt, Float32) isa Float32 +@test GO.signed_area(pt) isa Float64 +@test GO.area(pt, Float32) isa Float32 +@test GO.area(mpt) == GO.signed_area(mpt) == LG.area(mpt) == 0 +@test GO.area(empty_mpt) == LG.area(empty_mpt) == 0 @test GO.area(l1) == GO.signed_area(l1) == LG.area(l1) == 0 +@test GO.area(empty_l) == LG.area(empty_l) == 0 +@test GO.area(ml1) == GO.signed_area(ml1) == LG.area(ml1) == 0 +@test GO.area(empty_ml) == LG.area(empty_ml) == 0 @test GO.area(r1) == GO.signed_area(r1) == LG.area(r1) == 0 +@test GO.area(empty_r) == LG.area(empty_r) == 0 # Polygons have non-zero area # CCW polygons have positive signed area @test GO.area(p1) == GO.signed_area(p1) == LG.area(p1) @test GO.signed_area(p1) > 0 +# Float32 calculations +@test GO.area(p1) isa Float64 +@test GO.area(p1, Float32) isa Float32 # CW polygons have negative signed area a2 = LG.area(p2) @test GO.area(p2) == a2 @@ -41,5 +66,19 @@ a2 = LG.area(p2) a4 = LG.area(p4) @test GO.area(p4) == a4 @test GO.signed_area(p4) == -a4 +# Empty polygon +@test GO.area(empty_p) == LG.area(empty_p) == 0 +@test GO.signed_area(empty_p) == 0 + # Multipolygon calculations work @test GO.area(mp1) == a2 + a4 +@test GO.area(mp1, Float32) isa Float32 +# Empty multipolygon +@test GO.area(empty_mp) == LG.area(empty_mp) == 0 + + +# Geometry collection summed area +@test GO.area(c) == LG.area(c) +@test GO.area(c, Float32) isa Float32 +# Empty collection +@test GO.area(empty_c) == LG.area(empty_c) == 0 diff --git a/test/methods/distance.jl b/test/methods/distance.jl index e4c6499e1..7d66b13f8 100644 --- a/test/methods/distance.jl +++ b/test/methods/distance.jl @@ -10,6 +10,8 @@ pt9 = LG.Point([3.5, 3.1]) pt10 = LG.Point([10.0, 10.0]) pt11 = LG.Point([2.5, 7.0]) +mpt1 = LG.MultiPoint([pt1, pt2, pt3]) + l1 = LG.LineString([[0.0, 0.0], [0.0, 5.0], [5.0, 5.0]]) r1 = LG.LinearRing([[0.0, 0.0], [0.0, 5.0], [5.0, 5.0], [0.0, 0.0]]) @@ -23,12 +25,17 @@ p2 = LG.Polygon(r5) mp1 = LG.MultiPolygon([p1, p2]) +c1 = LG.GeometryCollection([pt1, r1, p1]) + # Point and Point # Distance from point to same point @test GO.distance(pt1, pt1) == LG.distance(pt1, pt1) # Distance from point to different point @test GO.distance(pt1, pt2) ≈ GO.distance(pt2, pt1) ≈ LG.distance(pt1, pt2) +# Return types +@test GO.distance(pt1, pt1) isa Float64 +@test GO.distance(pt1, pt1, Float32) isa Float32 # Point and Line @@ -40,6 +47,9 @@ mp1 = LG.MultiPolygon([p1, p2]) @test GO.distance(pt3, l1) ≈ GO.distance(l1, pt3) ≈ LG.distance(pt3, l1) # Point closer to one segment than another @test GO.distance(pt4, l1) ≈ GO.distance(l1, pt4) ≈ LG.distance(pt4, l1) +# Return types +@test GO.distance(pt1, l1) isa Float64 +@test GO.distance(pt1, l1, Float32) isa Float32 # Point and Ring @@ -75,6 +85,14 @@ mp1 = LG.MultiPolygon([p1, p2]) @test GO.distance(pt8, p1) ≈ LG.distance(pt8, p1) @test GO.signed_distance(pt8, p1) ≈ LG.distance(pt8, p1) @test GO.distance(pt9, p1) ≈ LG.distance(pt9, p1) +# Return types +@test GO.distance(pt1, p1) isa Float64 +@test GO.distance(pt1, p1, Float32) isa Float32 + +# Point and MultiPoint +@test GO.distance(pt4, mpt1) == LG.distance(pt4, mpt1) +@test GO.distance(pt4, mpt1) isa Float64 +@test GO.distance(pt4, mpt1, Float32) isa Float32 # Point and MultiPolygon # Point outside of either polygon @@ -88,3 +106,6 @@ mp1 = LG.MultiPolygon([p1, p2]) @test GO.signed_distance(pt11, mp1) ≈ -(min(LG.distance(pt11, r2), LG.distance(pt11, r3), LG.distance(pt11, r4), LG.distance(pt11, r5))) +# Point and Geometry Collection +@test GO.distance(pt1, c1) == LG.distance(pt1, c1) +