From 9d24c8d07a7034a6d97a611b5c69122a99562464 Mon Sep 17 00:00:00 2001 From: Skylar A Gering Date: Sat, 23 Dec 2023 17:13:02 -0800 Subject: [PATCH] Sg/area dist (#27) * 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 * Start cleanup of signed_area and signed_distance * Finish area and distance updates * Start area tests * Finish area and distance tests * Remove try file * Simplify docs --- src/GeometryOps.jl | 4 +- src/methods/area.jl | 127 +++++++++++++++ src/methods/distance.jl | 275 +++++++++++++++++++++++++++++++++ src/methods/signed_area.jl | 107 ------------- src/methods/signed_distance.jl | 101 ------------ test/methods/area.jl | 45 ++++++ test/methods/distance.jl | 90 +++++++++++ test/methods/signed_area.jl | 6 - test/runtests.jl | 3 +- 9 files changed, 541 insertions(+), 217 deletions(-) create mode 100644 src/methods/area.jl create mode 100644 src/methods/distance.jl delete mode 100644 src/methods/signed_area.jl delete mode 100644 src/methods/signed_distance.jl create mode 100644 test/methods/area.jl create mode 100644 test/methods/distance.jl delete mode 100644 test/methods/signed_area.jl diff --git a/src/GeometryOps.jl b/src/GeometryOps.jl index 9e19dd553..470e4b0c2 100644 --- a/src/GeometryOps.jl +++ b/src/GeometryOps.jl @@ -20,8 +20,8 @@ include("primitives.jl") include("utils.jl") include("methods/bools.jl") -include("methods/signed_distance.jl") -include("methods/signed_area.jl") +include("methods/distance.jl") +include("methods/area.jl") include("methods/centroid.jl") include("methods/intersects.jl") include("methods/contains.jl") diff --git a/src/methods/area.jl b/src/methods/area.jl new file mode 100644 index 000000000..23446e726 --- /dev/null +++ b/src/methods/area.jl @@ -0,0 +1,127 @@ +# # Area and signed area + +export area, signed_area + +#= +## What is area? What is signed area? + +Area is the amount of space occupied by a two-dimensional figure. It is always +a positive value. Signed area is simply the integral over the exterior path of +a polygon, minus the sum of integrals over its interior holes. It is signed such +that a clockwise path has a positive area, and a counterclockwise path has a +negative area. The area is the absolute value of the signed area. + +To provide an example, consider this rectangle: +```@example rect +using GeometryOps +using GeometryOps.GeometryBasics +using Makie + +rect = Polygon([Point(0,0), Point(0,1), Point(1,1), Point(1,0), Point(0, 0)]) +f, a, p = poly(rect; axis = (; aspect = DataAspect())) +``` +This is clearly a rectangle, etc. But now let's look at how the points look: +```@example rect +lines!(a, rect; color = 1:length(coordinates(rect))+1) +f +``` +The points are ordered in a clockwise fashion, which means that the signed area +is negative. If we reverse the order of the points, we get a postive area. + +## Implementation + +This is the GeoInterface-compatible implementation. First, we implement a +wrapper method that dispatches to the correct implementation based on the +geometry trait. This is also used in the implementation, since it's a lot less +work! + +Note that area (and signed area) are zero for all points and curves, even +if the curves are closed like with a linear ring. Also note that signed area +really only makes sense for polygons, given with a multipolygon can have several +polygons each with a different orientation and thus the absolute value of the +signed area might not be the area. This is why signed area is only implemented +for polygons. +=# + +""" + area(geom)::Real + +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 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. +""" +area(geom) = area(GI.trait(geom), geom) + +""" + signed_area(geom)::Real + +Returns the signed area of the geometry, based on winding order. This is +computed slighly differently for different geometries: + - The signed area of a point is always zero. + - The signed area of a curve is always zero. + - The signed area of a polygon is computed with the shoelace formula and is + positive if the polygon coordinates wind clockwise and negative if + 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) + +# Curves +area(::CT, curve) where CT <: GI.AbstractCurveTrait = + zero(typeof(GI.x(GI.getpoint(curve, 1)))) + +signed_area(trait::CT, curve) where CT <: GI.AbstractCurveTrait = + area(trait, curve) + +# Polygons +area(trait::GI.PolygonTrait, geom) = abs(signed_area(trait, geom)) + +function signed_area(::GI.PolygonTrait, poly) + s_area = _signed_area(GI.getexterior(poly)) + area = abs(s_area) + # Remove hole areas from total + for hole in GI.gethole(poly) + area -= abs(_signed_area(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))) + +#= +Helper function: + +Calculates the signed area of a given curve. This is equivalent to integrating +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 + np = GI.npoint(geom) + 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 +end \ No newline at end of file diff --git a/src/methods/distance.jl b/src/methods/distance.jl new file mode 100644 index 000000000..c14c295f0 --- /dev/null +++ b/src/methods/distance.jl @@ -0,0 +1,275 @@ +# # Distance and signed distance + +export distance, signed_distance + +#= +## What is distance? What is signed distance? + +Distance is the distance of a point to another geometry. This is always a +positive number. If a point is inside of geometry, so on a curve or inside of a +polygon, the distance will be zero. Signed distance is mainly used for polygons +and multipolygons. If a point is outside of a geometry, signed distance has the +same value as distance. However, points within the geometry have a negative +distance representing the distance of a point to the closest boundary. +Therefore, for all "non-filled" geometries, like curves, the distance will +either be postitive or 0. + +To provide an example, consider this rectangle: +```@example rect +using GeometryOps +using GeometryOps.GeometryBasics +using Makie + +rect = Polygon([Point(0,0), Point(0,1), Point(1,1), Point(1,0), Point(0, 0)]) +point_in = Point(0.5, 0.5) +point_out = Point(0.5, 1.5) +f, a, p = poly(rect; axis = (; aspect = DataAspect())) +scatter!(f, point_in) +scatter!(f, point_out) +f +``` +This is clearly a rectangle with one point inside and one point outside. The +points are both an equal distance to the polygon. The distance to point_in is +negative while the distance to point_out is positive. +```@example rect +distance(point_in, poly) # == 0 +signed_distance(point_in, poly) # < 0 +signed_distance(point_out, poly) # > 0 +``` + +## Implementation + +This is the GeoInterface-compatible implementation. First, we implement a +wrapper method that dispatches to the correct implementation based on the +geometry trait. This is also used in the implementation, since it's a lot less +work! + +Distance and signed distance are only implemented for points to other geometries +right now. This could be extended to include distance from other geometries in +the future. + +The distance calculated is the Euclidean distance using the Pythagorean theorem. +Also note that singed_distance only makes sense for "filled-in" shapes, like +polygons, so it isn't implemented for curves. +=# + +""" + distance(point, geom)::Real + +Calculates the ditance from the geometry `g1` to the `point`. The distance +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 + point to the closest segment of the linestring. + - The distance from a point to a linear ring is the minimum distance from + the point to the closest segment of the linear ring. + - 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. +""" +distance(point, geom) = distance( + GI.trait(point), point, + GI.trait(geom), geom, +) + +""" + signed_distance(point, geom)::Real + +Calculates the signed distance from the geometry `geom` to the given point. +Points within `geom` have a negative signed distance, and points outside of +`geom` have a positive signed distance. + - The signed distance from a point to a point, line, linestring, or linear + ring is equal to the distance between the two. + - The signed distance from a point to a polygon is negative if the point is + 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. +""" +signed_distance(point, geom) = signed_distance( + 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) + +# Point-Point +distance(::GI.PointTrait, point, ::GI.PointTrait, geom) = + euclid_distance(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 + +# Point-Line +distance(::GI.PointTrait, point, ::GI.LineTrait, geom) = + _distance_line(point, GI.getpoint(geom, 1), GI.getpoint(geom, 2)) + +# Point-LineString +distance(::GI.PointTrait, point, ::GI.LineStringTrait, geom) = + _distance_curve(point, geom, close_curve = false) + +# Point-LinearRing +distance(::GI.PointTrait, point, ::GI.LinearRingTrait, geom) = + _distance_curve(point, geom, close_curve = true) + +# Point-Polygon +function distance(::GI.PointTrait, point, ::GI.PolygonTrait, geom) + T = typeof(GI.x(point)) + GI.within(point, geom) && return zero(T) + return _distance_polygon(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 +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 +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)) + 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 given their x and y values. +Base.@propagate_inbounds _euclid_distance(x1, y1, x2, y2) = + 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) + x0, y0 = GeoInterface.x(p0), GeoInterface.y(p0) + x1, y1 = GeoInterface.x(p1), GeoInterface.y(p1) + x2, y2 = GeoInterface.x(p2), GeoInterface.y(p2) + + xfirst, yfirst, xlast, ylast = x1 < x2 ? + (x1, y1, x2, y2) : (x2, y2, x1, y1) + + #= + Vectors from first point to last point (v) and from first point to point of + interest (w) to find the projection of w onto v to find closest point + =# + v = (xlast - xfirst, ylast - yfirst) + w = (x0 - xfirst, y0 - yfirst) + + c1 = sum(w .* v) + if c1 <= 0 # p0 is closest to first endpoint + return _euclid_distance(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) + end + + b2 = c1 / c2 # projection fraction + return _euclid_distance(x0, y0, xfirst + (b2 * v[1]), yfirst + (b2 * v[2])) +end + + +#= +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 + 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)) + 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) + min_dist = dist < min_dist ? dist : min_dist + p1 = p2 + end + return min_dist +end + +#= +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) + @inbounds for hole in GI.gethole(poly) + dist = _distance_curve(point, hole; close_curve = true) + min_dist = dist < min_dist ? dist : min_dist + end + return min_dist +end + + diff --git a/src/methods/signed_area.jl b/src/methods/signed_area.jl deleted file mode 100644 index c485d66c4..000000000 --- a/src/methods/signed_area.jl +++ /dev/null @@ -1,107 +0,0 @@ -# # Signed area - -export signed_area - -# ## What is signed area? - -# Signed area is simply the integral over the exterior path of a polygon, -# minus the sum of integrals over its interior holes. - -# It is signed such that a clockwise path has a positive area, and a -# counterclockwise path has a negative area. - -# To provide an example, consider this rectangle: -# ```@example rect -# using GeometryOps -# using GeometryOps.GeometryBasics -# using Makie -# -# rect = Polygon([Point(0,0), Point(0,1), Point(1,1), Point(1,0), Point(0, 0)]) -# f, a, p = poly(rect; axis = (; aspect = DataAspect())) -# ``` -# This is clearly a rectangle, etc. But now let's look at how the points look: -# ```@example rect -# lines!(a, rect; color = 1:length(coordinates(rect))+1) -# f -# ``` -# The points are ordered in a clockwise fashion, which means that the signed area -# is negative. If we reverse the order of the points, we get a postive area. - -# ## Implementation - -# This is the GeoInterface-compatible implementation. - -# First, we implement a wrapper method that dispatches to the correct -# implementation based on the geometry trait. -# -# This is also used in the implementation, since it's a lot less work! -""" - signed_area(geom)::Real - -Returns the signed area of the geometry, based on winding order. -""" -signed_area(x) = signed_area(GI.trait(x), x) - -# TODOS here: -# 1. This could conceivably be multithreaded. How to indicate that it should be so? -# 2. What to do for corner cases (nan point, etc)? -function signed_area(::Union{GI.LineStringTrait,GI.LinearRingTrait}, geom) - # Basically, we integrate the area under the line string, which gives us - # the signed area. - point₁ = GI.getpoint(geom, 1) - point₂ = GI.getpoint(geom, 2) - area = GI.x(point₁) * GI.y(point₂) - GI.y(point₁) * GI.x(point₂) - for point in GI.getpoint(geom) - # Advance the point buffers by 1 point - point₁ = point₂ - point₂ = point - # Accumulate the area into `area` - area += GI.x(point₁) * GI.y(point₂) - GI.y(point₁) * GI.x(point₂) - end - area /= 2 - return area -end - -# This subtracts the -function signed_area(::GI.PolygonTrait, geom) - s_area = signed_area(GI.getexterior(geom)) - area = abs(s_area) - for hole in GI.gethole(geom) - area -= abs(signed_area(hole)) - end - return area * sign(s_area) -end - -signed_area(::GI.MultiPolygonTrait, geom) = sum((signed_area(poly) for poly in GI.getpolygon(geom))) - -# This should _theoretically_ work for anything, but I haven't actually tested yet! - -# Below is the original GeometryBasics implementation: - -# # ```julia -# function signed_area(a::Point{2, T}, b::Point{2, T}, c::Point{2, T}) where T -# return ((b[1] - a[1]) * (c[2] - a[2]) - (c[1] - a[1]) * (b[2] - a[2])) / 2 -# end -# -# function signed_area(points::AbstractVector{<: Point{2, T}}) where {T} -# area = sum((points[i][1] * points[i+1][2] - points[i][2] * points[i+1][1] for i in 1:(length(points)-1))) / 2.0 -# end -# -# function signed_area(ls::GeometryBasics.LineString) -# # coords = GeometryBasics.decompose(Point2f, ls) -# return sum((p1[1] * p2[2] - p1[2] * p2[1] for (p1, p2) in ls)) / 2.0#signed_area(coords) -# end -# -# function signed_area(poly::GeometryBasics.Polygon{2}) -# s_area = signed_area(poly.exterior) -# area = abs(s_area) -# for hole in poly.interiors -# area -= abs(signed_area(hole)) -# end -# return area * sign(s_area) -# end -# -# # WARNING: this may not do what you expect, since it's -# # sensitive to winding order. Use GeoInterface.area instead. -# signed_area(mp::MultiPolygon) = sum(signed_area.(mp.polygons)) -# ``` diff --git a/src/methods/signed_distance.jl b/src/methods/signed_distance.jl deleted file mode 100644 index c0ccb0e20..000000000 --- a/src/methods/signed_distance.jl +++ /dev/null @@ -1,101 +0,0 @@ -# # Signed distance - -export signed_distance - -# TODO: clean this up. It already supports GeoInterface. - -Base.@propagate_inbounds euclid_distance(p1, p2) = sqrt((GeoInterface.x(p2)-GeoInterface.x(p1))^2 + (GeoInterface.y(p2)-GeoInterface.y(p1))^2) -euclid_distance(x1, y1, x2, y2) = sqrt((x2-x1)^2 + (y2-y1)^2) - - - -" Distance from p0 to the line segment formed by p1 and p2. Implementation from Turf.jl." -function _distance(p0, p1, p2) - x0, y0 = GeoInterface.x(p0), GeoInterface.y(p0) - x1, y1 = GeoInterface.x(p1), GeoInterface.y(p1) - x2, y2 = GeoInterface.x(p2), GeoInterface.y(p2) - - if x1 < x2 - xfirst, yfirst = x1, y1 - xlast, ylast = x2, y2 - else - xfirst, yfirst = x2, y2 - xlast, ylast = x1, y1 - end - - v = (xlast - xfirst, ylast - yfirst) - w = (x0 - xfirst, y0 - yfirst) - - c1 = sum(w .* v) - if c1 <= 0 - return euclid_distance(x0, y0, xfirst, yfirst) - end - - c2 = sum(v .* v) - - if c2 <= c1 - return euclid_distance(x0, y0, xlast, ylast) - end - - b2 = c1 / c2 - - return euclid_distance(x0, y0, xfirst + (b2 * v[1]), yfirst + (b2 * v[2])) -end - - -function _distance(linestring, xy) - mindist = typemax(Float64) - N = GeoInterface.npoint(linestring) - @assert N ≥ 3 - p1 = GeoInterface.getpoint(linestring, 1) - p2 = p1 - - for point_ind in 2:N - p2 = GeoInterface.getpoint(linestring, point_ind) - newdist = _distance(xy, p1, p2) - if newdist < mindist - mindist = newdist - end - p1 = p2 - end - - return mindist -end - -function signed_distance(::GeoInterface.PolygonTrait, poly, x, y) - - xy = (x, y) - mindist = _distance(GeoInterface.getexterior(poly), xy) - - @inbounds for hole in GeoInterface.gethole(poly) - newdist = _distance(hole, xy) - if newdist < mindist - mindist = newdist - end - end - - if GeoInterface.contains(poly, GeoInterface.convert(Base.parentmodule(typeof(poly)), (x, y))) - return mindist - else - return -mindist - end -end - -function signed_distance(::GeoInterface.MultiPolygonTrait, multipoly, x, y) - distances = signed_distance.(GeoInterface.getpolygon(multipoly), x, y) - max_val, max_ind = findmax(distances) - return max_val -end - - -""" - signed_distance(geom, x::Real, y::Real)::Float64 - -Calculates the signed distance from the geometry `geom` to the point -defined by `(x, y)`. Points within `geom` have a negative distance, -and points outside of `geom` have a positive distance. - -If `geom` is a MultiPolygon, then this function returns the maximum distance -to any of the polygons in `geom`. -""" -signed_distance(geom, x, y) = signed_distance(GeoInterface.geomtrait(geom), geom, x, y) diff --git a/test/methods/area.jl b/test/methods/area.jl new file mode 100644 index 000000000..b69e2797e --- /dev/null +++ b/test/methods/area.jl @@ -0,0 +1,45 @@ +pt = LG.Point([0.0, 0.0]) +l1 = LG.LineString([[0.0, 0.0], [0.5, 0.5], [1.0, 0.5]]) +r1 = LG.LinearRing([[0.0, 0.0], [1.0, 0.0], [1.0, 2.0], [0.0, 0.0]]) +p1 = LG.Polygon([ + [[10.0, 0.0], [30.0, 0.0], [30.0, 20.0], [10.0, 20.0], [10.0, 0.0]], +]) +p2 = LG.Polygon([ + [[10.0, 0.0], [10.0, 20.0], [30.0, 20.0], [30.0, 0.0], [10.0, 0.0]], + [[15.0, 1.0], [15.0, 11.0], [25.0, 11.0], [25.0, 1.0], [15.0, 1.0]], +]) +p3 = LG.Polygon([ + [[10.0, 0.0], [10.0, 20.0], [30.0, 20.0], [30.0, 0.0], [10.0, 0.0]], + [[15.0, 1.0], [25.0, 1.0], [25.0, 11.0], [15.0, 11.0], [15.0, 1.0]], +]) +p4 = LG.Polygon([ + [ + [0.0, 5.0], [2.0, 2.0], [5.0, 2.0], [2.0, -2.0], [5.0, -5.0], + [0.0, -2.0], [-5.0, -5.0], [-2.0, -2.0], [-5.0, 2.0], [-2.0, 2.0], + [0.0, 5.0], + ], +]) +mp1 = LG.MultiPolygon([p2, p4]) + + +# Points, lines, and rings have zero area +@test GO.area(pt) == GO.signed_area(pt) == LG.area(pt) == 0 +@test GO.area(l1) == GO.signed_area(l1) == LG.area(l1) == 0 +@test GO.area(r1) == GO.signed_area(r1) == LG.area(r1) == 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 +# CW polygons have negative signed area +a2 = LG.area(p2) +@test GO.area(p2) == a2 +@test GO.signed_area(p2) == -a2 +# Winding order of holes doesn't affect sign of signed area +@test GO.signed_area(p3) == -a2 +# Concave polygon correctly calculates area +a4 = LG.area(p4) +@test GO.area(p4) == a4 +@test GO.signed_area(p4) == -a4 +# Multipolygon calculations work +@test GO.area(mp1) == a2 + a4 diff --git a/test/methods/distance.jl b/test/methods/distance.jl new file mode 100644 index 000000000..e4c6499e1 --- /dev/null +++ b/test/methods/distance.jl @@ -0,0 +1,90 @@ +pt1 = LG.Point([0.0, 0.0]) +pt2 = LG.Point([0.0, 1.0]) +pt3 = LG.Point([2.5, 2.5]) +pt4 = LG.Point([3.0, 3.0]) +pt5 = LG.Point([5.1, 5.0]) +pt6 = LG.Point([3.0, 1.0]) +pt7 = LG.Point([0.1, 4.9]) +pt8 = LG.Point([2.0, 1.1]) +pt9 = LG.Point([3.5, 3.1]) +pt10 = LG.Point([10.0, 10.0]) +pt11 = LG.Point([2.5, 7.0]) + +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]]) +r2 = LG.LinearRing([[0.0, 0.0], [0.0, 5.0], [5.0, 5.0], [5.0, 0.0], [0.0, 0.0]]) +r3 = LG.LinearRing([[1.0, 1.0], [3.0, 2.0], [4.0, 1.0], [1.0, 1.0]]) +r4 = LG.LinearRing([[4.0, 3.0], [3.0, 3.0], [4.0, 4.0], [4.0, 3.0]]) +r5 = LG.LinearRing([[0.0, 6.0], [2.5, 8.0], [5.0, 6.0], [0.0, 6.0]]) + +p1 = LG.Polygon(r2, [r3, r4]) +p2 = LG.Polygon(r5) + +mp1 = LG.MultiPolygon([p1, p2]) + +# 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) + +# Point and Line + +#Point on line vertex +@test GO.distance(pt1, l1) == GO.distance(l1, pt1) == LG.distance(pt1, l1) +# Point on line edge +@test GO.distance(pt2, l1) == GO.distance(l1, pt2) == LG.distance(pt2, l1) +# Point equidistant from both segments +@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) + +# Point and Ring + +# Point on linear ring +@test GO.distance(pt1, r1) == LG.distance(pt1, r1) +@test GO.distance(pt3, r1) == LG.distance(pt3, r1) +# Point outside of linear ring +@test GO.distance(pt5, r1) ≈ LG.distance(pt5, r1) +# Point inside of hole created by linear ring +@test GO.distance(pt3, r1) ≈ LG.distance(pt3, r1) +@test GO.distance(pt4, r1) ≈ LG.distance(pt4, r1) + +# Point and Polygon +# Point on polygon exterior edge +@test GO.distance(pt1, p1) == LG.distance(pt1, p1) +@test GO.signed_distance(pt1, p1) == 0 +@test GO.distance(pt2, p1) == LG.distance(pt2, p1) +# Point on polygon hole edge +@test GO.distance(pt4, p1) == LG.distance(pt4, p1) +@test GO.signed_distance(pt4, p1) == 0 +@test GO.distance(pt6, p1) == LG.distance(pt6, p1) +# Point inside of polygon +@test GO.distance(pt3, p1) == LG.distance(pt3, p1) +@test GO.signed_distance(pt3, p1) ≈ + -(min(LG.distance(pt3, r2), LG.distance(pt3, r3), LG.distance(pt3, r4))) +@test GO.distance(pt7, p1) == LG.distance(pt7, p1) +@test GO.signed_distance(pt7, p1) ≈ + -(min(LG.distance(pt7, r2), LG.distance(pt7, r3), LG.distance(pt7, r4))) +# Point outside of polyon exterior +@test GO.distance(pt5, p1) ≈ LG.distance(pt5, p1) +@test GO.signed_distance(pt5, p1) ≈ LG.distance(pt5, p1) +# Point inside of polygon hole +@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) + +# Point and MultiPolygon +# Point outside of either polygon +@test GO.distance(pt5, mp1) ≈ LG.distance(pt5, mp1) +@test GO.distance(pt10, mp1) ≈ LG.distance(pt10, mp1) +# Point within one polygon +@test GO.distance(pt3, mp1) == LG.distance(pt3, mp1) +@test GO.signed_distance(pt3, mp1) ≈ + -(min(LG.distance(pt3, r2), LG.distance(pt3, r3), LG.distance(pt3, r4), LG.distance(pt3, r5))) +@test GO.distance(pt11, mp1) == LG.distance(pt11, mp1) +@test GO.signed_distance(pt11, mp1) ≈ + -(min(LG.distance(pt11, r2), LG.distance(pt11, r3), LG.distance(pt11, r4), LG.distance(pt11, r5))) + diff --git a/test/methods/signed_area.jl b/test/methods/signed_area.jl deleted file mode 100644 index 2abae1930..000000000 --- a/test/methods/signed_area.jl +++ /dev/null @@ -1,6 +0,0 @@ -@testset "ArchGDAL" begin - # Yoinked from ArchGDAL tests - p1 = AG.fromWKT("POLYGON((0 0, 10 0, 10 10, 0 10, 0 0))") - @test GI.area(p1) ≈ abs(GeometryOps.signed_area(p1)) - @test GeometryOps.signed_area(p1) > 0 # test that the signed area is positive -end \ No newline at end of file diff --git a/test/runtests.jl b/test/runtests.jl index ee2065017..d077acf2d 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -20,7 +20,8 @@ const GO = GeometryOps @testset "Centroid" begin include("methods/centroid.jl") end @testset "Equals" begin include("methods/equals.jl") end @testset "Intersect" begin include("methods/intersects.jl") end - @testset "Signed Area" begin include("methods/signed_area.jl") end + @testset "Area" begin include("methods/area.jl") end + @testset "Distance" begin include("methods/distance.jl") end @testset "Overlaps" begin include("methods/overlaps.jl") end # Transformations @testset "Reproject" begin include("transformations/reproject.jl") end