diff --git a/.github/workflows/CI.yml b/.github/workflows/CI.yml index 87ee77d79..84b450519 100644 --- a/.github/workflows/CI.yml +++ b/.github/workflows/CI.yml @@ -35,6 +35,15 @@ jobs: - uses: julia-actions/cache@v1 - uses: julia-actions/julia-buildpkg@v1 - uses: julia-actions/julia-runtest@v1 + - uses: julia-actions/julia-processcoverage@v1 + with: + directories: src,ext + - name: Upload coverage reports to Codecov + uses: codecov/codecov-action@v4.0.1 + with: + files: lcov.info + token: ${{ secrets.CODECOV_TOKEN }} + slug: JuliaGeo/GeometryOps.jl docs: name: Documentation runs-on: ubuntu-latest @@ -55,6 +64,7 @@ jobs: install-package: false env: GITHUB_TOKEN: ${{ secrets.GITHUB_TOKEN }} + DOCUMENTER_KEY: ${{ secrets.DOCUMENTER_KEY }} - run: | julia --project=docs -e ' using Documenter: DocMeta, doctest diff --git a/Project.toml b/Project.toml index 5d22ae65f..3e21575c6 100644 --- a/Project.toml +++ b/Project.toml @@ -1,10 +1,11 @@ name = "GeometryOps" uuid = "3251bfac-6a57-4b6d-aa61-ac1fef2975ab" authors = ["Anshul Singhvi and contributors"] -version = "0.1.10" +version = "0.1.11" [deps] CoordinateTransformations = "150eb455-5306-5404-9cee-2592286d6298" +DelaunayTriangulation = "927a84f5-c5f4-47a5-9785-b46e178433df" ExactPredicates = "429591f6-91af-11e9-00e2-59fbe8cec110" GeoInterface = "cf35fbd7-0cd7-5166-be24-54bfbe79505f" GeometryBasics = "5c1252a2-5f33-56bf-86c9-59e7332b4326" @@ -28,6 +29,7 @@ GeometryOpsProjExt = "Proj" [compat] CoordinateTransformations = "0.5, 0.6" +DelaunayTriangulation = "1.0.4" ExactPredicates = "2.2.8" DimensionalData = "0.27" FlexiJoins = "0.1.30" diff --git a/README.md b/README.md index c3e9f9244..4c6fd6d53 100644 --- a/README.md +++ b/README.md @@ -1,7 +1,7 @@ GeometryOps.jl ![Lifecycle:Experimental](https://img.shields.io/badge/Lifecycle-Experimental-339999) -[![Stable](https://img.shields.io/badge/docs-stable-blue.svg)](https://JuliaGeo.github.io/GeometryOps.jl/v0.1.10/) +[![Stable](https://img.shields.io/badge/docs-stable-blue.svg)](https://JuliaGeo.github.io/GeometryOps.jl/stable/) [![Dev](https://img.shields.io/badge/docs-dev-blue.svg)](https://JuliaGeo.github.io/GeometryOps.jl/dev/) [![Build Status](https://github.com/JuliaGeo/GeometryOps.jl/actions/workflows/CI.yml/badge.svg?branch=main)](https://github.com/JuliaGeo/GeometryOps.jl/actions/workflows/CI.yml?query=branch%3Amain) diff --git a/benchmarks/benchmarks.jl b/benchmarks/benchmarks.jl index bd28d63df..af5f15ed9 100644 --- a/benchmarks/benchmarks.jl +++ b/benchmarks/benchmarks.jl @@ -11,7 +11,7 @@ import GeometryOps as GO, GeoInterface as GI, GeometryBasics as GB, LibGEOS as LG -import GeoJSON +import GeoJSON, NaturalEarth # In order to benchmark, we'll actually use the new [Chairmarks.jl](https://github.com/lilithhafner/Chairmarks.jl), # since it's significantly faster than BenchmarkTools. To keep benchmarks organized, though, we'll still use BenchmarkTools' # `BenchmarkGroup` structure. @@ -184,7 +184,7 @@ usa_difference_suite = usa_poly_suite["difference"] = BenchmarkGroup(["title:USA usa_intersection_suite = usa_poly_suite["intersection"] = BenchmarkGroup(["title:USA intersection", "subtitle:Tested on CONUS"]) usa_union_suite = usa_poly_suite["union"] = BenchmarkGroup(["title:USA union", "subtitle:Tested on CONUS"]) -fc = GeoJSON.read(read(download("https://rawcdn.githack.com/nvkelso/natural-earth-vector/ca96624a56bd078437bca8184e78163e5039ad19/geojson/ne_10m_admin_0_countries.geojson"))) +fc = NaturalEarth.naturalearth("admin_0_countries", 10) usa_multipoly = fc.geometry[findfirst(==("United States of America"), fc.NAME)] |> x -> GI.convert(LG, x) |> LG.makeValid |> GO.tuples usa_poly = GI.getgeom(usa_multipoly, findmax(GO.area.(GI.getgeom(usa_multipoly)))[2]) # isolate the poly with the most area @@ -194,7 +194,7 @@ f, a, p = plot(usa_poly; label = "Original", axis = (; aspect = DataAspect())); axislegend(a) f -# Now, we get to benchmarking: +# Now, we get to benchmarking: usa_o_lg, usa_o_go = lg_and_go(usa_poly) diff --git a/docs/src/tutorials/creating_geometry.md b/docs/src/tutorials/creating_geometry.md index a7838b8d9..dfcd2abf6 100644 --- a/docs/src/tutorials/creating_geometry.md +++ b/docs/src/tutorials/creating_geometry.md @@ -186,7 +186,7 @@ Read the land `MultiPolygon`s as a `GeoJSON.FeatureCollection`. land_geo = GeoJSON.read(land_path) ```` -We then need to create a figure with a `GeoAxis` that can handle the projection between `source` and `destinaton` CRS. For GeoMakie, `source` is the CRS of the input and `dest` is the CRS you want to visualize in. +We then need to create a figure with a `GeoAxis` that can handle the projection between `source` and `destination` CRS. For GeoMakie, `source` is the CRS of the input and `dest` is the CRS you want to visualize in. ````@example creating_geometry fig = Figure(size=(1000, 500)); @@ -267,7 +267,7 @@ y = r .* (k + 1) .* sin.(ϴ) .- r .* sin.((k + 1) .* ϴ); ring4 = GI.LinearRing(Point.(zip(x, y))) ```` -But this time when we create the `Polygon` we beed to specify the `CRS` at the time of creation, making it a geospatial polygon +But this time when we create the `Polygon` we need to specify the `CRS` at the time of creation, making it a geospatial polygon ````@example creating_geometry geopoly1 = GI.Polygon([ring4], crs = source_crs1) diff --git a/docs/src/tutorials/geodesic_paths.md b/docs/src/tutorials/geodesic_paths.md new file mode 100644 index 000000000..f6f4bf653 --- /dev/null +++ b/docs/src/tutorials/geodesic_paths.md @@ -0,0 +1,17 @@ +# Geodesic paths + +Geodesic paths are paths computed on an ellipsoid, as opposed to a plane. + +```@example geodesic +import GeometryOps as GO, GeoInterface as GI +using CairoMakie, GeoMakie + + +IAH = (-95.358421, 29.749907) +AMS = (4.897070, 52.377956) + + +fig, ga, _cp = lines(GeoMakie.coastlines(); axis = (; type = GeoAxis)) +lines!(ga, GO.segmentize(GO.GeodesicSegments(; max_distance = 100_000), GI.LineString([IAH, AMS])); color = Makie.wong_colors()[2]) +fig +``` \ No newline at end of file diff --git a/ext/GeometryOpsLibGEOSExt/simple_overrides.jl b/ext/GeometryOpsLibGEOSExt/simple_overrides.jl index 0836469f8..824534ec7 100644 --- a/ext/GeometryOpsLibGEOSExt/simple_overrides.jl +++ b/ext/GeometryOpsLibGEOSExt/simple_overrides.jl @@ -64,3 +64,17 @@ function GO.intersects(::GEOS, geom_a, geom_b) return LG.intersects(GI.convert(LG, geom_a), GI.convert(LG, geom_b)) end +# ## Convex hull +function GO.convex_hull(::GEOS, geoms) + return LG.convexhull( + LG.MultiPoint( + collect( + GO.flatten( + x -> GI.convert(LG, x), + GI.PointTrait, + geoms + ) + ) + ) + ) +end \ No newline at end of file diff --git a/ext/GeometryOpsProjExt/reproject.jl b/ext/GeometryOpsProjExt/reproject.jl index aecd0646b..616e13fd0 100644 --- a/ext/GeometryOpsProjExt/reproject.jl +++ b/ext/GeometryOpsProjExt/reproject.jl @@ -14,7 +14,7 @@ function reproject(geom; end # If its still nothing, error - isnothing(source_crs) && throw(ArgumentError("geom has no crs attatched. Pass a `source_crs` keyword")) + isnothing(source_crs) && throw(ArgumentError("geom has no crs attached. Pass a `source_crs` keyword")) # Otherwise reproject reproject(geom, source_crs, target_crs; kw...) diff --git a/src/GeometryOps.jl b/src/GeometryOps.jl index 1aced0e82..af6628c70 100644 --- a/src/GeometryOps.jl +++ b/src/GeometryOps.jl @@ -4,13 +4,14 @@ module GeometryOps using GeoInterface using GeometryBasics -import Tables using LinearAlgebra, Statistics + +import Tables import GeometryBasics.StaticArrays +import DelaunayTriangulation # for convex hull and triangulation import ExactPredicates import Base.@kwdef - -using GeoInterface.Extents: Extents +import GeoInterface.Extents: Extents const GI = GeoInterface const GB = GeometryBasics @@ -28,6 +29,7 @@ include("methods/area.jl") include("methods/barycentric.jl") include("methods/buffer.jl") include("methods/centroid.jl") +include("methods/convex_hull.jl") include("methods/distance.jl") include("methods/equals.jl") include("methods/clipping/predicates.jl") diff --git a/src/methods/angles.jl b/src/methods/angles.jl index c7af8732b..ca130c895 100644 --- a/src/methods/angles.jl +++ b/src/methods/angles.jl @@ -39,7 +39,7 @@ This is computed differently for different geometries: - The angles of a point is an empty vector. - The angles of a single line segment is an empty vector. - The angles of a linestring or linearring is a vector of angles formed by the curve. - - The angles of a polygin is a vector of vectors of angles formed by each ring. + - The angles of a polygon is a vector of vectors of angles formed by each ring. - The angles of a multi-geometry collection is a vector of the angles of each of the sub-geometries as defined above. @@ -83,7 +83,7 @@ function _angles(::Type{T}, ::GI.LinearRingTrait, geom; interior = true) where T end #= The angles of a polygon is a vector of polygon angles. Note that if there are holes -within the polyogn, the angles will be listed after the exterior ring angles in order of the +within the polygon, the angles will be listed after the exterior ring angles in order of the holes. All angles, including the hole angles, are interior angles of the polygon.=# function _angles(::Type{T}, ::GI.PolygonTrait, geom) where T angles = _angles(T, GI.LinearRingTrait(), GI.getexterior(geom); interior = true) diff --git a/src/methods/area.jl b/src/methods/area.jl index ae4b00c5f..b98024bb8 100644 --- a/src/methods/area.jl +++ b/src/methods/area.jl @@ -29,7 +29,7 @@ lines!( f ``` The points are ordered in a counterclockwise fashion, which means that the signed area -is negative. If we reverse the order of the points, we get a postive area. +is negative. If we reverse the order of the points, we get a positive area. ```@example rect GO.signed_area(rect) # -1.0 ``` @@ -76,7 +76,7 @@ end signed_area(geom, [T = Float64])::T Returns the signed area of a single geometry, based on winding order. -This is computed slighly differently for different geometries: +This is computed slightly differently for different geometries: - The signed area of a point is always zero. - The signed area of a curve is always zero. diff --git a/src/methods/centroid.jl b/src/methods/centroid.jl index 9dd0e4167..207203f47 100644 --- a/src/methods/centroid.jl +++ b/src/methods/centroid.jl @@ -45,7 +45,7 @@ LineString or LinearRing when they are closed, for example as the interior hole of a polygon. The helper functions centroid_and_length and centroid_and_area are made -availible just in case the user also needs the area or length to decrease +available just in case the user also needs the area or length to decrease repeat computation. =# """ diff --git a/src/methods/clipping/clipping_processor.jl b/src/methods/clipping/clipping_processor.jl index aaf938d95..97484404c 100644 --- a/src/methods/clipping/clipping_processor.jl +++ b/src/methods/clipping/clipping_processor.jl @@ -74,10 +74,10 @@ poly_a, including its intersection points with poly_b. The information stored in PolyNode is needed for clipping using the Greiner-Hormann clipping algorithm. Note: After calling this function, a_list is not fully formed because the neighboring -indicies of the intersection points in b_list still need to be updated. Also we still have +indices of the intersection points in b_list still need to be updated. Also we still have not update the entry and exit flags for a_list. -The a_idx_list is a list of the indicies of intersection points in a_list. The value at +The a_idx_list is a list of the indices of intersection points in a_list. The value at index i of a_idx_list is the location in a_list where the ith intersection point lies. =# function _build_a_list(::Type{T}, poly_a, poly_b; exact) where T @@ -179,7 +179,7 @@ creates a vector of PolyNodes to represent poly_b. The information stored in eac is needed for clipping using the Greiner-Hormann clipping algorithm. Note: after calling this function, b_list is not fully updated. The entry/exit flags still -need to be updated. However, the neightbor value in a_list is now updated. +need to be updated. However, the neighbor value in a_list is now updated. =# function _build_b_list(::Type{T}, a_idx_list, a_list, n_b_intrs, poly_b) where T # Sort intersection points by insertion order in b_list @@ -609,8 +609,8 @@ _get_poly_type(::Type{T}) where T = #= _find_non_cross_orientation(a_list, b_list, a_poly, b_poly; exact) -For polygns with no crossing intersection points, either one polygon is inside of another, -or they are seperate polygons with no intersection (other than an edge or point). +For polygons with no crossing intersection points, either one polygon is inside of another, +or they are separate polygons with no intersection (other than an edge or point). Return two booleans that represent if a is inside b (potentially with shared edges / points) and visa versa if b is inside of a. @@ -643,6 +643,7 @@ function _add_holes_to_polys!(::Type{T}, return_polys, hole_iterator, remove_pol for i in 1:n_polys n_new_per_poly = 0 for curr_hole in Iterators.map(tuples, hole_iterator) # loop through all holes + curr_hole = _linearring(curr_hole) # loop through all pieces of original polygon (new pieces added to end of list) for j in Iterators.flatten((i:i, (n_polys + 1):(n_polys + n_new_per_poly))) curr_poly = return_polys[j] @@ -655,7 +656,7 @@ function _add_holes_to_polys!(::Type{T}, return_polys, hole_iterator, remove_pol append!(remove_poly_idx, falses(n_new_pieces)) n_new_per_poly += n_new_pieces end - if !on_ext && !out_ext # hole is completly within exterior + if !on_ext && !out_ext # hole is completely within exterior push!(curr_poly.geom, new_hole) else # hole is partially within and outside of polygon's exterior new_polys = difference(curr_poly_ext, new_hole_poly, T; target=GI.PolygonTrait()) @@ -669,7 +670,7 @@ function _add_holes_to_polys!(::Type{T}, return_polys, hole_iterator, remove_pol n_new_per_poly += n_new_polys end end - # polygon is completly within hole + # polygon is completely within hole elseif coveredby(curr_poly_ext, GI.Polygon(StaticArrays.SVector(curr_hole))) remove_poly_idx[j] = true end @@ -687,7 +688,7 @@ end The new hole is combined with any existing holes in curr_poly. The holes can be combined into a larger hole if they are intersecting. If this happens, then the new, combined hole is -returned with the orignal holes making up the new hole removed from curr_poly. Additionally, +returned with the original holes making up the new hole removed from curr_poly. Additionally, if the combined holes form a ring, the interior is added to the return_polys as a new polygon piece. Additionally, holes leftover after combination will be checked for it they are in the "main" polygon or in one of these new pieces and moved accordingly. @@ -750,7 +751,7 @@ function _remove_collinear_points!(polys, remove_idx, poly_a, poly_b) continue else p3 = p - # check if p2 is approximatly on the edge formed by p1 and p3 - remove if so + # check if p2 is approximately on the edge formed by p1 and p3 - remove if so if Predicates.orient(p1, p2, p3; exact = _False()) == 0 remove_idx[i - 1] = true end diff --git a/src/methods/clipping/coverage.jl b/src/methods/clipping/coverage.jl index 597fb470f..cf8762274 100644 --- a/src/methods/clipping/coverage.jl +++ b/src/methods/clipping/coverage.jl @@ -4,7 +4,7 @@ export coverage ## What is coverage? Coverage is the amount of geometry area within a bounding box defined by the minimum and -maximum x and y-coordiantes of that bounding box, or an Extent containing that information. +maximum x and y-coordinates of that bounding box, or an Extent containing that information. To provide an example, consider this rectangle: ```@example rect diff --git a/src/methods/clipping/cut.jl b/src/methods/clipping/cut.jl index 1303ceb7e..be3f0a855 100644 --- a/src/methods/clipping/cut.jl +++ b/src/methods/clipping/cut.jl @@ -98,7 +98,7 @@ of cut geometry in Vector{Vector{Tuple}} format. Note: degenerate cases where intersection points are vertices do not work right now. =# function _cut(::Type{T}, geom, line, geom_list, intr_list, n_intr_pts; exact) where T - # Sort and catagorize the intersection points + # Sort and categorize the intersection points sort!(intr_list, by = x -> geom_list[x].fracs[2]) _flag_ent_exit!(GI.LineTrait(), line, geom_list; exact) # Add first point to output list diff --git a/src/methods/clipping/difference.jl b/src/methods/clipping/difference.jl index 604e4606b..d121e56c5 100644 --- a/src/methods/clipping/difference.jl +++ b/src/methods/clipping/difference.jl @@ -83,7 +83,7 @@ function _difference( end end end - # Remove uneeded collinear points on same edge + # Remove unneeded collinear points on same edge _remove_collinear_points!(polys, remove_idx, poly_a, poly_b) return polys end @@ -147,7 +147,7 @@ end #= Multipolygon with multipolygon difference - note that all intersection regions between sub-polygons of `multipoly_a` and sub-polygons of `multipoly_b` will be removed from the corresponding sub-polygon of `multipoly_a`. Unless specified with `fix_multipoly = nothing`, -`multipolygon_a` will be validated using the given (defauly is `UnionIntersectingPolygons()`) +`multipolygon_a` will be validated using the given (default is `UnionIntersectingPolygons()`) correction. =# function _difference( target::TraitTarget{GI.PolygonTrait}, ::Type{T}, @@ -169,7 +169,7 @@ function _difference( else difference(GI.MultiPolygon(polys), poly_b; target, fix_multipoly) end - #= One multipoly_a has been completly covered (and thus removed) there is no need to + #= One multipoly_a has been completely covered (and thus removed) there is no need to continue taking the difference =# isempty(polys) && break end diff --git a/src/methods/clipping/intersection.jl b/src/methods/clipping/intersection.jl index 405be4428..5892bd968 100644 --- a/src/methods/clipping/intersection.jl +++ b/src/methods/clipping/intersection.jl @@ -5,7 +5,7 @@ export intersection, intersection_points Enum LineOrientation Enum for the orientation of a line with respect to a curve. A line can be `line_cross` (crossing over the curve), `line_hinge` (crossing the endpoint of the curve), -`line_over` (colinear with the curve), or `line_out` (not interacting with the curve). +`line_over` (collinear with the curve), or `line_out` (not interacting with the curve). """ @enum LineOrientation line_cross=1 line_hinge=2 line_over=3 line_out=4 @@ -86,7 +86,7 @@ function _intersection( hole_iterator = Iterators.flatten((GI.gethole(poly_a), GI.gethole(poly_b))) _add_holes_to_polys!(T, polys, hole_iterator, remove_idx; exact) end - # Remove uneeded collinear points on same edge + # Remove unneeded collinear points on same edge _remove_collinear_points!(polys, remove_idx, poly_a, poly_b) return polys end @@ -128,7 +128,7 @@ function _intersection( end #= Multipolygon with polygon intersection is equivalent to taking the intersection of the -poylgon with the multipolygon and thus simply switches the order of operations and calls the +polygon with the multipolygon and thus simply switches the order of operations and calls the above method. =# _intersection( target::TraitTarget{GI.PolygonTrait}, ::Type{T}, @@ -197,7 +197,7 @@ intersection_points(geom_a, geom_b, ::Type{T} = Float64) where T <: AbstractFloa _intersection_points(T, GI.trait(geom_a), geom_a, GI.trait(geom_b), geom_b) -#= Calculates the list of intersection points between two geometries, inlcuding line +#= Calculates the list of intersection points between two geometries, including line segments, line strings, linear rings, polygons, and multipolygons. =# function _intersection_points(::Type{T}, ::GI.AbstractTrait, a, ::GI.AbstractTrait, b; exact = _True()) where T # Initialize an empty list of points @@ -384,7 +384,7 @@ end #= If lines defined by (a1, a2) and (b1, b2) meet at one point that is not an endpoint of either segment, they form a crossing intersection with a singular intersection point. That -point is caculated by finding the fractional distance along each segment the point occurs +point is calculated by finding the fractional distance along each segment the point occurs at (α, β). If the point is too close to an endpoint to be distinct, the point shares a value with the endpoint, but with a non-zero and non-one fractional value. If the intersection point calculated is outside of the envelope of the two segments due to floating point error, @@ -407,9 +407,9 @@ function _find_cross_intersection(::Type{T}, a1, a2, b1, b2, a_ext, b_ext) where β = _clamped_frac(Δbax * Δay - Δbay * Δax, a_cross_b, eps(T)) #= Intersection will be where a1 + α * Δa = b1 + β * Δb. However, due to floating point - innacurracies, α and β calculations may yeild different intersection points. Average + inaccuracies, α and β calculations may yield different intersection points. Average both points together to minimize difference from real value, as long as segment isn't - vertical or horizontal as this will almost certianly lead to the point being outside the + vertical or horizontal as this will almost certainly lead to the point being outside the envelope due to floating point error. Also note that floating point limitations could make intersection be endpoint if α≈0 or α≈1.=# x = if Δax == 0 diff --git a/src/methods/clipping/union.jl b/src/methods/clipping/union.jl index 8e1d785e4..e4b4bda8e 100644 --- a/src/methods/clipping/union.jl +++ b/src/methods/clipping/union.jl @@ -62,9 +62,9 @@ function _union( if n_pieces == 0 # no crossing points, determine if either poly is inside the other a_in_b, b_in_a = _find_non_cross_orientation(a_list, b_list, ext_a, ext_b; exact) if a_in_b - push!(polys, GI.Polygon([tuples(ext_b)])) + push!(polys, GI.Polygon([_linearring(tuples(ext_b))])) elseif b_in_a - push!(polys, GI.Polygon([tuples(ext_a)])) + push!(polys, GI.Polygon([_linearring(tuples(ext_a))])) else push!(polys, tuples(poly_a)) push!(polys, tuples(poly_b)) @@ -83,7 +83,7 @@ function _union( if GI.nhole(poly_a) != 0 || GI.nhole(poly_b) != 0 _add_union_holes!(polys, a_in_b, b_in_a, poly_a, poly_b; exact) end - # Remove uneeded collinear points on same edge + # Remove unneeded collinear points on same edge _remove_collinear_points!(polys, [false], poly_a, poly_b) return polys end @@ -124,6 +124,7 @@ function _add_union_holes!(polys, a_in_b, b_in_a, poly_a, poly_b; exact) current_poly = n_a_holes > 0 ? ext_poly_b : poly_a # Loop over all holes in both original polygons for (i, ih) in enumerate(Iterators.flatten((GI.gethole(poly_a), GI.gethole(poly_b)))) + ih = _linearring(ih) in_ext, _, _ = _line_polygon_interactions(ih, curr_exterior_poly; exact, closed_line = true) if !in_ext #= if the hole isn't in the overlapping region between the two polygons, add @@ -160,7 +161,7 @@ function _add_union_holes_contained_polys!(polys, interior_poly, exterior_poly; in_ih, on_ih, out_ih = _line_polygon_interactions(ext_int_ring, poly_ih; exact, closed_line = true) if in_ih # at least part of interior polygon exterior is within the ith hole if !on_ih && !out_ih - #= interior polygon is completly within the ith hole - polygons aren't + #= interior polygon is completely within the ith hole - polygons aren't touching and do not actually form a union =# polys[1] = tuples(interior_poly) push!(polys, tuples(exterior_poly)) @@ -235,7 +236,7 @@ function _union( return polys end -#= Multipolygon with polygon union is equivalent to taking the union of the poylgon with the +#= Multipolygon with polygon union is equivalent to taking the union of the polygon with the multipolygon and thus simply switches the order of operations and calls the above method. =# _union( target::TraitTarget{GI.PolygonTrait}, ::Type{T}, diff --git a/src/methods/convex_hull.jl b/src/methods/convex_hull.jl new file mode 100644 index 000000000..429b5281a --- /dev/null +++ b/src/methods/convex_hull.jl @@ -0,0 +1,131 @@ +#= +# Convex hull + +The [_convex hull_](https://en.wikipedia.org/wiki/Convex_hull) of a set of points is the smallest [**convex**](https://en.wikipedia.org/wiki/Convex_set) polygon that contains all the points. + +GeometryOps.jl provides a number of methods for computing the convex hull of a set of points, usually +linked to other Julia packages. + +For now, we expose one algorithm, [MonotoneChainMethod](@ref), which uses the [DelaunayTriangulation.jl](https://github.com/JuliaGeometry/DelaunayTriangulation.jl) +package. The `GEOS()` interface also supports convex hulls. + +Future work could include other algorithms, such as [Quickhull.jl](https://github.com/augustt198/Quickhull.jl), or similar, via package extensions. + + +## Example + +### Simple hull +```@example simple +import GeometryOps as GO, GeoInterface as GI +using CairoMakie # to plot + +points = randn(GO.Point2f, 100) +f, a, p = plot(points; label = "Points") +hull_poly = GO.convex_hull(points) +lines!(a, hull_poly; label = "Convex hull", color = Makie.wong_colors()[2]) +axislegend(a) +f +``` + +## Convex hull of the USA +```@example usa +import GeometryOps as GO, GeoInterface as GI +using CairoMakie # to plot +using NaturalEarth # for data + +all_adm0 = naturalearth("admin_0_countries", 110) +usa = all_adm0.geometry[findfirst(==("USA"), all_adm0.ADM0_A3)] +f, a, p = lines(usa) +lines!(a, GO.convex_hull(usa); color = Makie.wong_colors()[2]) +f +``` + +## Investigating the winding order + +The winding order of the monotone chain method is counterclockwise, +while the winding order of the GEOS method is clockwise. + +GeometryOps' convexity detection says that the GEOS hull is convex, +while the monotone chain method hull is not. However, they are both going +over the same points (we checked), it's just that the winding order is different. + +In reality, both sets are convex, but we need to fix the GeometryOps convexity detector +([`isconcave`](@ref))! + +We may also decide at a later date to change the returned winding order of the polygon, but +most algorithms are robust to that, and you can always [`fix`](@ref) it... + +```@example windingorder +import GeoInterface as GI, GeometryOps as GO, LibGEOS as LG +using CairoMakie # to plot + +points = rand(Point2{Float64}, 100) +go_hull = GO.convex_hull(GO.MonotoneChainMethod(), points) +lg_hull = GO.convex_hull(GO.GEOS(), points) + +fig = Figure() +a1, p1 = lines(fig[1, 1], go_hull; color = 1:GI.npoint(go_hull), axis = (; title = "MonotoneChainMethod()")) +a2, p2 = lines(fig[2, 1], lg_hull; color = 1:GI.npoint(lg_hull), axis = (; title = "GEOS()")) +cb = Colorbar(fig[1:2, 2], p1; label = "Vertex number") +fig +``` + +## Implementation + +=# + +""" + convex_hull([method], geometries) + +Compute the convex hull of the points in `geometries`. +Returns a `GI.Polygon` representing the convex hull. + +Note that the polygon returned is wound counterclockwise +as in the Simple Features standard by default. If you +choose GEOS, the winding order will be inverted. + +!!! warning + This interface only computes the 2-dimensional convex hull! + + For higher dimensional hulls, use the relevant package (Qhull.jl, Quickhull.jl, or similar). +""" +function convex_hull end + +""" + MonotoneChainMethod() + +This is an algorithm for the [`convex_hull`](@ref) function. + +Uses [`DelaunayTriangulation.jl`](https://github.com/JuliaGeometry/DelaunayTriangulation.jl) to compute the convex hull. +This is a pure Julia algorithm which provides an optimal Delaunay triangulation. + +See also [`convex_hull`](@ref) +""" +struct MonotoneChainMethod end + +# GrahamScanMethod, etc. can be implemented in GO as well, if someone wants to. +# If we add an extension on Quickhull.jl, then that would be another algorithm. + +convex_hull(geometries) = convex_hull(MonotoneChainMethod(), geometries) + +# TODO: have this respect the CRS by pulling it out of `geometries`. +function convex_hull(::MonotoneChainMethod, geometries) + # Extract all points as tuples. We have to collect and allocate + # here, because DelaunayTriangulation only accepts vectors of + # point-like geoms. + + # Cleanest would be to use the iterable from GO.flatten directly, + # but that would require us to implement the convex hull algorithm + # directly. + + # TODO: create a specialized method that extracts only the information + # required, GeometryBasics points can be passed through directly. + points = collect(flatten(tuples, GI.PointTrait, geometries)) + # Compute the convex hull using DelTri (shorthand for DelaunayTriangulation.jl). + hull = DelaunayTriangulation.convex_hull(points) + # Convert the result to a `GI.Polygon` and return it. + # View would be more efficient here, but re-allocating + # is cleaner. + point_vec = DelaunayTriangulation.get_points(hull)[DelaunayTriangulation.get_vertices(hull)] + return GI.Polygon([GI.LinearRing(point_vec)]) +end diff --git a/src/methods/distance.jl b/src/methods/distance.jl index ef470330f..0333f2916 100644 --- a/src/methods/distance.jl +++ b/src/methods/distance.jl @@ -12,7 +12,7 @@ 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. +either be positive or 0. To provide an example, consider this rectangle: ```@example rect diff --git a/src/methods/equals.jl b/src/methods/equals.jl index 15761d63e..cf0efe11d 100644 --- a/src/methods/equals.jl +++ b/src/methods/equals.jl @@ -23,7 +23,7 @@ lines!(GI.getpoint(l2), color = :orange) scatter!(GI.getpoint(l2), color = :orange) f ``` -We can see that the two lines do not share a commen set of points and edges in +We can see that the two lines do not share a common set of points and edges in the plot, so they are not equal: ```@example equals GO.equals(l1, l2) # returns false @@ -44,7 +44,7 @@ same order. The winding order also doesn't have to be the same to represent the same geometry. This requires checking every point against every other point in the two geometries we are comparing. Also, some geometries must be "closed" like polygons and linear rings. These will be assumed to be closed, even if they -don't have a repeated last point explicity written in the coordinates. +don't have a repeated last point explicitly written in the coordinates. Additionally, geometries and multi-geometries can be equal if the multi-geometry only includes that single geometry. =# @@ -210,7 +210,7 @@ end )::Bool Two lines/linestrings are equal if they share the same set of points going -along the curve. Note that lines/linestrings aren't closed by defintion. +along the curve. Note that lines/linestrings aren't closed by definition. """ equals( ::Union{GI.LineTrait, GI.LineStringTrait}, l1, @@ -224,7 +224,7 @@ equals( )::Bool A line/linestring and a linear ring are equal if they share the same set of -points going along the curve. Note that lines aren't closed by defintion, but +points going along the curve. Note that lines aren't closed by definition, but rings are, so the line must have a repeated last point to be equal """ equals( @@ -239,7 +239,7 @@ equals( )::Bool A linear ring and a line/linestring are equal if they share the same set of -points going along the curve. Note that lines aren't closed by defintion, but +points going along the curve. Note that lines aren't closed by definition, but rings are, so the line must have a repeated last point to be equal """ equals( diff --git a/src/methods/geom_relations/contains.jl b/src/methods/geom_relations/contains.jl index 79ec841e1..35d0c6e60 100644 --- a/src/methods/geom_relations/contains.jl +++ b/src/methods/geom_relations/contains.jl @@ -5,8 +5,8 @@ export contains #= ## What is contains? -The contains function checks if a given geometry completly contains another -geometry, or in other words, that the second geometry is completly within the +The contains function checks if a given geometry completely contains another +geometry, or in other words, that the second geometry is completely within the first. This requires that the two interiors intersect and that the interior and boundary of the second geometry is not in the exterior of the first geometry. diff --git a/src/methods/geom_relations/coveredby.jl b/src/methods/geom_relations/coveredby.jl index a74cd6568..ff5ed7f06 100644 --- a/src/methods/geom_relations/coveredby.jl +++ b/src/methods/geom_relations/coveredby.jl @@ -160,7 +160,7 @@ _coveredby( ) #= Linestring is coveredby a polygon if all interior and boundary points of the -line are in the polygon interior or on its edges, inlcuding hole edges. =# +line are in the polygon interior or on its edges, including hole edges. =# _coveredby( ::Union{GI.LineTrait, GI.LineStringTrait}, g1, ::GI.PolygonTrait, g2, @@ -203,7 +203,7 @@ _coveredby( ) #= Linearring is coveredby a polygon if all vertices and edges of the ring are -in the polygon interior or on the polygon edges, inlcuding hole edges. =# +in the polygon interior or on the polygon edges, including hole edges. =# _coveredby( ::GI.LinearRingTrait, g1, ::GI.PolygonTrait, g2, diff --git a/src/methods/geom_relations/covers.jl b/src/methods/geom_relations/covers.jl index a8e067a10..47e99dead 100644 --- a/src/methods/geom_relations/covers.jl +++ b/src/methods/geom_relations/covers.jl @@ -5,7 +5,7 @@ export covers #= ## What is covers? -The covers function checks if a given geometry completly covers another +The covers function checks if a given geometry completely covers another geometry. For this to be true, the "contained" geometry's interior and boundaries must be covered by the "covering" geometry's interior and boundaries. The interiors do not need to overlap. diff --git a/src/methods/geom_relations/disjoint.jl b/src/methods/geom_relations/disjoint.jl index 58774d03d..dfb06a838 100644 --- a/src/methods/geom_relations/disjoint.jl +++ b/src/methods/geom_relations/disjoint.jl @@ -171,7 +171,7 @@ _disjoint( ) #= Geometry is disjoint from a linestring if the line's interior and boundary -points don't intersect with the geometrie's interior and boundary points. =# +points don't intersect with the geometry's interior and boundary points. =# _disjoint( trait1::Union{GI.LinearRingTrait, GI.PolygonTrait}, g1, trait2::Union{GI.LineTrait, GI.LineStringTrait}, g2, diff --git a/src/methods/geom_relations/geom_geom_processors.jl b/src/methods/geom_relations/geom_geom_processors.jl index 9df829974..a0b9d7e99 100644 --- a/src/methods/geom_relations/geom_geom_processors.jl +++ b/src/methods/geom_relations/geom_geom_processors.jl @@ -96,12 +96,12 @@ If cross_allow is true, segments of the line and curve can be disjoint. If in_require is true, the interiors of the line and curve must meet in at least one point. -If on_require is true, the bounday of one of the two geometries can meet the +If on_require is true, the boundary of one of the two geometries can meet the interior or boundary of the other geometry in at least one point. If out_require is true, there must be at least one point of the given line that is exterior of the curve. -If the point is in an "allowed" location and meets all requirments, return true. +If the point is in an "allowed" location and meets all requirements, return true. Else, return false. If closed_line is true, line is treated as a closed line where the first and @@ -125,7 +125,7 @@ function _inner_line_curve_process( closed_line = false, closed_curve = false, exact, ) - # Set up requirments + # Set up requirements in_req_met = !in_require on_req_met = !on_require out_req_met = !out_require @@ -152,7 +152,7 @@ function _inner_line_curve_process( # If segments are co-linear if seg_val == line_over !over_allow && return false - # at least one point in, meets requirments + # at least one point in, meets requirements in_req_met = true point_val = _point_segment_orientation(l_start, c_start, c_end) # If entire segment isn't covered, consider remaining section @@ -208,7 +208,7 @@ function _inner_line_curve_process( end #= If entire segment (le to ls) isn't covered by segment (cs to ce), find remaining section -part of section outside of cs to ce. If completly covered, increase segment index i. =# +part of section outside of cs to ce. If completely covered, increase segment index i. =# function _find_new_seg(i, ls, le, cs, ce) break_off = true if _point_segment_orientation(le, cs, ce) != point_out @@ -251,7 +251,7 @@ If on_require is true, the line must have at least one point on the polygon'same If out_require is true, the line must have at least one point outside of the polygon. -If the point is in an "allowed" location and meets all requirments, return true. +If the point is in an "allowed" location and meets all requirements, return true. Else, return false. If closed_line is true, line is treated as a closed line where the first and @@ -298,7 +298,7 @@ function _inner_line_polygon_process( !out_allow && return false out_req_met = true end - if on_hole # hole bounday is polygon boundary + if on_hole # hole boundary is polygon boundary !on_allow && return false on_req_met = true end @@ -329,7 +329,7 @@ If on_require is true, one of the polygon's must have at least one boundary If out_require is true, the first polygon must have at least one interior point outside of the second polygon. -If the point is in an "allowed" location and meets all requirments, return true. +If the point is in an "allowed" location and meets all requirements, return true. Else, return false. =# @inline function _polygon_polygon_process(poly1, poly2; kw...) @@ -482,7 +482,7 @@ for each scenario. Note that this uses the Algorithm by Hao and Sun (2018): https://doi.org/10.3390/sym10100477 -Paper seperates orientation of point and edge into 26 cases. For each case, it +Paper separates orientation of point and edge into 26 cases. For each case, it is either a case where the point is on the edge (returns on), where a ray from the point (x, y) to infinity along the line y = y cut through the edge (k += 1), or the ray does not pass through the edge (do nothing and continue). If the ray @@ -643,7 +643,7 @@ Returns a tuple of booleans: (in_poly, on_poly, out_poly). If in_poly is true, some of the lines interior points interact with the polygon interior points. If in_poly is true, endpoints of either the line intersect with the polygon or - the line interacts with the polygon boundary, including hole bounaries. + the line interacts with the polygon boundary, including hole boundaries. If out_curve is true, at least one segments of the line is outside the polygon, including inside of holes. diff --git a/src/methods/geom_relations/overlaps.jl b/src/methods/geom_relations/overlaps.jl index 4bbdbeb21..8cb966a0b 100644 --- a/src/methods/geom_relations/overlaps.jl +++ b/src/methods/geom_relations/overlaps.jl @@ -11,7 +11,7 @@ contained, within, or equal to the other. Note that this means it is impossible for a single point to overlap with a single point and a line only overlaps with another line if only a section of -each line is colinear. +each line is collinear. To provide an example, consider these two lines: ```@example overlaps @@ -42,12 +42,12 @@ implementation based on the geometry trait. This is also used in the implementation, since it's a lot less work! Note that that since only elements of the same dimension can overlap, any two -geometries with traits that are of different dimensions autmoatically can +geometries with traits that are of different dimensions automatically can return false. For geometries with the same trait dimension, we must make sure that they share a point, an edge, or area for points, lines, and polygons/multipolygons -respectivly, without being contained. +respectively, without being contained. =# """ @@ -116,7 +116,7 @@ end """ overlaps(::GI.LineTrait, line1, ::GI.LineTrait, line)::Bool -If the lines overlap, meaning that they are colinear but each have one endpoint +If the lines overlap, meaning that they are collinear but each have one endpoint outside of the other line, return true. Else false. """ overlaps(::GI.LineTrait, line1, ::GI.LineTrait, line) = @@ -212,7 +212,7 @@ function overlaps( return false end -#= If the edges overlap, meaning that they are colinear but each have one endpoint +#= If the edges overlap, meaning that they are collinear but each have one endpoint outside of the other edge, return true. Else false. =# function _overlaps( (a1, a2)::Edge, diff --git a/src/methods/geom_relations/touches.jl b/src/methods/geom_relations/touches.jl index 1d9018bc4..43a08e7a1 100644 --- a/src/methods/geom_relations/touches.jl +++ b/src/methods/geom_relations/touches.jl @@ -8,7 +8,7 @@ export touches The touches function checks if one geometry touches another geometry. In other words, the interiors of the two geometries don't interact, but one of the geometries must have a boundary point that interacts with either the other -geometies interior or boundary. +geometry's interior or boundary. To provide an example, consider these two lines: @@ -61,7 +61,7 @@ const TOUCHES_EXACT = (exact = _False(),) Return `true` if the first geometry touches the second geometry. In other words, the two interiors cannot interact, but one of the geometries must have a -boundary point that interacts with either the other geometies interior or +boundary point that interacts with either the other geometry's interior or boundary. ## Examples @@ -128,8 +128,8 @@ _touches( # # Lines touching geometries -#= Linestring touches another line if at least one bounday point interacts with -the bounday of interior of the other line, but the interiors don't interact. =# +#= Linestring touches another line if at least one boundary point interacts with +the boundary of interior of the other line, but the interiors don't interact. =# _touches( ::Union{GI.LineTrait, GI.LineStringTrait}, g1, ::Union{GI.LineTrait, GI.LineStringTrait}, g2, @@ -181,14 +181,14 @@ _touches( ) = _touches(trait2, g2, trait1, g1) #= Linearring cannot touch another linear ring since they are both exclusively -made up of interior points and no bounday points =# +made up of interior points and no boundary points =# _touches( ::GI.LinearRingTrait, g1, ::GI.LinearRingTrait, g2, ) = false #= Linearring touches a polygon if at least one of the points of the ring -interact with the polygon bounday and non are in the polygon interior. =# +interact with the polygon boundary and non are in the polygon interior. =# _touches( ::GI.LinearRingTrait, g1, ::GI.PolygonTrait, g2, @@ -203,8 +203,8 @@ _touches( # # Polygons touch geometries -#= Polygon touches a curve if at least one of the curve bounday points interacts -with the polygon's bounday and no curve points interact with the interior.=# +#= Polygon touches a curve if at least one of the curve boundary points interacts +with the polygon's boundary and no curve points interact with the interior.=# _touches( trait1::GI.PolygonTrait, g1, trait2::GI.AbstractCurveTrait, g2 diff --git a/src/methods/geom_relations/within.jl b/src/methods/geom_relations/within.jl index 5b240b819..8d2e34b16 100644 --- a/src/methods/geom_relations/within.jl +++ b/src/methods/geom_relations/within.jl @@ -40,14 +40,14 @@ implementation based on the geometry trait. Each of these calls a method in the geom_geom_processors file. The methods in this file determine if the given geometries meet a set of criteria. For the -`within` function and arguments g1 and g2, this criteria is as follows: - - points of g1 are allowed to be in the interior of g2 (either through +`within` function and arguments `g1` and `g2`, this criteria is as follows: + - points of `g1` are allowed to be in the interior of `g2` (either through overlap or crossing for lines) - - points of g1 are allowed to be on the boundary of g2 - - points of g1 are not allowed to be in the exterior of g2 - - at least one point of g1 is required to be in the interior of g2 - - no points of g1 are required to be on the boundary of g2 - - no points of g1 are required to be in the exterior of g2 + - points of `g1` are allowed to be on the boundary of `g2` + - points of `g1` are not allowed to be in the exterior of `g2` + - at least one point of `g1` is required to be in the interior of `g2` + - no points of `g1` are required to be on the boundary of `g2` + - no points of `g1` are required to be in the exterior of `g2` The code for the specific implementations is in the geom_geom_processors file. =# diff --git a/src/methods/polygonize.jl b/src/methods/polygonize.jl index b238f4c63..49e7c5f6c 100644 --- a/src/methods/polygonize.jl +++ b/src/methods/polygonize.jl @@ -35,7 +35,7 @@ f Now, we can use the `polygonize` function to convert the raster data into polygons. For this particular example, we chose a range of z-values between 0.8 and 3.2, -which would provide two distinct polyogns with holes. +which would provide two distinct polygons with holes. ```@example polygonize polygons = polygonize(xs, ys, 0.8 .< zs .< 3.2) @@ -117,12 +117,12 @@ function _polygonize(f::Base.Callable, xs::AbstractRange, ys::AbstractRange, A:: xhalf = step(xs) / 2 yhalf = step(ys) / 2 # Make bounds ranges first to avoid floating point error making gaps or overlaps - xbounds = first(xs) - xhalf : step(xs) : last(xs) + xhalf - ybounds = first(ys) - yhalf : step(ys) : last(ys) + yhalf + xbounds = range(first(xs) - xhalf; step = step(xs), length = length(xs) + 1) + ybounds = range(first(ys) - yhalf; step = step(ys), length = length(ys) + 1) Tx = eltype(xbounds) Ty = eltype(ybounds) - xvec = similar(Vector{Tuple{Tx,Tx}}, xs) - yvec = similar(Vector{Tuple{Ty,Ty}}, ys) + xvec = similar(Vector{Tuple{Tx,Tx}}, length(xs)) + yvec = similar(Vector{Tuple{Ty,Ty}}, length(ys)) for (xind, i) in enumerate(eachindex(xvec)) xvec[i] = xbounds[xind], xbounds[xind+1] end @@ -334,7 +334,7 @@ function _polygonize(f, xs::AbstractVector{T}, ys::AbstractVector{T}, A::Abstrac assigned_holes == length(holes) || @warn "Not all holes were assigned to polygons, $(length(holes) - assigned_holes) where missed from $(length(holes)) holes and $(length(polygons)) polygons" if isempty(polygons) - # TODO: this really should return an emtpty MultiPolygon but + # TODO: this really should return an empty MultiPolygon but # GeoInterface wrappers cant do that yet, which is not ideal... @warn "No polgons found, check your data or try another function for `f`" return nothing diff --git a/src/primitives.jl b/src/primitives.jl index 01014c1cf..3a6af61a9 100644 --- a/src/primitives.jl +++ b/src/primitives.jl @@ -223,7 +223,7 @@ end if calc_extent isa _True # Calculate the extent of the features extent = mapreduce(GI.extent, Extents.union, features) - # Return a FeatureCollection with features, crs and caculated extent + # Return a FeatureCollection with features, crs and calculated extent return GI.FeatureCollection(features; crs, extent) else # Return a FeatureCollection with features and crs @@ -241,10 +241,10 @@ end if calc_extent isa _True # Calculate the extent of the geometry extent = GI.extent(geometry) - # Return a new Feature with the new geometry and calculated extent, but the oroginal properties and crs + # Return a new Feature with the new geometry and calculated extent, but the original properties and crs return GI.Feature(geometry; properties, crs, extent) else - # Return a new Feature with the new geometry, but the oroginal properties and crs + # Return a new Feature with the new geometry, but the original properties and crs return GI.Feature(geometry; properties, crs) end end @@ -260,15 +260,28 @@ end geoms = _maptasks(apply_to_geom, 1:GI.ngeom(geom), threaded) return _apply_inner(geom, geoms, crs, calc_extent) end +@inline function _apply(f::F, target::TraitTarget{<:PointTrait}, trait::GI.PolygonTrait, geom; + crs=GI.crs(geom), calc_extent=_False(), threaded +)::(GI.geointerface_geomtype(trait)) where F + # We need to force rebuilding a LinearRing not a LineString + geoms = _maptasks(1:GI.ngeom(geom), threaded) do i + lr = GI.getgeom(geom, i) + points = map(GI.getgeom(lr)) do p + _apply(f, target, p; crs, calc_extent, threaded=_False()) + end + _linearring(_apply_inner(lr, points, crs, calc_extent)) + end + return _apply_inner(geom, geoms, crs, calc_extent) +end function _apply_inner(geom, geoms, crs, calc_extent::_True) # Calculate the extent of the sub geometries extent = mapreduce(GI.extent, Extents.union, geoms) # Return a new geometry of the same trait as `geom`, - # holding tnew `geoms` with `crs` and calcualted extent + # holding the new `geoms` with `crs` and calculated extent return rebuild(geom, geoms; crs, extent) end function _apply_inner(geom, geoms, crs, calc_extent::_False) - # Return a new geometryof the same trait as `geom`, holding the new `geoms` with `crs` + # Return a new geometry of the same trait as `geom`, holding the new `geoms` with `crs` return rebuild(geom, geoms; crs) end # Fail loudly if we hit PointTrait without running `f` @@ -308,7 +321,7 @@ end @inline _applyreduce(f::F, op::O, target, geom; threaded, init) where {F, O} = _applyreduce(f, op, target, GI.trait(geom), geom; threaded, init) -# Maybe use threads recucing over arrays +# Maybe use threads reducing over arrays @inline function _applyreduce(f::F, op::O, target, ::Nothing, A::AbstractArray; threaded, init) where {F, O} applyreduce_array(i) = _applyreduce(f, op, target, A[i]; threaded=_False(), init) _mapreducetasks(applyreduce_array, op, eachindex(A), threaded; init) @@ -455,7 +468,7 @@ Reconstruct `geom` from an iterable of component objects that match its structur All objects in `components` must have the same `GeoInterface.trait`. -Ususally used in combination with `flatten`. +Usually used in combination with `flatten`. """ function reconstruct(geom, components) obj, iter = _reconstruct(geom, components) @@ -537,7 +550,7 @@ using Base.Threads: nthreads, @threads, @spawn # Threading utility, modified Mason Protters threading PSA -# run `f` over ntasks, where f recieves an AbstractArray/range +# run `f` over ntasks, where f receives an AbstractArray/range # of linear indices @inline function _maptasks(f::F, taskrange, threaded::_True)::Vector where F ntasks = length(taskrange) @@ -568,7 +581,7 @@ Base.@assume_effects :foldable @inline function _maptasks(f::F, taskrange, threa end # Threading utility, modified Mason Protters threading PSA -# run `f` over ntasks, where f recieves an AbstractArray/range +# run `f` over ntasks, where f receives an AbstractArray/range # of linear indices # # WARNING: this will not work for mean/median - only ops diff --git a/src/transformations/correction/intersecting_polygons.jl b/src/transformations/correction/intersecting_polygons.jl index e8858d45f..2223f6b89 100644 --- a/src/transformations/correction/intersecting_polygons.jl +++ b/src/transformations/correction/intersecting_polygons.jl @@ -5,7 +5,7 @@ export UnionIntersectingPolygons #= If the sub-polygons of a multipolygon are intersecting, this makes them invalid according to specification. Each sub-polygon of a multipolygon being disjoint (other than by a single -point) is a requirment for a valid multipolygon. However, different libraries may achieve +point) is a requirement for a valid multipolygon. However, different libraries may achieve this in different ways. For example, taking the union of all sub-polygons of a multipolygon will create a new @@ -141,4 +141,4 @@ function (::DiffIntersectingPolygons)(::GI.MultiPolygonTrait, multipoly) keepat!(diff_multipoly.geom, keep_idx) end return diff_multipoly -end \ No newline at end of file +end diff --git a/src/transformations/reproject.jl b/src/transformations/reproject.jl index 1b9dd6855..ab3707374 100644 --- a/src/transformations/reproject.jl +++ b/src/transformations/reproject.jl @@ -27,12 +27,12 @@ geometries, wrapping views of a `Vector{Proj.Point{D}}`, where `D` is the dimens ## Arguments - `geometry`: Any GeoInterface.jl compatible geometries. -- `source_crs`: the source coordinate referece system, as a GeoFormatTypes.jl object or a string. -- `target_crs`: the target coordinate referece system, as a GeoFormatTypes.jl object or a string. +- `source_crs`: the source coordinate reference system, as a GeoFormatTypes.jl object or a string. +- `target_crs`: the target coordinate reference system, as a GeoFormatTypes.jl object or a string. If these a passed as keywords, `transform` will take priority. Without it `target_crs` is always needed, and `source_crs` is -needed if it is not retreivable from the geometry with `GeoInterface.crs(geometry)`. +needed if it is not retrievable from the geometry with `GeoInterface.crs(geometry)`. ## Keywords diff --git a/src/transformations/simplify.jl b/src/transformations/simplify.jl index 35ccabdb1..8a36a48cc 100644 --- a/src/transformations/simplify.jl +++ b/src/transformations/simplify.jl @@ -37,10 +37,11 @@ We benchmark these methods against LibGEOS's `simplify` implementation, which us using BenchmarkTools, Chairmarks, GeoJSON, CairoMakie import GeometryOps as GO, LibGEOS as LG, GeoInterface as GI using CoordinateTransformations +using NaturalEarth import Main: plot_trials # hide lg_and_go(geometry) = (GI.convert(LG, geometry), GO.tuples(geometry)) # Load in the Natural Earth admin GeoJSON, then extract the USA's geometry -fc = GeoJSON.read(read(download("https://rawcdn.githack.com/nvkelso/natural-earth-vector/ca96624a56bd078437bca8184e78163e5039ad19/geojson/ne_10m_admin_0_countries.geojson"))) +fc = NaturalEarth.naturalearth("admin_0_countries", 10) usa_multipoly = fc.geometry[findfirst(==("United States of America"), fc.NAME)] |> x -> GI.convert(LG, x) |> LG.makeValid |> GO.tuples include(joinpath(dirname(dirname(pathof(GO))), "test", "data", "polygon_generation.jl")) @@ -134,7 +135,7 @@ or nested vectors or a table of these. [`RadialDistance`](@ref), [`DouglasPeucker`](@ref), or [`VisvalingamWhyatt`](@ref) algorithms are available, -listed in order of increasing quality but decreaseing performance. +listed in order of increasing quality but decreasing performance. `PoinTrait` and `MultiPointTrait` are returned unchanged. @@ -199,7 +200,7 @@ simplify( ) = _simplify(DouglasPeucker(; kw...), data; prefilter_alg, calc_extent, threaded, crs) -#= For each algorithm, apply simplication to all curves, multipoints, and +#= For each algorithm, apply simplification to all curves, multipoints, and points, reconstructing everything else around them. =# function _simplify(alg::Union{SimplifyAlg, GEOS}, data; prefilter_alg=nothing, kw...) simplifier(geom) = _simplify(GI.trait(geom), alg, geom; prefilter_alg) @@ -245,7 +246,7 @@ Simplifies geometries by removing points less than $SIMPLIFY_ALG_KEYWORDS - `tol`: the minimum distance between points. -Note: user input `tol` is squared to avoid uneccesary computation in algorithm. +Note: user input `tol` is squared to avoid unnecessary computation in algorithm. """ @kwdef struct RadialDistance <: SimplifyAlg number::Union{Int64,Nothing} = nothing @@ -284,7 +285,7 @@ Simplifies geometries by removing points below `tol` distance from the line between its neighboring points. $DOUGLAS_PEUCKER_KEYWORDS -Note: user input `tol` is squared to avoid uneccesary computation in algorithm. +Note: user input `tol` is squared to avoid unnecessary computation in algorithm. """ @kwdef struct DouglasPeucker <: SimplifyAlg number::Union{Int64,Nothing} = nothing @@ -304,7 +305,7 @@ end function _simplify(alg::DouglasPeucker, points::Vector, preserve_endpoint) npoints = length(points) npoints <= MIN_POINTS && return points - # Determine stopping critetia + # Determine stopping criteria max_points = if !isnothing(alg.tol) npoints else @@ -395,7 +396,7 @@ function _simplify(alg::DouglasPeucker, points::Vector, preserve_endpoint) end #= find maximum distance of any point between the start_idx and end_idx to the line formed -by conencting the points at start_idx and end_idx. Note that the first index of maximum +by connecting the points at start_idx and end_idx. Note that the first index of maximum value will be used, which might cause differences in results from other algorithms.=# function _find_max_squared_dist(points, start_idx, end_idx) max_idx = start_idx @@ -422,7 +423,7 @@ distance from the line between its neighboring points. $SIMPLIFY_ALG_KEYWORDS - `tol`: the minimum area of a triangle made with a point and its neighboring points. -Note: user input `tol` is doubled to avoid uneccesary computation in algorithm. +Note: user input `tol` is doubled to avoid unnecessary computation in algorithm. """ @kwdef struct VisvalingamWhyatt <: SimplifyAlg number::Union{Int,Nothing} = nothing diff --git a/src/transformations/transform.jl b/src/transformations/transform.jl index b09e85482..ed06b44fa 100644 --- a/src/transformations/transform.jl +++ b/src/transformations/transform.jl @@ -31,7 +31,7 @@ re.SVector{2, Float64}[[4.5, 3.5], [6.5, 5.5], [8.5, 7.5], [4.5, 3.5]], nothing, rraysCore.SVector{2, Float64}[[6.5, 5.5], [8.5, 7.5], [9.5, 8.5], [6.5, 5.5]], nothing, nothing)], nothing, nothing) ``` -With Rotations.jl you need to actuall multiply the Rotation +With Rotations.jl you need to actually multiply the Rotation by the `SVector` point, which is easy using an anonymous function. ```julia diff --git a/src/utils.jl b/src/utils.jl index 82a1e664d..bf815db2e 100644 --- a/src/utils.jl +++ b/src/utils.jl @@ -126,3 +126,6 @@ function _point_in_extent(p, extent::Extents.Extent) (x1, x2), (y1, y2) = extent.X, extent.Y return x1 ≤ GI.x(p) ≤ x2 && y1 ≤ GI.y(p) ≤ y2 end + +_linearring(geom::GI.LineString) = GI.LinearRing(parent(geom); extent=geom.extent, crs=geom.crs) +_linearring(geom::GI.LinearRing) = geom diff --git a/test/data/polygon_generation.jl b/test/data/polygon_generation.jl index 6fdf3bc74..1242f9bea 100644 --- a/test/data/polygon_generation.jl +++ b/test/data/polygon_generation.jl @@ -24,7 +24,7 @@ Inputs: Output: Vector{Vector{Vector{T}}} representing polygon coordinates Note: - Check your outputs! No guarentee that the polygon's aren't self-intersecting + Check your outputs! No guarantee that the polygon's aren't self-intersecting """ function generate_random_poly( x, diff --git a/test/extensions/flexijoins.jl b/test/extensions/flexijoins.jl index 80e0432ac..0d6f4b287 100644 --- a/test/extensions/flexijoins.jl +++ b/test/extensions/flexijoins.jl @@ -1,22 +1,23 @@ -import GeometryOps as GO, GeoInterface as GI -using FlexiJoins, DataFrames +using Test +using FlexiJoins +using DataFrames +import GeometryOps as GO +import GeoInterface as GI +using ..TestHelpers + +points = GI.MultiPoint(tuple.(rand(100), rand(100))) pl = GI.Polygon([GI.LinearRing([(0, 0), (1, 0), (1, 1), (0, 0)])]) pu = GI.Polygon([GI.LinearRing([(0, 0), (0, 1), (1, 1), (0, 0)])]) -poly_df = DataFrame(geometry = [pl, pu], color = [:red, :blue]) - -points = tuple.(rand(100), rand(100)) -points_df = DataFrame(geometry = points) -@testset "Basic integration" begin - - @test_nowarn joined_df = FlexiJoins.innerjoin((poly_df, points_df), by_pred(:geometry, GO.contains, :geometry)) +@testset_implementations "Polygon DataDrame" begin + points_df = DataFrame(geometry=collect(GI.getpoint($points))) + poly_df = DataFrame(geometry=[$pl, $pu], color=[:red, :blue]) # Test that the join happened correctly joined_df = FlexiJoins.innerjoin((poly_df, points_df), by_pred(:geometry, GO.contains, :geometry)) - @test all(GO.contains.((pl,), joined_df.geometry_1[joined_df.color .== :red])) - @test all(GO.contains.((pu,), joined_df.geometry_1[joined_df.color .== :blue])) + @test all(GO.contains.(($pl,), joined_df.geometry_1[joined_df.color .== :red])) + @test all(GO.contains.(($pu,), joined_df.geometry_1[joined_df.color .== :blue])) # Test that within also works @test_nowarn joined_df = FlexiJoins.innerjoin((points_df, poly_df), by_pred(:geometry, GO.within, :geometry)) - end diff --git a/test/helpers.jl b/test/helpers.jl new file mode 100644 index 000000000..abe6ff8fc --- /dev/null +++ b/test/helpers.jl @@ -0,0 +1,126 @@ +module TestHelpers + +using Test, GeoInterface, ArchGDAL, GeometryBasics, LibGEOS + +export @test_implementations, @testset_implementations + +const TEST_MODULES = [GeoInterface, ArchGDAL, GeometryBasics, LibGEOS] + +# Monkey-patch GeometryBasics to have correct methods. +# TODO: push this up to GB! + +@eval GeometryBasics begin + # MultiGeometry ncoord implementations + GeoInterface.ncoord(::GeoInterface.MultiPolygonTrait, ::GeometryBasics.MultiPolygon{N}) where N = N + GeoInterface.ncoord(::GeoInterface.MultiLineStringTrait, ::GeometryBasics.MultiLineString{N}) where N = N + GeoInterface.ncoord(::GeoInterface.MultiPointTrait, ::GeometryBasics.MultiPoint{N}) where N = N + # LinearRing and LineString confusion + GeometryBasics.geointerface_geomtype(::GeoInterface.LinearRingTrait) = LineString + function GeoInterface.convert(::Type{LineString}, ::GeoInterface.LinearRingTrait, geom) + return GeoInterface.convert(LineString, GeoInterface.LineStringTrait(), geom) # forward to the linestring conversion method + end + # Line interface + GeometryBasics.geointerface_geomtype(::GeoInterface.LineTrait) = Line + function GeoInterface.convert(::Type{Line}, ::GeoInterface.LineTrait, geom) + p1, p2 = GeoInterface.getpoint(geom) + return Line(GeoInterface.convert(Point, GeoInterface.PointTrait(), p1), GeoInterface.convert(Point, GeoInterface.PointTrait(), p2)) + end + # GeometryCollection interface - currently just a large Union + const _ALL_GB_GEOM_TYPES = Union{Point, Line, LineString, Polygon, MultiPolygon, MultiLineString, MultiPoint} + GeometryBasics.geointerface_geomtype(::GeoInterface.GeometryCollectionTrait) = Vector{_ALL_GB_GEOM_TYPES} + function GeoInterface.convert(::Type{Vector{_ALL_GB_GEOM_TYPES}}, ::GeoInterface.GeometryCollectionTrait, geoms) + return _ALL_GB_GEOM_TYPES[GeoInterface.convert(GeometryBasics, g) for g in GeoInterface.getgeom(geoms)] + end +end + + +# Macro to run a block of `code` for multiple modules, +# using GeoInterface.convert for each var in `args` +macro test_implementations(code::Expr) + _test_implementations_inner(TEST_MODULES, code) +end +macro test_implementations(modules::Union{Expr,Vector}, code::Expr) + _test_implementations_inner(modules, code) +end + +function _test_implementations_inner(modules::Union{Expr,Vector}, code::Expr) + vars = Dict{Symbol,Symbol}() + code1 = _quasiquote!(code, vars) + modules1 = modules isa Expr ? modules.args : modules + tests = Expr(:block) + + for mod in modules1 + expr = Expr(:block) + for (var, genkey) in pairs(vars) + push!(expr.args, :($genkey = $GeoInterface.convert($mod, $var))) + end + push!(expr.args, :(@test $code1)) + push!(tests.args, expr) + end + + return esc(tests) +end + +# Macro to run a block of `code` for multiple modules, +# using GeoInterface.convert for each var in `args` +macro testset_implementations(code::Expr) + _testset_implementations_inner("", TEST_MODULES, code) +end +macro testset_implementations(arg, code::Expr) + if arg isa String || arg isa Expr && arg.head == :string + _testset_implementations_inner(arg, TEST_MODULES, code) + else + _testset_implementations_inner("", arg, code) + end +end +macro testset_implementations(title, modules::Union{Expr,Vector}, code::Expr) + _testset_implementations_inner(title, modules, code) +end + +function _testset_implementations_inner(title, modules::Union{Expr,Vector}, code::Expr) + vars = Dict{Symbol,Symbol}() + code1 = _quasiquote!(code, vars) + modules1 = modules isa Expr ? modules.args : modules + testsets = Expr(:block) + + for mod in modules1 + expr = Expr(:block) + for (var, genkey) in pairs(vars) + push!(expr.args, :($genkey = $GeoInterface.convert($mod, $var))) + end + # Manually define the testset macrocall and all string interpolation + testset = Expr( + :macrocall, + Symbol("@testset"), + LineNumberNode(@__LINE__, @__FILE__), + Expr(:string, mod, " ", title), + code1 + ) + push!(expr.args, testset) + push!(testsets.args, expr) + end + + return esc(testsets) +end + +# Taken from BenchmarkTools.jl +_quasiquote!(ex, vars) = ex +function _quasiquote!(ex::Expr, vars::Dict) + if ex.head === :($) + v = ex.args[1] + gen = if v isa Symbol + haskey(vars, v) ? vars[v] : gensym(v) + else + gensym() + end + vars[v] = gen + return v + elseif ex.head !== :quote + for i in 1:length(ex.args) + ex.args[i] = _quasiquote!(ex.args[i], vars) + end + end + return ex +end + +end diff --git a/test/methods/angles.jl b/test/methods/angles.jl index 0f6c55160..4004642a5 100644 --- a/test/methods/angles.jl +++ b/test/methods/angles.jl @@ -1,5 +1,9 @@ using Test -import GeoInterface as GI, GeometryOps as GO, LibGEOS as LG +import GeoInterface as GI +import GeometryOps as GO +import GeometryBasics as GB +import LibGEOS as LG +using ..TestHelpers pt1 = GI.Point((0.0, 0.0)) mpt1 = GI.MultiPoint([pt1, pt1]) @@ -9,39 +13,43 @@ concave_coords = [(0.0, 0.0), (0.0, 1.0), (-1.0, 1.0), (-1.0, 2.0), (2.0, 2.0), l2 = GI.LineString(concave_coords) l3 = GI.LineString(concave_coords[1:(end - 1)]) r1 = GI.LinearRing(concave_coords) -r2 = GI.LinearRing(concave_coords[1:(end - 1)]) -r3 = GI.LinearRing([(1.0, 1.0), (1.0, 1.5), (1.5, 1.5), (1.5, 1.0), (1.0, 1.0)]) +r2 = GI.LinearRing([(1.0, 1.0), (1.0, 1.5), (1.5, 1.5), (1.5, 1.0), (1.0, 1.0)]) concave_angles = [90.0, 270.0, 90.0, 90.0, 90.0, 90.0] -p1 = GI.Polygon([r3]) +p1 = GI.Polygon([r2]) p2 = GI.Polygon([[(0.0, 0.0), (0.0, 4.0), (3.0, 0.0), (0.0, 0.0)]]) p3 = GI.Polygon([[(-3.0, -2.0), (0.0,0.0), (5.0, 0.0), (-3.0, -2.0)]]) p4 = GI.Polygon([r1]) -p5 = GI.Polygon([r1, r3]) +p5 = GI.Polygon([r1, r2]) mp1 = GI.MultiPolygon([p2, p3]) c1 = GI.GeometryCollection([pt1, l2, p2]) -# Points and lines -@test isempty(GO.angles(pt1)) -@test isempty(GO.angles(mpt1)) -@test isempty(GO.angles(l1)) - -# LineStrings and Linear Rings -@test all(isapprox.(GO.angles(l2), concave_angles, atol = 1e-3)) -@test all(isapprox.(GO.angles(l3), concave_angles[2:(end - 1)], atol = 1e-3)) -@test all(isapprox.(GO.angles(r1), concave_angles, atol = 1e-3)) -@test all(isapprox.(GO.angles(r2), concave_angles, atol = 1e-3)) - -# Polygons -p2_angles = [90.0, 36.8699, 53.1301] -p3_angles = [19.6538, 146.3099, 14.0362] -@test all(isapprox.(GO.angles(p1), [90.0 for _ in 1:4], atol = 1e-3)) -@test all(isapprox.(GO.angles(p2), p2_angles, atol = 1e-3)) -@test all(isapprox.(GO.angles(p3), p3_angles, atol = 1e-3)) -@test all(isapprox.(GO.angles(p4), concave_angles, atol = 1e-3)) -@test all(isapprox.(GO.angles(p5), vcat(concave_angles, [270.0 for _ in 1:4]), atol = 1e-3)) - -# Multi-geometries -@test all(isapprox.(GO.angles(mp1), [p2_angles; p3_angles], atol = 1e-3)) -@test all(isapprox.(GO.angles(c1), [concave_angles; p2_angles], atol = 1e-3)) +# Line is not a widely available geometry type +@testset_implementations "line angles" [GB, GI] begin + @test isempty(GO.angles($l1)) +end + +@testset_implementations "angles" begin + # Points and lines + @test isempty(GO.angles($pt1)) + @test isempty(GO.angles($mpt1)) + + # LineStrings and Linear Rings + @test all(isapprox.(GO.angles($l2), concave_angles, atol = 1e-3)) + @test all(isapprox.(GO.angles($l3), concave_angles[2:(end - 1)], atol = 1e-3)) + @test all(isapprox.(GO.angles($r1), concave_angles, atol = 1e-3)) + + # Polygons + p2_angles = [90.0, 36.8699, 53.1301] + p3_angles = [19.6538, 146.3099, 14.0362] + @test all(isapprox.(GO.angles($p1), [90.0 for _ in 1:4], atol = 1e-3)) + @test all(isapprox.(GO.angles($p2), p2_angles, atol = 1e-3)) + @test all(isapprox.(GO.angles($p3), p3_angles, atol = 1e-3)) + @test all(isapprox.(GO.angles($p4), concave_angles, atol = 1e-3)) + @test all(isapprox.(GO.angles($p5), vcat(concave_angles, [270.0 for _ in 1:4]), atol = 1e-3)) + + # Multi-geometries + @test all(isapprox.(GO.angles($mp1), [p2_angles; p3_angles], atol = 1e-3)) + @test all(isapprox.(GO.angles(c1), [concave_angles; p2_angles], atol = 1e-3)) +end diff --git a/test/methods/area.jl b/test/methods/area.jl index c1b4a495e..3b2c7ea01 100644 --- a/test/methods/area.jl +++ b/test/methods/area.jl @@ -1,5 +1,8 @@ using Test -import GeoInterface as GI, GeometryOps as GO, LibGEOS as LG +import GeoInterface as GI +import GeometryOps as GO +import LibGEOS as LG +using ..TestHelpers pt = LG.Point([0.0, 0.0]) empty_pt = LG.readgeom("POINT EMPTY") @@ -33,55 +36,64 @@ p4 = LG.Polygon([ 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]) +c = LG.GeometryCollection([p1, p2, r1]) +c_with_epty_l = 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 +@testset_implementations "That handle empty geoms" [LG] begin + @test GO.area($empty_pt) == LG.area($empty_pt) == 0 + @test GO.area($empty_mpt) == LG.area($empty_mpt) == 0 + @test GO.area($empty_l) == LG.area($empty_l) == 0 + @test GO.area($empty_ml) == LG.area($empty_ml) == 0 + @test GO.area($empty_r) == LG.area($empty_r) == 0 + # Empty polygon + @test GO.signed_area($empty_p) == 0 + @test GO.area($empty_p) == LG.area($empty_p) == 0 + # Empty multipolygon + @test GO.area($empty_mp) == LG.area($empty_mp) == 0 + # Empty collection + @test GO.area(c_with_epty_l) == LG.area(c_with_epty_l) + @test GO.area(c_with_epty_l, Float32) isa Float32 + @test GO.area(empty_c) == LG.area(empty_c) == 0 +end -# 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 -@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 -# Empty polygon -@test GO.area(empty_p) == LG.area(empty_p) == 0 -@test GO.signed_area(empty_p) == 0 +@testset "With GeometryCollection" begin + # Geometry collection summed area + @test GO.area(c) == LG.area(c) + @test GO.area(c, Float32) isa Float32 +end -# 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 +@testset_implementations "all" begin + # Points, lines, and rings have zero area + @test GO.area($pt) == GO.signed_area($pt) == LG.area($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($l1) == GO.signed_area($l1) == LG.area($l1) == 0 + @test GO.area($ml1) == GO.signed_area($ml1) == LG.area($ml1) == 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 + # 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 + @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 -# 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 \ No newline at end of file + # Multipolygon calculations work + @test GO.area($mp1) == a2 + a4 + @test GO.area($mp1, Float32) isa Float32 +end diff --git a/test/methods/barycentric.jl b/test/methods/barycentric.jl index 75319a88d..a6faa9bd9 100644 --- a/test/methods/barycentric.jl +++ b/test/methods/barycentric.jl @@ -1,30 +1,34 @@ using Test import GeoInterface as GI, GeometryOps as GO, LibGEOS as LG -using GeometryOps, GeoInterface, GeometryBasics +using GeometryOps, GeoInterface, GeometryBasics, barycentric_coordinates, MeanValue -@testset "Mean value coordinates" begin - @testset "Preserving return type" begin - @test barycentric_coordinates(MeanValue(), Point2{BigFloat}[(0,0), (1,0), (0,1)], Point2{BigFloat}(1,1)) isa Vector{BigFloat} - @test barycentric_coordinates(MeanValue(), Point2{BigFloat}[(0,0), (1,0), (0,1)], Point2f(1,1)) isa Vector{BigFloat} # keep the most precise type - @test barycentric_coordinates(MeanValue(), Point2{Float64}[(0,0), (1,0), (0,1)], Point2{Float64}(1,1)) isa Vector{Float64} - @test barycentric_coordinates(MeanValue(), Point2{Float32}[(0,0), (1,0), (0,1)], Point2{Float32}(1,1)) isa Vector{Float32} - end - @testset "Triangle coordinates" begin - # Test that the barycentric coordinates for (0,0), (1,0), (0,1), and (0.5,0.5) are (0.5,0.25,0.25) - @test all(barycentric_coordinates(MeanValue(), Point2f[(0,0), (1,0), (0,1)], Point2f(0.5,0.5)) .== (0.5,0.25,0.25)) - # Test that the barycentric coordinates for (0,0), (1,0), (0,1), and (1,1) are (-1,1,1) - @test all(barycentric_coordinates(MeanValue(), Point2f[(0,0), (1,0), (0,1)], Point2f(1,1)) .≈ (-1,1,1)) - # Test that calculations with different number types (in this case Float64) are also accurate - @test all(barycentric_coordinates(MeanValue(), Point2{Float64}[(0,0), (1,0), (0,1)], Point2{Float64}(1,1)) .≈ (-1,1,1)) +t1, p1 = GI.LineString([(0.,0.), (1.,0.), (0.,1.)]), (0.5, 0.5) +t2, p2 = GI.LineString([(0.,0.), (1.,0.), (0.,1.)]), (1., 1.) +t3, p3 = GI.LineString(GeometryBasics.Point2d[(0,0), (1,0), (0,1)]), (1., 1.) + +@testset "Triangle coordinates" begin + # Test that the barycentric coordinates for (0,0), (1,0), (0,1), and (0.5,0.5) are (0.5,0.25,0.25) + @test all(barycentric_coordinates(MeanValue(), t1, p1) .== (0.5,0.25,0.25)) + # Test that the barycentric coordinates for (0,0), (1,0), (0,1), and (1,1) are (-1,1,1) + @test all(barycentric_coordinates(MeanValue(), t2, p2) .≈ (-1,1,1)) + # Test that calculations with different number types (in this case Float64) are also accurate + @test all(barycentric_coordinates(MeanValue(), Point2{Float64}, Point2{Float64}(1,1)) .≈ (-1,1,1)) +end + +@testset "Preserving return type" begin + @test barycentric_coordinates(MeanValue(), Point2{BigFloat}[(0,0), (1,0), (0,1)], Point2{BigFloat}(1,1)) isa Vector{BigFloat} + @test barycentric_coordinates(MeanValue(), Point2{BigFloat}[(0,0), (1,0), (0,1)], Point2f(1,1)) isa Vector{BigFloat} # keep the most precise type + @test barycentric_coordinates(MeanValue(), Point2{Float64}[(0,0), (1,0), (0,1)], Point2{Float64}(1,1)) isa Vector{Float64} + @test barycentric_coordinates(MeanValue(), Point2{Float32}[(0,0), (1,0), (0,1)], Point2{Float32}(1,1)) isa Vector{Float32} +end + +@testset "Tests for helper methods" begin + @testset "`t_value`" begin + @test GO.t_value(Point2f(0,0), Point2f(1,0), 1, 1) == 0 + @test GO.t_value(Point2f(0, 1), Point2f(1, 0), 1, 2) == -0.5f0 end - @testset "Tests for helper methods" begin - @testset "`t_value`" begin - @test GeometryOps.t_value(Point2f(0,0), Point2f(1,0), 1, 1) == 0 - @test GeometryOps.t_value(Point2f(0, 1), Point2f(1, 0), 1, 2) == -0.5f0 - end - @testset "`_det`" begin - @test GeometryOps._det((1,0), (0,1)) == 1f0 - @test GeometryOps._det(Point2f(1, 2), Point2f(3, 4)) == -2.0f0 - end + @testset "`_det`" begin + @test GO._det((1,0), (0,1)) == 1f0 + @test GO._det(Point2f(1, 2), Point2f(3, 4)) == -2.0f0 end end diff --git a/test/methods/centroid.jl b/test/methods/centroid.jl index 431f4732b..224a96da5 100644 --- a/test/methods/centroid.jl +++ b/test/methods/centroid.jl @@ -1,116 +1,125 @@ using Test -import GeoInterface as GI, - GeometryOps as GO, - LibGEOS as LG, - ArchGDAL as AG +import GeoInterface as GI +import GeometryOps as GO +import LibGEOS as LG +import ArchGDAL as AG +using ..TestHelpers @testset "Lines/Rings" begin - # Basic line string - l1 = LG.LineString([[0.0, 0.0], [10.0, 0.0], [10.0, 10.0]]) - c1, len1 = GO.centroid_and_length(l1) - c1_from_LG = LG.centroid(l1) - @test c1[1] ≈ GI.x(c1_from_LG) - @test c1[2] ≈ GI.y(c1_from_LG) - @test len1 ≈ 20.0 - # Spiral line string + l1 = LG.LineString([[0.0, 0.0], [10.0, 0.0], [10.0, 10.0]]) l2 = LG.LineString([[0.0, 0.0], [2.5, -2.5], [-5.0, -3.0], [-4.0, 6.0], [10.0, 10.0], [12.0, -14.56]]) - c2, len2 = GO.centroid_and_length(l2) - c2_from_LG = LG.centroid(l2) - @test c2[1] ≈ GI.x(c2_from_LG) - @test c2[2] ≈ GI.y(c2_from_LG) - @test len2 ≈ 59.3090856788928 - - # Test that non-closed line strings throw an error for centroid_and_area - @test_throws AssertionError GO.centroid_and_area(l2) - # Basic linear ring - note that this still uses weighting by length r1 = LG.LinearRing([[0.0, 0.0], [3456.0, 7894.0], [6291.0, 1954.0], [0.0, 0.0]]) - c3 = GO.centroid(r1) - c3_from_LG = LG.centroid(r1) - @test c3[1] ≈ GI.x(c3_from_LG) - @test c3[2] ≈ GI.y(c3_from_LG) + + @testset_implementations "Basic line string" begin + c1_from_LG = LG.centroid($l1) + c1, len1 = GO.centroid_and_length($l1) + @test c1[1] ≈ GI.x(c1_from_LG) + @test c1[2] ≈ GI.y(c1_from_LG) + @test len1 ≈ 20.0 + end + + @testset_implementations "Spiral line string" begin + c2_from_LG = LG.centroid($l2) + c2, len2 = GO.centroid_and_length($l2) + @test c2[1] ≈ GI.x(c2_from_LG) + @test c2[2] ≈ GI.y(c2_from_LG) + @test len2 ≈ 59.3090856788928 + + # Test that non-closed line strings throw an error for centroid_and_area + @test_throws AssertionError GO.centroid_and_area(l2) + end + + @testset_implementations "Basic linear ring" begin + c3_from_LG = LG.centroid($r1) + c3 = GO.centroid($r1) + @test c3[1] ≈ GI.x(c3_from_LG) + @test c3[2] ≈ GI.y(c3_from_LG) + end end -@testset "Polygons" begin - # Basic rectangle - p1 = AG.fromWKT("POLYGON((0 0, 10 0, 10 10, 0 10, 0 0))") - c1 = GO.centroid(p1) +# Basic rectangle +p1 = AG.fromWKT("POLYGON((0 0, 10 0, 10 10, 0 10, 0 0))") +# Concave c-shaped polygon +p2 = LG.Polygon([[ + [11.0, 0.0], [11.0, 3.0], [14.0, 3.0], [14.0, 2.0], [12.0, 2.0], + [12.0, 1.0], [14.0, 1.0], [14.0, 0.0], [11.0, 0.0], +]]) +# Randomly generated polygon with lots of sides +p3 = LG.Polygon([[ + [14.567, 8.974], [13.579, 8.849], [12.076, 8.769], [11.725, 8.567], + [11.424, 6.451], [10.187, 7.712], [8.187, 6.795], [8.065, 8.096], + [6.827, 8.287], [7.1628, 8.9221], [5.8428, 10.468], [7.987, 11.734], + [7.989, 12.081], [8.787, 11.930], [7.568, 13.926], [9.330, 13.340], + [9.6817, 13.913], [10.391, 12.222], [12.150, 12.032], [14.567, 8.974], +]]) +# Polygon with one hole +p4 = LG.Polygon([ + [[0.0, 0.0], [10.0, 0.0], [10.0, 10.0], [0.0, 10.0], [0.0, 0.0]], + [[2.3, 2.7], [2.5, 4.9], [4.1, 5.2], [4.2, 1.9], [2.3, 2.7]], +]) +# Polygon with two holes +p5 = LG.Polygon([ + [[-10.0, -10.0], [-2.0, 0.0], [6.0, -10.0], [-10.0, -10.0]], + [[-8.0, -8.0], [-8.0, -7.0], [-4.0, -7.0], [-4.0, -8.0], [-8.0, -8.0]], + [[-3.0,-9.0], [3.0, -9.0], [3.0, -8.5], [-3.0, -8.5], [-3.0, -9.0]], +]) +# Same polygon as P5 but using a GeoInterface polygon +p6 = GI.Polygon([ + [(-10.0, -10.0), (-2.0, 0.0), (6.0, -10.0), (-10.0, -10.0)], + [(-8.0, -8.0), (-8.0, -7.0), (-4.0, -7.0), (-4.0, -8.0), (-8.0, -8.0)], + [(-3.0, -9.0), (3.0, -9.0), (3.0, -8.5), (-3.0, -8.5), (-3.0, -9.0)], +]) + +@testset_implementations "Polygons" begin + c1 = GO.centroid($p1) c1 .≈ (5, 5) @test GI.x(c1) ≈ 5 @test GI.y(c1) ≈ 5 - # Concave c-shaped polygon - p2 = LG.Polygon([[ - [11.0, 0.0], [11.0, 3.0], [14.0, 3.0], [14.0, 2.0], [12.0, 2.0], - [12.0, 1.0], [14.0, 1.0], [14.0, 0.0], [11.0, 0.0], - ]]) - c2, area2 = GO.centroid_and_area(p2) - c2_from_LG = LG.centroid(p2) + c2, area2 = GO.centroid_and_area($p2) + c2_from_LG = LG.centroid($p2) @test c2[1] ≈ GI.x(c2_from_LG) @test c2[2] ≈ GI.y(c2_from_LG) @test area2 ≈ LG.area(p2) - # Randomly generated polygon with lots of sides - p3 = LG.Polygon([[ - [14.567, 8.974], [13.579, 8.849], [12.076, 8.769], [11.725, 8.567], - [11.424, 6.451], [10.187, 7.712], [8.187, 6.795], [8.065, 8.096], - [6.827, 8.287], [7.1628, 8.9221], [5.8428, 10.468], [7.987, 11.734], - [7.989, 12.081], [8.787, 11.930], [7.568, 13.926], [9.330, 13.340], - [9.6817, 13.913], [10.391, 12.222], [12.150, 12.032], [14.567, 8.974], - ]]) - c3, area3 = GO.centroid_and_area(p3) - c3_from_LG = LG.centroid(p3) + c3, area3 = GO.centroid_and_area($p3) + c3_from_LG = LG.centroid($p3) @test c3[1] ≈ GI.x(c3_from_LG) @test c3[2] ≈ GI.y(c3_from_LG) - @test area3 ≈ LG.area(p3) + @test area3 ≈ LG.area($p3) - # Polygon with one hole - p4 = LG.Polygon([ - [[0.0, 0.0], [10.0, 0.0], [10.0, 10.0], [0.0, 10.0], [0.0, 0.0]], - [[2.3, 2.7], [2.5, 4.9], [4.1, 5.2], [4.2, 1.9], [2.3, 2.7]], - ]) - c4, area4 = GO.centroid_and_area(p4) - c4_from_LG = LG.centroid(p4) + c4, area4 = GO.centroid_and_area($p4) + c4_from_LG = LG.centroid($p4) @test c4[1] ≈ GI.x(c4_from_LG) @test c4[2] ≈ GI.y(c4_from_LG) - @test area4 ≈ LG.area(p4) + @test area4 ≈ LG.area($p4) - # Polygon with two holes - p5 = LG.Polygon([ - [[-10.0, -10.0], [-2.0, 0.0], [6.0, -10.0], [-10.0, -10.0]], - [[-8.0, -8.0], [-8.0, -7.0], [-4.0, -7.0], [-4.0, -8.0], [-8.0, -8.0]], - [[-3.0,-9.0], [3.0, -9.0], [3.0, -8.5], [-3.0, -8.5], [-3.0, -9.0]], - ]) - c5 = GO.centroid(p5) - c5_from_LG = LG.centroid(p5) + c5 = GO.centroid($p5) + c5_from_LG = LG.centroid($p5) @test c5[1] ≈ GI.x(c5_from_LG) @test c5[2] ≈ GI.y(c5_from_LG) - - # Same polygon as P5 but using a GeoInterface polygon - p6 = GI.Polygon([ - [(-10.0, -10.0), (-2.0, 0.0), (6.0, -10.0), (-10.0, -10.0)], - [(-8.0, -8.0), (-8.0, -7.0), (-4.0, -7.0), (-4.0, -8.0), (-8.0, -8.0)], - [(-3.0, -9.0), (3.0, -9.0), (3.0, -8.5), (-3.0, -8.5), (-3.0, -9.0)], - ]) - c6 = GO.centroid(p6) + c6 = GO.centroid($p6) @test all(c5 .≈ c6) end -@testset "MultiPolygons" begin - # Combine poylgons made above - m1 = LG.MultiPolygon([ - [ - [[11.0, 0.0], [11.0, 3.0], [14.0, 3.0], [14.0, 2.0], [12.0, 2.0], - [12.0, 1.0], [14.0, 1.0], [14.0, 0.0], [11.0, 0.0]], - ], - [ - [[0.0, 0.0], [10.0, 0.0], [10.0, 10.0], [0.0, 10.0], [0.0, 0.0]], - [[2.3, 2.7], [2.5, 4.9], [4.1, 5.2], [4.2, 1.9], [2.3, 2.7]], - ] - ]) - c1, area1 = GO.centroid_and_area(m1) - c1_from_LG = LG.centroid(m1) + +# Combine polygons made above +m1 = LG.MultiPolygon([ + [ + [[11.0, 0.0], [11.0, 3.0], [14.0, 3.0], [14.0, 2.0], [12.0, 2.0], + [12.0, 1.0], [14.0, 1.0], [14.0, 0.0], [11.0, 0.0]], + ], + [ + [[0.0, 0.0], [10.0, 0.0], [10.0, 10.0], [0.0, 10.0], [0.0, 0.0]], + [[2.3, 2.7], [2.5, 4.9], [4.1, 5.2], [4.2, 1.9], [2.3, 2.7]], + ] +]) + +@testset_implementations "MultiPolygons" begin + c1, area1 = GO.centroid_and_area($m1) + c1_from_LG = LG.centroid($m1) @test c1[1] ≈ GI.x(c1_from_LG) @test c1[2] ≈ GI.y(c1_from_LG) - @test area1 ≈ LG.area(m1) + @test area1 ≈ LG.area($m1) end diff --git a/test/methods/clipping/coverage.jl b/test/methods/clipping/coverage.jl index 01189687a..7127bb4c6 100644 --- a/test/methods/clipping/coverage.jl +++ b/test/methods/clipping/coverage.jl @@ -1,64 +1,80 @@ using Test -import GeoInterface as GI, GeometryOps as GO, LibGEOS as LG -using GeoInterface -using GeoInterface.Extents +import GeoInterface as GI +import GeometryBasics as GB +import GeometryOps as GO +import LibGEOS as LG +using ..TestHelpers cell_extremes = (0.0, 20.0, 0.0, 20.0) -cell_extent = Extents.Extent(X = (0.0, 20.0), Y = (0.0, 20.0)) +cell_extent = GI.Extents.Extent(X = (0.0, 20.0), Y = (0.0, 20.0)) cell_poly = GI.Polygon([[(0.0, 0.0), (0.0, 20.0), (20.0, 20.0), (20.0, 0.0), (0.0, 0.0)]]) cell_area = 400.0 -# # Points, lines, curves pt1 = (1.0, 0.0) -@test GO.coverage(pt1, cell_extremes...) == 0.0 mpt1 = GI.MultiPoint([pt1, (25.0, 0.0), (100.0, 50.0)]) -@test GO.coverage(mpt1, cell_extremes...) == 0.0 -l1 = GI.LineString([(-1.0, 5.0), (5.0, 10.0), (10.0, -2.0), (22, 5.0)]) -@test GO.coverage(l1, cell_extremes...) == 0.0 +l1 = GI.LineString([(-1.0, 5.0), (5.0, 10.0), (10.0, -2.0), (22.0, 5.0)]) + +@testset_implementations "Points, lines, curves" begin + @test GO.coverage($pt1, cell_extremes...) == 0.0 + @test GO.coverage($mpt1, cell_extremes...) == 0.0 + @test GO.coverage($l1, cell_extremes...) == 0.0 +end # # Polygons -# polygon is the same as the cell (input is extent) p1 = GI.Polygon([[(0.0, 0.0), (0.0, 20.0), (20.0, 20.0), (20.0, 0.0), (0.0, 0.0)]]) -@test GO.coverage(p1, cell_extent) == cell_area -# polygon is the same as the cell (input is min/max values) -@test GO.coverage(p1, cell_extremes...) == cell_area -# polygon is bigger than the cell -p2 = GI.Polygon([[(-10, -10.0), (-10.0, 30.0), (30.0, 30.0), (300.0, -10.0), (-10.0, -10.0)]]) -@test GO.coverage(p2, cell_extremes...) == cell_area -# polygon is completly inside of cell +p2 = GI.Polygon([[(-10.0, -10.0), (-10.0, 30.0), (30.0, 30.0), (300.0, -10.0), (-10.0, -10.0)]]) +p2b = GI.Polygon([[(-10, -10.0), (-10.0, 30.0), (30.0, 30.0), (300.0, -10.0), (-10.0, -10.0)]]) +# polygon is completely inside of cell p3 = GI.Polygon([[(5.0, 5.0), (5.0, 15.0), (15.0, 15.0), (15.0, 5.0), (5.0, 5.0)]]) -@test GO.coverage(p3, cell_extremes...) ≈ LG.area(LG.intersection(cell_poly, p3)) -# polygon exits cell through one edge p4 = GI.Polygon([[(5.0, 5.0), (5.0, 25.0), (15.0, 25.0), (15.0, 5.0), (5.0, 5.0)]]) -@test GO.coverage(p4, cell_extremes...) ≈ LG.area(LG.intersection(cell_poly, p4)) p5 = GI.Polygon([[(5.0, 5.0), (5.0, 25.0), (25.0, 25.0), (25.0, 5.0), (5.0, 5.0)]]) -@test GO.coverage(p5, cell_extremes...) ≈ LG.area(LG.intersection(cell_poly, p5)) -# polygon exits cell through multiple edges (north and east) p6 = GI.Polygon([[(20.8826, 6.4239), (15.9663, 2.3014), (8.6078, 2.0995), (2.6849, 6.4088), (0.8449, 12.7452), (3.0813, 19.1654), (9.1906, 23.2520), (15.5835, 22.9101), (20.9143, 18.5933), (20.8826, 6.4239)]]) -@test GO.coverage(p6, cell_extremes...) ≈ LG.area(LG.intersection(p6, cell_poly)) -# polygon exits cell through multiple edges (west and east) p7 = GI.Polygon([[(-5.0, 10.0), (-5.0, 25.0), (25.0, 25.0), (25.0, 10.0), (-5.0, 10.0)]]) -@test GO.coverage(p7, cell_extremes...) ≈ LG.area(LG.intersection(p7, cell_poly)) -# non-convex polygon split into two pieces (same wall) p8 = GI.Polygon([[(-10.0, 15.0), (10.0, 15.0), (10.0, 12.0), (-5.0, 12.0), (-5.0, 9.0), (10.0, 9.0), (10.0, 6.0), (-10.0, 6.0), (-10.0, 15.0)]]) -@test GO.coverage(p8, cell_extremes...) ≈ LG.area(LG.intersection(p8, cell_poly)) -# counter-clockwise polygon p9 = GI.Polygon([[(-10.0, 15.0), (-10.0, 6.0), (10.0, 6.0), (10.0, 9.0), (-5.0, 9.0), (-5.0, 12.0), (10.0, 12.0), (10.0, 15.0), (-10.0, 15.0),]]) -@test GO.coverage(p9, cell_extremes...) ≈ LG.area(LG.intersection(p9, cell_poly)) -# polygon with a hole p10 = GI.Polygon([[(0.0, 0.0), (0.0, 20.0), (20.0, 20.0), (20.0, 0.0), (0.0, 0.0)], [(10.0, 10.0), (10.0, 15.0), (15.0, 15.0), (15.0, 10.0), (10.0, 10.0)], [(7.0, 7.0), (5.0, 7.0), (5.0, 5.0), (7.0, 7.0)], ]) -@test GO.coverage(p10, cell_extremes...) ≈ LG.area(LG.intersection(p10, cell_poly)) -# non-convex polygon split into two pieces (differnet walls) +# non-convex polygon split into two pieces (different walls) p11 = GI.Polygon([[(-10.0, 15.0), (-10.0, -10.0), (15.0, -10.0), (15.0, 5.0), (10.0, 5.0), (10.0, -5.0), (-5.0, -5.0), (-5.0, 10.0), (5.0, 10.0), (5.0, 15.0), (-10.0, 15.0)]]) -@test GO.coverage(p11, cell_extremes...) ≈ LG.area(LG.intersection(p11, cell_poly)) -# # MultiPolygons +@testset_implementations "Polygons" [GI, GB, LG] begin + # polygon is the same as the cell (input is min/max values) + @test GO.coverage($p1, cell_extent) == cell_area + @test GO.coverage($p1, cell_extremes...) == cell_area + @test GO.coverage($p10, cell_extremes...) ≈ LG.area(LG.intersection($p10, cell_poly)) + # polygon is the same as the cell (input is extent) + @test GO.coverage($p1, cell_extent) == cell_area + # polygon is the same as the cell (input is min/max values) + @test GO.coverage($p1, cell_extremes...) == cell_area + # polygon is bigger than the cell + @test GO.coverage($p2, cell_extremes...) == cell_area + # polygon is bigger than the cell + @test_implementations GO.coverage($p2b, cell_extremes...) == cell_area + # polygon is completly inside of cell + @test GO.coverage($p3, cell_extremes...) ≈ LG.area(LG.intersection(cell_poly, $p3)) + # polygon exits cell through one edge + @test GO.coverage($p4, cell_extremes...) ≈ LG.area(LG.intersection(cell_poly, $p4)) + @test GO.coverage($p5, cell_extremes...) ≈ LG.area(LG.intersection(cell_poly, $p5)) + # polygon exits cell through multiple edges (north and east) + @test GO.coverage($p6, cell_extremes...) ≈ LG.area(LG.intersection($p6, cell_poly)) + # polygon exits cell through multiple edges (west and east) + @test GO.coverage($p7, cell_extremes...) ≈ LG.area(LG.intersection($p7, cell_poly)) + # non-convex polygon split into two pieces (same wall) + @test GO.coverage($p8, cell_extremes...) ≈ LG.area(LG.intersection($p8, cell_poly)) + # counter-clockwise polygon + @test GO.coverage($p9, cell_extremes...) ≈ LG.area(LG.intersection($p9, cell_poly)) + # polygon with a hole + @test GO.coverage($p10, cell_extremes...) ≈ LG.area(LG.intersection($p10, cell_poly)) + # non-convex polygon split into two pieces (differnet walls) + @test GO.coverage($p11, cell_extremes...) ≈ LG.area(LG.intersection($p11, cell_poly)) +end + +# TODO: MultiPolygons and GeometryCollection mp1 = GI.MultiPolygon([p1, p2]) == 2cell_area -gc1 = GI.GeometryCollection([pt1, l1, p1]) == cell_area \ No newline at end of file +gc1 = GI.GeometryCollection([pt1, l1, p1]) == cell_area diff --git a/test/methods/clipping/cut.jl b/test/methods/clipping/cut.jl index b2e0630e2..7530bde2a 100644 --- a/test/methods/clipping/cut.jl +++ b/test/methods/clipping/cut.jl @@ -1,7 +1,10 @@ using Test +import GeoInterface as GI +import GeometryOps as GO +using ..TestHelpers -import GeoInterface as GI, GeometryOps as GO - +# TODO this should also work with a LineString +# Some geom packages don't have `Line` l1 = GI.Line([(5.0, -5.0), (5.0, 15.0)]) l2 = GI.Line([(-1.0, 8.0), (2.0, 11.0)]) l3 = GI.Line([(-1.0, 6.0), (11.0, 6.0)]) @@ -13,9 +16,8 @@ p1 = GI.Polygon([r1]) p2 = GI.Polygon([r1, r2]) p3 = GI.Polygon([[(0.0, 0.0), (0.0, 5.0), (2.5, 7.5), (5.0, 5.0), (7.5, 7.5), (10.0, 5.0), (10.0, 0.0), (0.0, 0.0)]]) -@testset "Cut Polygon" begin - # Cut convex polygon - cut_polys = GO.cut(p1, l1) +@testset_implementations "Cut convex polygon" begin + cut_polys = GO.cut($p1, l1) @test all(GO.equals.( cut_polys, [ @@ -24,7 +26,7 @@ p3 = GI.Polygon([[(0.0, 0.0), (0.0, 5.0), (2.5, 7.5), (5.0, 5.0), (7.5, 7.5), (1 ], )) - cut_polys = GO.cut(p1, l2) + cut_polys = GO.cut($p1, l2) @test all(GO.equals.( cut_polys, [ @@ -34,7 +36,7 @@ p3 = GI.Polygon([[(0.0, 0.0), (0.0, 5.0), (2.5, 7.5), (5.0, 5.0), (7.5, 7.5), (1 )) # Cut convex polygon with hole - cut_polys = GO.cut(p2, l1) + cut_polys = GO.cut($p2, l1) @test all(GO.equals.( cut_polys, [ @@ -42,7 +44,7 @@ p3 = GI.Polygon([[(0.0, 0.0), (0.0, 5.0), (2.5, 7.5), (5.0, 5.0), (7.5, 7.5), (1 GI.Polygon([[(5.0, 0.0), (10.0, 0.0), (10.0, 10.0), (5.0, 10.0), (5.0, 4.0), (8.0, 4.0), (8.0, 2.0), (5.0, 2.0), (5.0, 0.0)]]), ], )) - cut_polys = GO.cut(p2, l3) + cut_polys = GO.cut($p2, l3) @test all(GO.equals.( cut_polys, [ @@ -52,7 +54,7 @@ p3 = GI.Polygon([[(0.0, 0.0), (0.0, 5.0), (2.5, 7.5), (5.0, 5.0), (7.5, 7.5), (1 )) # Cut concave polygon into three pieces - cut_polys = GO.cut(p3, l3) + cut_polys = GO.cut($p3, l3) @test all(GO.equals.( cut_polys, [ @@ -63,8 +65,8 @@ p3 = GI.Polygon([[(0.0, 0.0), (0.0, 5.0), (2.5, 7.5), (5.0, 5.0), (7.5, 7.5), (1 )) # Line doesn't cut through polygon - cut_polys = GO.cut(p1, l4) - @test all(GO.equals.(cut_polys, [p1])) - cut_polys = GO.cut(p1, l5) - @test all(GO.equals.(cut_polys, [p1])) -end \ No newline at end of file + cut_polys = GO.cut($p1, l4) + @test all(GO.equals.(cut_polys, [$p1])) + cut_polys = GO.cut($p1, l5) + @test all(GO.equals.(cut_polys, [$p1])) +end diff --git a/test/methods/clipping/intersection_points.jl b/test/methods/clipping/intersection_points.jl index 28c5a00f4..4cfc7074a 100644 --- a/test/methods/clipping/intersection_points.jl +++ b/test/methods/clipping/intersection_points.jl @@ -1,5 +1,8 @@ using Test -import GeoInterface as GI, GeometryOps as GO, LibGEOS as LG +import GeoInterface as GI +import GeometryOps as GO +import LibGEOS as LG +using ..TestHelpers l1 = GI.LineString([(90000.0, 1000.0), (90000.0, 22500.0), (95000.0, 22500.0), (95000.0, 1000.0), (90000.0, 1000.0)]) l2 = GI.LineString([(90000.0, 7500.0), (107500.0, 27500.0), (112500.0, 27500.0), (95000.0, 7500.0), (90000.0, 7500.0)]) @@ -10,25 +13,22 @@ l6 = GI.LineString([(0.0, 25000.0), (0.0, 29000.0), (20000.0, 29000.0), (20000.0 p1, p2 = GI.Polygon([l1]), GI.Polygon([l2]) -# Three intersection points -GO_l1_l2_mp = GI.MultiPoint(GO.intersection_points(l1, l2)) -LG_l1_l2_mp = GI.MultiPoint(collect(GI.getpoint(LG.intersection(l1, l2)))) -@test GO.equals(GO_l1_l2_mp, LG_l1_l2_mp) +@testset_implementations begin + # Three intersection points + LG_l1_l2_mp = GI.MultiPoint(collect(GI.getpoint(LG.intersection($l1, $l2)))) + @test GO.equals(GI.MultiPoint(GO.intersection_points($l1, $l2)), LG_l1_l2_mp) -# Four intersection points with large intersection -GO_l3_l4_mp = GI.MultiPoint(GO.intersection_points(l3, l4)) -LG_l3_l4_mp = GI.MultiPoint(collect(GI.getpoint(LG.intersection(l3, l4)))) -@test GO.equals(GO_l3_l4_mp, LG_l3_l4_mp) + # Four intersection points with large intersection + LG_l3_l4_mp = GI.MultiPoint(collect(GI.getpoint(LG.intersection($l3, $l4)))) + @test GO.equals(GI.MultiPoint(GO.intersection_points($l3, $l4)), LG_l3_l4_mp) -# Four intersection points with very small intersection -GO_l5_l6_mp = GI.MultiPoint(GO.intersection_points(l5, l6)) -LG_l5_l6_mp = GI.MultiPoint(collect(GI.getpoint(LG.intersection(l5, l6)))) -@test GO.equals(GO_l5_l6_mp, LG_l5_l6_mp) + # Four intersection points with very small intersection + LG_l5_l6_mp = GI.MultiPoint(collect(GI.getpoint(LG.intersection($l5, $l6)))) + @test GO.equals(GI.MultiPoint(GO.intersection_points($l5, $l6)), LG_l5_l6_mp) -# Test that intersection points between lines and polygons is equivalent -GO_p1_p2_mp = GI.MultiPoint(GO.intersection_points(p1, p2)) -@test GO.equals(GO_p1_p2_mp, GO_l1_l2_mp) + # Test that intersection points between lines and polygons is equivalent + @test GO.equals(GI.MultiPoint(GO.intersection_points($p1, $p2)), GI.MultiPoint(GO.intersection_points($l1, $l2))) -# No intersection points between polygon and line -GO_p1_l6 = GO.intersection_points(p1, l6) -@test isempty(GO_p1_l6) + # No intersection points between polygon and line + @test isempty(GO.intersection_points($p1, $l6)) +end diff --git a/test/methods/clipping/polygon_clipping.jl b/test/methods/clipping/polygon_clipping.jl index 17a441deb..432650e94 100644 --- a/test/methods/clipping/polygon_clipping.jl +++ b/test/methods/clipping/polygon_clipping.jl @@ -1,5 +1,8 @@ using Test -import GeometryOps as GO, GeoInterface as GI, LibGEOS as LG +import GeometryOps as GO +import GeoInterface as GI +import LibGEOS as LG +using ..TestHelpers # Test of polygon clipping p1 = GI.Polygon([[(0.0, 0.0), (5.0, 5.0), (10.0, 0.0), (5.0, -5.0), (0.0, 0.0)]]) @@ -121,12 +124,12 @@ test_pairs = [ (p13, p14, "p13", "p14", "Polygons whose intersection is two distinct regions"), (p15, p16, "p15", "p16", "First polygon in second polygon"), (p16, p15, "p16", "p15", "Second polygon in first polygon"), - (p17, p18, "p17", "p18", "First polygon with a hole (hole completly in second polygon), second without a hole"), - (p18, p17, "p18", "p17", "First polygon with no hole, second with a hole (hole completly in first polygon)"), - (p19, p20, "p19", "p20", "First polygon with a hole (hole not completly in second polygon), second without a hole"), - (p42, p20, "p42", "p20", "First polygon with two holes (holes not completly in second polygon), second without a hole"), - (p20, p19, "p20", "p19", "First polygon with no hole, second with a hole (hole not completly in first polygon)"), - (p20, p42, "p20", "p42", "First polygon with no holes, second with two holes (holes not completly in second polygon)"), + (p17, p18, "p17", "p18", "First polygon with a hole (hole completely in second polygon), second without a hole"), + (p18, p17, "p18", "p17", "First polygon with no hole, second with a hole (hole completely in first polygon)"), + (p19, p20, "p19", "p20", "First polygon with a hole (hole not completely in second polygon), second without a hole"), + (p42, p20, "p42", "p20", "First polygon with two holes (holes not completely in second polygon), second without a hole"), + (p20, p19, "p20", "p19", "First polygon with no hole, second with a hole (hole not completely in first polygon)"), + (p20, p42, "p20", "p42", "First polygon with no holes, second with two holes (holes not completely in second polygon)"), (p21, p22, "p21", "p22", "Polygons form cross, splitting each other"), (p23, p24, "p23", "p24", "Polygons are both donuts with intersecting holes"), (p25, p26, "p25", "p26", "Polygons both have two holes that intersect in various ways"), @@ -136,8 +139,8 @@ test_pairs = [ (p33, p34, "p33", "p34", "One polygon inside of the other, sharing an edge"), (p33, p35, "p33", "p35", "Polygons outside of one another, sharing an edge"), (p33, p36, "p33", "p36", "All intersection points are V-intersections as defined by Foster"), - (p33, p37, "p33", "p37", "Polygons are completly disjoint (no holes)"), - (p38, p39, "p38", "p39", "Polygons are completly disjoint (both have one hole)"), + (p33, p37, "p33", "p37", "Polygons are completely disjoint (no holes)"), + (p38, p39, "p38", "p39", "Polygons are completely disjoint (both have one hole)"), (p40, p41, "p40", "p41", "Two overlapping polygons with three total holes in overlap region"), (p42, p43, "p42", "p43", "First polygon 2 holes, second polygon 3 holes. Holes do not overlap"), (p43, p42, "p43", "p42", "First polygon 3 holes, second polygon 2 holes. Holes do not overlap"), @@ -202,9 +205,9 @@ end # Test clipping functions and print error message if tests fail function test_clipping(GO_f, LG_f, f_name) for (p1, p2, sg1, sg2, sdesc) in test_pairs - pass_test = compare_GO_LG_clipping(GO_f, LG_f, p1, p2) - @test pass_test - !pass_test && println("\n↑ TEST INFO: $sg1 $f_name $sg2 - $sdesc \n\n") + @testset_implementations "$sg1 $f_name $sg2 - $sdesc" begin + @test compare_GO_LG_clipping(GO_f, LG_f, $p1, $p2) + end end return end diff --git a/test/methods/convex_hull.jl b/test/methods/convex_hull.jl new file mode 100644 index 000000000..7caeb0c3d --- /dev/null +++ b/test/methods/convex_hull.jl @@ -0,0 +1,45 @@ +using Test +import GeoInterface as GI +import GeometryOps as GO +import LibGEOS as LG +using ..TestHelpers + +@testset "Basic example" begin + points = tuple.(rand(100), rand(100)) + hull = GO.convex_hull(points) + @test !GO.isconcave(hull) skip=true # TODO: fix + @test !GO.isclockwise(GI.getexterior(hull)) # exterior should be ccw + @test all(x -> GO.covers(hull, x), points) # but the orientation is right + # Test against LibGEOS, by testing that all the points are the same + # This is robust against winding order and starting at a different point, etc. + @test isempty( + setdiff( + collect(GO.flatten(GO.tuples, GI.PointTrait, hull)), + collect(GO.flatten(GO.tuples, GI.PointTrait, GO.convex_hull(GO.GEOS(), points))) + ) + ) +end + +@testset "Duplicated points" begin + points = tuple.(rand(100), rand(100)) + @test_nowarn hull = GO.convex_hull(vcat(points, points)) + single_hull = GO.convex_hull(points) + double_hull = GO.convex_hull(vcat(points, points)) + + @test GO.equals(GI.getexterior(single_hull), GI.getexterior(double_hull)) + @test !GO.isconcave(double_hull) skip=true # TODO: fix + @test !GO.isclockwise(GI.getexterior(double_hull)) # exterior should be ccw + @test all(x -> GO.covers(single_hull, x), points) + @test all(x -> GO.covers(double_hull, x), points) + # Test against LibGEOS, by testing that all the points are the same + # This is robust against winding order and starting at a different point, etc. + @test isempty( + setdiff( + collect(GO.flatten(GO.tuples, GI.PointTrait, double_hull)), + collect(GO.flatten(GO.tuples, GI.PointTrait, GO.convex_hull(GO.GEOS(), points))) + ) + ) +end + +# The rest of the tests are handled in DelaunayTriangulation.jl, this is simply testing +# that the methods work as expected. diff --git a/test/methods/distance.jl b/test/methods/distance.jl index fee0fba9a..a2a48c865 100644 --- a/test/methods/distance.jl +++ b/test/methods/distance.jl @@ -1,5 +1,9 @@ using Test -import GeoInterface as GI, GeometryOps as GO, LibGEOS as LG +import ArchGDAL as AG +import GeoInterface as GI +import GeometryOps as GO +import LibGEOS as LG +using ..TestHelpers pt1 = LG.Point([0.0, 0.0]) pt2 = LG.Point([0.0, 1.0]) @@ -30,85 +34,88 @@ 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 - -#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) -# Return types -@test GO.distance(pt1, l1) isa Float64 -@test GO.distance(pt1, l1, Float32) isa Float32 - -# 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) -# 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 -@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))) - -# Point and Geometry Collection -@test GO.distance(pt1, c1) == LG.distance(pt1, c1) - +@testset_implementations "Where LinearRing exists" [LG, GI] begin + # 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) +end + +@testset_implementations "Where GeometryCollection exists" [LG, AG, GI] begin + @test GO.distance($pt1, c1) == LG.distance($pt1, c1) +end + +@testset_implementations "Point and Point" begin + # 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 +end + +@testset_implementations "Point and Line" begin + #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) + # Return types + @test GO.distance($pt1, $l1) isa Float64 + @test GO.distance($pt1, $l1, Float32) isa Float32 +end + +@testset_implementations "Point and Polygon" begin + # 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 polygon 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) + # Return types + @test GO.distance($pt1, $p1) isa Float64 + @test GO.distance($pt1, $p1, Float32) isa Float32 +end + +@testset_implementations "Point and MultiPoint" begin + @test GO.distance($pt4, $mpt1) == LG.distance(pt4, mpt1) + @test GO.distance($pt4, $mpt1) isa Float64 + @test GO.distance($pt4, $mpt1, Float32) isa Float32 +end + +@testset_implementations "Point and MultiPolygon" begin + # 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))) +end diff --git a/test/methods/equals.jl b/test/methods/equals.jl index 1fd725f36..38fa45c04 100644 --- a/test/methods/equals.jl +++ b/test/methods/equals.jl @@ -1,134 +1,149 @@ using Test -import GeoInterface as GI, GeometryOps as GO, LibGEOS as LG +import GeoInterface as GI +import GeometryBasics as GB +import GeometryOps as GO +import LibGEOS as LG +using ..TestHelpers -@testset "Points/MultiPoints" begin - p1 = LG.Point([0.0, 0.0]) - p2 = LG.Point([0.0, 1.0]) +p1 = LG.Point([0.0, 0.0]) +p2 = LG.Point([0.0, 1.0]) +mp1 = LG.MultiPoint([[0.0, 1.0], [2.0, 2.0]]) +mp2 = LG.MultiPoint([[0.0, 1.0], [2.0, 2.0], [3.0, 3.0]]) +mp3 = LG.MultiPoint([p2]) + +@testset_implementations "Points/MultiPoints" begin # Same points - @test GO.equals(p1, p1) == LG.equals(p1, p1) - @test GO.equals(p2, p2) == LG.equals(p2, p2) + @test GO.equals($p1, $p1) == LG.equals($p1, $p1) + @test GO.equals($p2, $p2) == LG.equals($p2, $p2) # Different points - @test GO.equals(p1, p2) == LG.equals(p1, p2) + @test GO.equals($p1, $p2) == LG.equals($p1, $p2) - mp1 = LG.MultiPoint([[0.0, 1.0], [2.0, 2.0]]) - mp2 = LG.MultiPoint([[0.0, 1.0], [2.0, 2.0], [3.0, 3.0]]) - mp3 = LG.MultiPoint([p2]) # Same points - @test GO.equals(mp1, mp1) == LG.equals(mp1, mp1) - @test GO.equals(mp2, mp2) == LG.equals(mp2, mp2) + @test GO.equals($mp1, $mp1) == LG.equals($mp1, $mp1) + @test GO.equals($mp2, $mp2) == LG.equals($mp2, $mp2) # Different points - @test GO.equals(mp1, mp2) == LG.equals(mp1, mp2) - @test GO.equals(mp1, p1) == LG.equals(mp1, p1) + @test GO.equals($mp1, $mp2) == LG.equals($mp1, $mp2) + @test GO.equals($mp1, $p1) == LG.equals($mp1, $p1) # Point and multipoint - @test GO.equals(p2, mp3) == LG.equals(p2, mp3) + @test GO.equals($p2, $mp3) == LG.equals($p2, $mp3) end -@testset "Lines/Rings" begin - l1 = LG.LineString([[0.0, 0.0], [0.0, 10.0]]) - l2 = LG.LineString([[0.0, -10.0], [0.0, 20.0]]) +l1 = LG.LineString([[0.0, 0.0], [0.0, 10.0]]) +l2 = LG.LineString([[0.0, -10.0], [0.0, 20.0]]) +l3 = LG.LineString([[3.0, 0.0], [8.0, 5.0], [13.0, 0.0], [8.0, -5.0], [3.0, 0.0]]) +r1 = LG.LinearRing([[0.0, 0.0], [5.0, 5.0], [10.0, 0.0], [5.0, -5.0], [0.0, 0.0]]) +r2 = LG.LinearRing([[3.0, 0.0], [8.0, 5.0], [13.0, 0.0], [8.0, -5.0], [3.0, 0.0]]) +r3 = GI.LinearRing([[3.0, 0.0], [8.0, 5.0], [13.0, 0.0], [8.0, -5.0]]) + +@testset_implementations "LineStrings" begin # Equal lines - @test GO.equals(l1, l1) == LG.equals(l1, l1) - @test GO.equals(l2, l2) == LG.equals(l2, l2) + @test GO.equals($l1, $l1) == LG.equals($l1, $l1) + @test GO.equals($l2, $l2) == LG.equals($l2, $l2) # Different lines - @test GO.equals(l1, l2) == GO.equals(l2, l1) == LG.equals(l1, l2) + @test GO.equals($l1, $l2) == GO.equals($l2, $l1) == LG.equals($l1, $l2) - r1 = LG.LinearRing([[0.0, 0.0], [5.0, 5.0], [10.0, 0.0], [5.0, -5.0], [0.0, 0.0]]) - r2 = LG.LinearRing([[3.0, 0.0], [8.0, 5.0], [13.0, 0.0], [8.0, -5.0], [3.0, 0.0]]) - r3 = GI.LinearRing([[3.0, 0.0], [8.0, 5.0], [13.0, 0.0], [8.0, -5.0]]) - l3 = LG.LineString([[3.0, 0.0], [8.0, 5.0], [13.0, 0.0], [8.0, -5.0], [3.0, 0.0]]) # Equal rings - @test GO.equals(r1, r1) == LG.equals(r1, r1) - @test GO.equals(r2, r2) == LG.equals(r2, r2) - # Test equal rings without closing point - @test GO.equals(r2, r3) - @test GO.equals(r3, l3) + @test GO.equals($r1, $r1) == LG.equals($r1, $r1) + @test GO.equals($r2, $r2) == LG.equals($r2, $r2) # Different rings - @test GO.equals(r1, r2) == GO.equals(r2, r1) == LG.equals(r1, r2) + @test GO.equals($r1, $r2) == GO.equals($r2, $r1) == LG.equals($r1, $r2) # Equal linear ring and line string - @test GO.equals(r2, l3) == LG.equals(r2, l3) + @test GO.equals($r2, $l3) == LG.equals($r2, $l3) # Equal line string and line - @test GO.equals(l1, GI.Line([(0.0, 0.0), (0.0, 10.0)])) + @test GO.equals($l1, GI.Line([(0.0, 0.0), (0.0, 10.0)])) end -@testset "Polygons/MultiPolygons" begin - pt1 = LG.Point([0.0, 0.0]) - r1 = GI.LinearRing([(0, 0), (0, 5), (5, 5), (5, 0), (0, 0)]) - p1 = GI.Polygon([[(0, 0), (0, 5), (5, 5), (5, 0), (0, 0)]]) - p2 = GI.Polygon([[(1, 1), (1, 6), (6, 6), (6, 1), (1, 1)]]) - 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], [15.0, 11.0], [25.0, 11.0], [25.0, 1.0], [15.0, 1.0]] - ] - ) - p4 = LG.Polygon( - [ - [[10.0, 0.0], [10.0, 20.0], [30.0, 20.0], [30.0, 0.0], [10.0, 0.0]], - [[16.0, 1.0], [16.0, 11.0], [25.0, 11.0], [25.0, 1.0], [16.0, 1.0]] - ] - ) - p5 = 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]], - [[11.0, 1.0], [11.0, 2.0], [12.0, 2.0], [12.0, 1.0], [11.0, 1.0]] - ] - ) - p6 = GI.Polygon([[(6, 6), (6, 1), (1, 1), (1, 6), (6, 6)]]) - p7 = GI.Polygon([[(6, 6), (1, 6), (1, 1), (6, 1), (6, 6)]]) - p8 = GI.Polygon([[(6, 6), (1, 6), (1, 1), (6, 1)]]) +# LibGEOS rejects rings that are not closed, and they are not eaual in GeometryBasics or ArchGDAL? +@testset_implementations "Rings" [GI] begin + # Test equal rings without closing point + @test GO.equals($r2, $r3) + @test GO.equals($r3, $l3) +end + +pt1 = LG.Point([0.0, 0.0]) +r1 = GI.LinearRing([(0, 0), (0, 5), (5, 5), (5, 0), (0, 0)]) +p1 = GI.Polygon([[(0, 0), (0, 5), (5, 5), (5, 0), (0, 0)]]) +p2 = GI.Polygon([[(1, 1), (1, 6), (6, 6), (6, 1), (1, 1)]]) +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], [15.0, 11.0], [25.0, 11.0], [25.0, 1.0], [15.0, 1.0]] + ] +) +p4 = LG.Polygon( + [ + [[10.0, 0.0], [10.0, 20.0], [30.0, 20.0], [30.0, 0.0], [10.0, 0.0]], + [[16.0, 1.0], [16.0, 11.0], [25.0, 11.0], [25.0, 1.0], [16.0, 1.0]] + ] +) +p5 = 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]], + [[11.0, 1.0], [11.0, 2.0], [12.0, 2.0], [12.0, 1.0], [11.0, 1.0]] + ] +) +p6 = GI.Polygon([[(6, 6), (6, 1), (1, 1), (1, 6), (6, 6)]]) +p7 = GI.Polygon([[(6, 6), (1, 6), (1, 1), (6, 1), (6, 6)]]) +p8 = GI.Polygon([[(6, 6), (1, 6), (1, 1), (6, 1)]]) +p9 = LG.Polygon( + [[ + [-53.57208251953125, 28.287451910503744], + [-53.33038330078125, 28.29228897739706], + [-53.34136962890625, 28.430052892335723], + [-53.57208251953125, 28.287451910503744], + ]] +) +m1 = LG.MultiPolygon([ + [[[0.0, 0.0], [0.0, 5.0], [5.0, 5.0], [5.0, 0.0], [0.0, 0.0]]], + [ + [[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]] + ] +]) +m2 = LG.MultiPolygon([ + [ + [[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]] + ], + [[[0.0, 0.0], [0.0, 5.0], [5.0, 5.0], [5.0, 0.0], [0.0, 0.0]]] +]) +m3 = LG.MultiPolygon([p3]) + +@testset_implementations "Polygons" begin # Point and polygon aren't equal - GO.equals(pt1, p1) == LG.equals(pt1, p1) + GO.equals($pt1, $p1) == LG.equals($pt1, $p1) # Linear ring and polygon aren't equal - @test GO.equals(r1, p1) == LG.equals(r1, p1) + @test GO.equals($r1, $p1) == LG.equals($r1, $p1) # Equal polygon - @test GO.equals(p1, p1) == LG.equals(p1, p1) - @test GO.equals(p2, p2) == LG.equals(p2, p2) + @test GO.equals($p1, $p1) == LG.equals($p1, $p1) + @test GO.equals($p2, $p2) == LG.equals($p2, $p2) # Equal but offset polygons - @test GO.equals(p2, p6) == LG.equals(p2, p6) + @test GO.equals($p2, $p6) == LG.equals($p2, $p6) # Equal but opposite winding orders - @test GO.equals(p2, p7) == LG.equals(p2, p7) - # Equal but without closing point (implied) - @test GO.equals(p7, p8) + @test GO.equals($p2, $p7) == LG.equals($p2, $p7) # Different polygons - @test GO.equals(p1, p2) == LG.equals(p1, p2) + @test GO.equals($p1, $p2) == LG.equals($p1, $p2) # Equal polygons with holes - @test GO.equals(p3, p3) == LG.equals(p3, p3) + @test GO.equals($p3, $p3) == LG.equals($p3, $p3) # Same exterior, different hole - @test GO.equals(p3, p4) == LG.equals(p3, p4) + @test GO.equals($p3, $p4) == LG.equals($p3, $p4) # Same exterior and first hole, has an extra hole - @test GO.equals(p3, p5) == LG.equals(p3, p5) - - p9 = LG.Polygon( - [[ - [-53.57208251953125, 28.287451910503744], - [-53.33038330078125, 28.29228897739706], - [-53.34136962890625, 28.430052892335723], - [-53.57208251953125, 28.287451910503744], - ]] - ) + @test GO.equals($p3, $p5) == LG.equals($p3, $p5) # Complex polygon - @test GO.equals(p9, p9) == LG.equals(p9, p9) + @test GO.equals($p9, $p9) == LG.equals($p9, $p9) +end - m1 = LG.MultiPolygon([ - [[[0.0, 0.0], [0.0, 5.0], [5.0, 5.0], [5.0, 0.0], [0.0, 0.0]]], - [ - [[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]] - ] - ]) - m2 = LG.MultiPolygon([ - [ - [[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]] - ], - [[[0.0, 0.0], [0.0, 5.0], [5.0, 5.0], [5.0, 0.0], [0.0, 0.0]]] - ]) +@testset_implementations "Unclosed Polygons" [GB, GI] begin + # Equal but without closing point (implied) + @test GO.equals($p7, $p8) +end + +@testset_implementations "MultiPolygons" begin # Equal multipolygon - @test GO.equals(m1, m1) == LG.equals(m1, m1) + @test GO.equals($m1, $m1) == LG.equals($m1, $m1) # Equal multipolygon with different order - @test GO.equals(m2, m2) == LG.equals(m2, m2) + @test GO.equals($m2, $m2) == LG.equals($m2, $m2) # Equal polygon to multipolygon - m3 = LG.MultiPolygon([p3]) - @test GO.equals(p1, m3) == LG.equals(p1, m3) -end \ No newline at end of file + @test GO.equals($p1, $m3) == LG.equals($p1, $m3) +end diff --git a/test/methods/geom_relations.jl b/test/methods/geom_relations.jl index ca0568551..0435c10fe 100644 --- a/test/methods/geom_relations.jl +++ b/test/methods/geom_relations.jl @@ -1,6 +1,7 @@ import GeometryOps as GO import GeoInterface as GI import LibGEOS as LG +using ..TestHelpers # Tests of DE-9IM Methods pt1 = LG.Point([0.0, 0.0]) @@ -125,13 +126,13 @@ test_pairs = [ (r1, r1, "r1", "r1", "Same rings"), (r2, r1, "r2", "r1", "Disjoint ring with one 'inside' of hole created"), (r3, r1, "r3", "r1", "Disjoint ring with one 'outside' of hole created"), - (r4, r1, "r4", "r1", "Rings share two sides and rest of sides dont touch"), + (r4, r1, "r4", "r1", "Rings share two sides and rest of sides don't touch"), (r1, r5, "r1", "r5", "Ring shares all edges with other ring, plus an extra loop"), (r1, r6, "r1", "r6", "Rings share just one vertex"), (r1, r7, "r1", "r7", "Rings cross one another"), (r4, p1, "r4", "p1", "Ring on boundary of polygon"), (r1, p1, "r1", "p1", "Ring on boundary and cutting through polygon"), - (r2, p1, "r2", "p1", "Ring on hole bounday"), + (r2, p1, "r2", "p1", "Ring on hole boundary"), (r6, p1, "r6", "p1", "Ring touches polygon at one vertex"), (r7, p1, "r7", "p1", "Ring crosses through polygon"), (r8, p1, "r8", "p1", "Ring inside of polygon"), @@ -154,19 +155,18 @@ test_pairs = [ (mpt1, mpt3, "mpt1", "mpt3", "No shared points"), (ml1, ml2, "ml1", "ml2", "Lines in ml1 cross and touch ml2"), (mp1, mp2, "mp1", "mp2", "Polygons in mp1 are inside hole and overlap"), - (gc1, ml1, "gc1", "ml1", "Make sure collection works with multi-geom"), + # (gc1, ml1, "gc1", "ml1", "Make sure collection works with multi-geom"), ] -function test_geom_relation(GO_f, LG_f, f_name; swap_points = false) +function test_geom_relation(GO_f, LG_f, f_name; swap_points=false) for (g1, g2, sg1, sg2, sdesc) in test_pairs - if swap_points - g1, g2 = g2, g1 - sg1, sg2 = sg2, sg1 + @testset "$sg1 $sg2 $sdesc" begin + if swap_points + @test_implementations (GO_f($g2, $g1) == LG_f($g2, $g1)) + else + @test_implementations GO_f($g1, $g2) == LG_f($g1, $g2) + end end - go_val = GO_f(g1, g2) - lg_val = LG_f(g1, g2) - @test go_val == lg_val - go_val != lg_val && println("\n↑ TEST INFO: $sg1 $f_name $sg2 - $sdesc \n\n") end end @@ -180,118 +180,124 @@ end @testset "Overlaps" begin - @testset "Points/MultiPoints" begin - p1 = LG.Point([0.0, 0.0]) - p2 = LG.Point([0.0, 1.0]) + p1 = LG.Point([0.0, 0.0]) + p2 = LG.Point([0.0, 1.0]) + + mp1 = LG.MultiPoint([[0.0, 1.0], [4.0, 4.0]]) + mp2 = LG.MultiPoint([[0.0, 1.0], [2.0, 2.0]]) + mp3 = LG.MultiPoint([[0.0, 1.0], [2.0, 2.0], [3.0, 3.0]]) + mp4 = LG.MultiPoint([ + [-36.05712890625, 26.480407161007275], + [-35.7220458984375, 27.137368359795584], + [-35.13427734375, 26.83387451505858], + [-35.4638671875, 27.254629577800063], + [-35.5462646484375, 26.86328062676624], + [-35.3924560546875, 26.504988828743404], + ]) + mp5 = GI.MultiPoint([ + [-35.4638671875, 27.254629577800063], + [-35.5462646484375, 26.86328062676624], + [-35.3924560546875, 26.504988828743404], + [-35.2001953125, 26.12091815959972], + [-34.9969482421875, 26.455820238459893], + ]) + + @testset_implementations "Points/MultiPoints" begin # Two points can't overlap - @test GO.overlaps(p1, p1) == LG.overlaps(p1, p2) - - mp1 = LG.MultiPoint([[0.0, 1.0], [4.0, 4.0]]) - mp2 = LG.MultiPoint([[0.0, 1.0], [2.0, 2.0]]) - mp3 = LG.MultiPoint([[0.0, 1.0], [2.0, 2.0], [3.0, 3.0]]) + @test GO.overlaps($p1, $p1) == LG.overlaps($p1, $p2) # No shared points, doesn't overlap - @test GO.overlaps(p1, mp1) == LG.overlaps(p1, mp1) + @test GO.overlaps($p1, $mp1) == LG.overlaps($p1, $mp1) # One shared point, does overlap - @test GO.overlaps(p2, mp1) == LG.overlaps(p2, mp1) + @test GO.overlaps($p2, $mp1) == LG.overlaps($p2, $mp1) # All shared points, doesn't overlap - @test GO.overlaps(mp1, mp1) == LG.overlaps(mp1, mp1) + @test GO.overlaps($mp1, $mp1) == LG.overlaps($mp1, $mp1) # Not all shared points, overlaps - @test GO.overlaps(mp1, mp2) == LG.overlaps(mp1, mp2) + @test GO.overlaps($mp1, $mp2) == LG.overlaps($mp1, $mp2) # One set of points entirely inside other set, doesn't overlap - @test GO.overlaps(mp2, mp3) == LG.overlaps(mp2, mp3) + @test GO.overlaps($mp2, $mp3) == LG.overlaps($mp2, $mp3) # Not all points shared, overlaps - @test GO.overlaps(mp1, mp3) == LG.overlaps(mp1, mp3) - - mp1 = LG.MultiPoint([ - [-36.05712890625, 26.480407161007275], - [-35.7220458984375, 27.137368359795584], - [-35.13427734375, 26.83387451505858], - [-35.4638671875, 27.254629577800063], - [-35.5462646484375, 26.86328062676624], - [-35.3924560546875, 26.504988828743404], - ]) - mp2 = GI.MultiPoint([ - [-35.4638671875, 27.254629577800063], - [-35.5462646484375, 26.86328062676624], - [-35.3924560546875, 26.504988828743404], - [-35.2001953125, 26.12091815959972], - [-34.9969482421875, 26.455820238459893], - ]) + @test GO.overlaps($mp1, $mp3) == LG.overlaps($mp1, $mp3) # Some shared points, overlaps - @test GO.overlaps(mp1, mp2) == LG.overlaps(mp1, mp2) - @test GO.overlaps(mp1, mp2) == GO.overlaps(mp2, mp1) + @test GO.overlaps($mp1, $mp2) == LG.overlaps($mp1, $mp2) + @test GO.overlaps($mp1, $mp2) == GO.overlaps($mp2, $mp1) end - @testset "Lines/Rings" begin - l1 = LG.LineString([[0.0, 0.0], [0.0, 10.0]]) - l2 = LG.LineString([[0.0, -10.0], [0.0, 20.0]]) - l3 = LG.LineString([[0.0, -10.0], [0.0, 3.0]]) - l4 = LG.LineString([[5.0, -5.0], [5.0, 5.0]]) + l1 = LG.LineString([[0.0, 0.0], [0.0, 10.0]]) + l2 = LG.LineString([[0.0, -10.0], [0.0, 20.0]]) + l3 = LG.LineString([[0.0, -10.0], [0.0, 3.0]]) + l4 = LG.LineString([[5.0, -5.0], [5.0, 5.0]]) + # Linear rings that intersect but don't overlap + r1 = LG.LinearRing([[0.0, 0.0], [0.0, 5.0], [5.0, 5.0], [5.0, 0.0], [0.0, 0.0]]) + r2 = LG.LinearRing([[1.0, 1.0], [1.0, 6.0], [6.0, 6.0], [6.0, 1.0], [1.0, 1.0]]) + + @testset_implementations "Lines/Rings" begin # Line can't overlap with itself - @test GO.overlaps(l1, l1) == LG.overlaps(l1, l1) + @test GO.overlaps($l1, $l1) == LG.overlaps($l1, $l1) # Line completely within other line doesn't overlap - @test GO.overlaps(l1, l2) == GO.overlaps(l2, l1) == LG.overlaps(l1, l2) + @test GO.overlaps($l1, $l2) == GO.overlaps($l2, $l1) == LG.overlaps($l1, $l2) # Overlapping lines - @test GO.overlaps(l1, l3) == GO.overlaps(l3, l1) == LG.overlaps(l1, l3) + @test GO.overlaps($l1, $l3) == GO.overlaps($l3, $l1) == LG.overlaps($l1, $l3) # Lines that don't touch - @test GO.overlaps(l1, l4) == LG.overlaps(l1, l4) - # Linear rings that intersect but don't overlap - r1 = LG.LinearRing([[0.0, 0.0], [0.0, 5.0], [5.0, 5.0], [5.0, 0.0], [0.0, 0.0]]) - r2 = LG.LinearRing([[1.0, 1.0], [1.0, 6.0], [6.0, 6.0], [6.0, 1.0], [1.0, 1.0]]) - @test GO.overlaps(r1, r2) == LG.overlaps(r1, r2) + @test GO.overlaps($l1, $l4) == LG.overlaps($l1, $l4) + @test LG.overlaps($r1, $r2) == LG.overlaps($r1, $r2) end - @testset "Polygons/MultiPolygons" begin - p1 = LG.Polygon([[[0.0, 0.0], [0.0, 5.0], [5.0, 5.0], [5.0, 0.0], [0.0, 0.0]]]) - p2 = LG.Polygon([ + p1 = LG.Polygon([[[0.0, 0.0], [0.0, 5.0], [5.0, 5.0], [5.0, 0.0], [0.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([[[1.0, 1.0], [1.0, 6.0], [6.0, 6.0], [6.0, 1.0], [1.0, 1.0]]]) + p4 = LG.Polygon([[[20.0, 5.0], [20.0, 10.0], [18.0, 10.0], [18.0, 5.0], [20.0, 5.0]]]) + p5 = LG.Polygon( + [[ + [-53.57208251953125, 28.287451910503744], + [-53.33038330078125, 28.29228897739706], + [-53.34136352890625, 28.430052892335723], + [-53.57208251953125, 28.287451910503744], + ]] + ) + m1 = LG.MultiPolygon([ + [[[0.0, 0.0], [0.0, 5.0], [5.0, 5.0], [5.0, 0.0], [0.0, 0.0]]], + [ [[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]] - ]) + ] + ]) + + @testset_implementations "Polygons/MultiPolygons" begin # Test basic polygons that don't overlap - @test GO.overlaps(p1, p2) == LG.overlaps(p1, p2) - @test !GO.overlaps(p1, (1, 1)) - @test !GO.overlaps((1, 1), p2) + @test GO.overlaps($p1, $p2) == LG.overlaps($p1, $p2) + @test !GO.overlaps($p1, (1, 1)) + @test !GO.overlaps((1, 1), $p2) - p3 = LG.Polygon([[[1.0, 1.0], [1.0, 6.0], [6.0, 6.0], [6.0, 1.0], [1.0, 1.0]]]) # Test basic polygons that overlap - @test GO.overlaps(p1, p3) == LG.overlaps(p1, p3) + @test GO.overlaps($p1, $p3) == LG.overlaps($p1, $p3) - p4 = LG.Polygon([[[20.0, 5.0], [20.0, 10.0], [18.0, 10.0], [18.0, 5.0], [20.0, 5.0]]]) # Test one polygon within the other - @test GO.overlaps(p2, p4) == GO.overlaps(p4, p2) == LG.overlaps(p2, p4) + @test GO.overlaps($p2, $p4) == GO.overlaps($p4, $p2) == LG.overlaps($p2, $p4) - p5 = LG.Polygon( - [[ - [-53.57208251953125, 28.287451910503744], - [-53.33038330078125, 28.29228897739706], - [-53.34136352890625, 28.430052892335723], - [-53.57208251953125, 28.287451910503744], - ]] - ) # Test equal polygons - @test GO.overlaps(p5, p5) == LG.overlaps(p5, p5) + @test GO.overlaps($p5, $p5) == LG.overlaps($p5, $p5) - # Test multipolygons - m1 = LG.MultiPolygon([ - [[[0.0, 0.0], [0.0, 5.0], [5.0, 5.0], [5.0, 0.0], [0.0, 0.0]]], - [ - [[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]] - ] - ]) # Test polygon that overlaps with multipolygon - @test GO.overlaps(m1, p3) == LG.overlaps(m1, p3) + @test GO.overlaps($m1, $p3) == LG.overlaps($m1, $p3) # Test polygon in hole of multipolygon, doesn't overlap - @test GO.overlaps(m1, p4) == LG.overlaps(m1, p4) + @test GO.overlaps($m1, $p4) == LG.overlaps($m1, $p4) end end + @testset "Crosses" begin line6 = GI.LineString([(1.0, 1.0), (1.0, 2.0), (1.0, 3.0), (1.0, 4.0)]) + line7 = GI.LineString([(1, 1), (1, 2), (1, 3), (1, 4)]) + line8 = GI.LineString([(1, 1), (1, 2), (1, 3), (1, 4)]) poly7 = GI.Polygon([[(-1.0, 2.0), (3.0, 2.0), (3.0, 3.0), (-1.0, 3.0), (-1.0, 2.0)]]) - @test GO.crosses(GI.LineString([(-2.0, 2.0), (4.0, 2.0)]), line6) == true - @test GO.crosses(GI.LineString([(0.5, 2.5), (1.0, 1.0)]), poly7) == true - @test GO.crosses(GI.MultiPoint([(1.0, 2.0), (12.0, 12.0)]), GI.LineString([(1, 1), (1, 2), (1, 3), (1, 4)])) == true - @test GO.crosses(GI.MultiPoint([(1.0, 0.0), (12.0, 12.0)]), GI.LineString([(1, 1), (1, 2), (1, 3), (1, 4)])) == false - @test GO.crosses(GI.LineString([(-2.0, 2.0), (-4.0, 2.0)]), poly7) == false -end \ No newline at end of file + @testset_implementations begin + @test GO.crosses(GI.LineString([(-2.0, 2.0), (4.0, 2.0)]), $line6) == true + @test GO.crosses(GI.LineString([(0.5, 2.5), (1.0, 1.0)]), $poly7) == true + @test GO.crosses(GI.MultiPoint([(1.0, 2.0), (12.0, 12.0)]), $line7) == true + @test GO.crosses(GI.MultiPoint([(1.0, 0.0), (12.0, 12.0)]), $line8) == false + @test GO.crosses(GI.LineString([(-2.0, 2.0), (-4.0, 2.0)]), $poly7) == false + end +end diff --git a/test/methods/orientation.jl b/test/methods/orientation.jl index cc5c7b5c6..fef53ebd4 100644 --- a/test/methods/orientation.jl +++ b/test/methods/orientation.jl @@ -1,34 +1,36 @@ using Test, GeometryOps -import GeoInterface as GI, GeometryOps as GO +import GeoInterface as GI +import GeometryOps as GO +using ..TestHelpers -@testset "Orientation" begin - line1 = GI.LineString([[9.170356, 45.477985], [9.164434, 45.482551], [9.166644, 45.484003]]) - line2 = GI.LineString([[9.169356, 45.477985], [9.163434, 45.482551], [9.165644, 45.484003]]) - line3 = GI.LineString([ - (-111.544189453125, 24.186847428521244), - (-110.687255859375, 24.966140159912975), - (-110.4510498046875, 24.467150664739002), - (-109.9951171875, 25.180087808990645) - ]) - line4 = GI.LineString([ - (-111.4617919921875, 24.05148034322011), - (-110.8795166015625, 24.681961205014595), - (-110.841064453125, 24.14174098050432), - (-109.97863769531249, 24.617057340809524) - ]) +line1 = GI.LineString([[9.170356, 45.477985], [9.164434, 45.482551], [9.166644, 45.484003]]) +line2 = GI.LineString([[9.169356, 45.477985], [9.163434, 45.482551], [9.165644, 45.484003]]) +line3 = GI.LineString([ + (-111.544189453125, 24.186847428521244), + (-110.687255859375, 24.966140159912975), + (-110.4510498046875, 24.467150664739002), + (-109.9951171875, 25.180087808990645) +]) +line4 = GI.LineString([ + (-111.4617919921875, 24.05148034322011), + (-110.8795166015625, 24.681961205014595), + (-110.841064453125, 24.14174098050432), + (-109.97863769531249, 24.617057340809524) +]) - # @test isparallel(line1, line2) == true - # @test isparallel(line3, line4) == false - - poly1 = GI.Polygon([[[0, 0], [1, 0], [1, 1], [0.5, 0.5], [0, 1], [0, 0]]]) - poly2 = GI.Polygon([[[0, 0], [0, 1], [1, 1], [1, 0], [0, 0]]]) +poly1 = GI.Polygon([[[0, 0], [1, 0], [1, 1], [0.5, 0.5], [0, 1], [0, 0]]]) +poly2 = GI.Polygon([[[0, 0], [0, 1], [1, 1], [1, 0], [0, 0]]]) - @test isconcave(poly1) == true - @test isconcave(poly2) == false +l1 = GI.LineString([[0, 0], [1, 1], [1, 0], [0, 0]]) +l2 = GI.LineString([[0, 0], [1, 0], [1, 1], [0, 0]]) - l1 = GI.LineString([[0, 0], [1, 1], [1, 0], [0, 0]]) - l2 = GI.LineString([[0, 0], [1, 0], [1, 1], [0, 0]]) +@testset_implementations "Orientation" begin + # @test isparallel(line1, line2) == true + # @test isparallel(line3, line4) == false + + @test isconcave($poly1) == true + @test isconcave($poly2) == false - @test isclockwise(l1) == true - @test isclockwise(l2) == false + @test isclockwise($l1) == true + @test isclockwise($l2) == false end diff --git a/test/methods/polygonize.jl b/test/methods/polygonize.jl index ea37b15ed..7833486d7 100644 --- a/test/methods/polygonize.jl +++ b/test/methods/polygonize.jl @@ -1,9 +1,47 @@ using Test -import GeometryOps as GO import GeometryOps: polygonize -import GeoInterface as GI import DimensionalData as DD -import OffsetArrays, Rasters +import Rasters, DimensionalData +using GeometryOps, GeoInterface, Test +using ..TestHelpers + +import GeometryOps as GO +import GeoInterface as GI + +@testset "Polygonize with xs and ys, without offsetarrays" begin + @test !(@isdefined OffsetArrays) # to make sure this isn't loaded somewhere else + data = rand(1:4, 100, 100) .== 1 + unitrange = 1:100 + steprange = 1:1:100 + steprangelen = range(1, 100; length = 100) + data_mp = polygonize(data) + for range in (unitrange, steprange, steprangelen) + data_mp_range = polygonize(range, range, data) + @test GO.equals(data_mp, data_mp_range) + end + # ideally we'd have a better test to make sure this returns what we think it does + data_mp_range200 = polygonize(2:2:200, 2:2:200, data) + @test length(GI.coordinates(data_mp_range200)) == length(GI.coordinates(data_mp)) + + # this is an example that could throw floating point error + range_floats = -1.333333333333343:0.041666666666666664:0.374999999999986 + data2 = rand(1:4, length(range_floats), length(range_floats)) .== 1 + data_mp_range_floats = polygonize(range_floats, range_floats, data2) +end + +import OffsetArrays + +@testset "Polygonize with xs and ys, with offsetarrays" begin + data = rand(1:4, 100, 100) .== 1 + unitrange = 1:100 + steprange = 1:1:100 + steprangelen = range(1, 100; length = 100) + data_mp = polygonize(data) + for range in (unitrange, steprange, steprangelen) + data_mp_range = polygonize(range, range, data) + @test GO.equals(data_mp, data_mp_range) + end +end # Missing holes throw a warning, so testing there are # no warnings in a range of randomisation is one way to test @@ -75,3 +113,4 @@ end @test GI.crs(rast_mp) == GI.crs(evil) end end + diff --git a/test/primitives.jl b/test/primitives.jl index 6661a4eea..837916683 100644 --- a/test/primitives.jl +++ b/test/primitives.jl @@ -1,9 +1,15 @@ using Test -import GeoInterface as GI, GeometryOps as GO, GeometryBasics as GB +import ArchGDAL as AG +import GeometryBasics as GB +import GeometryOps as GO +import GeoInterface as GI +import LibGEOS as LG import Proj -import Shapefile, DataFrames +import Shapefile +import DataFrames using Downloads: download +using ..TestHelpers pv1 = [(1, 2), (3, 4), (5, 6), (1, 2)] pv2 = [(3, 4), (5, 6), (6, 7), (3, 4)] @@ -13,13 +19,14 @@ poly = GI.Polygon([lr1, lr2]) @testset "apply" begin - flipped_poly = GO.apply(GI.PointTrait, poly) do p - (GI.y(p), GI.x(p)) + @testset_implementations "apply flip" begin + flipped_poly = GO.apply(GI.PointTrait, $poly) do p + (GI.y(p), GI.x(p)) + end + @test flipped_poly == GI.Polygon([GI.LinearRing([(2, 1), (4, 3), (6, 5), (2, 1)]), + GI.LinearRing([(4, 3), (6, 5), (7, 6), (4, 3)])]) end - @test flipped_poly == GI.Polygon([GI.LinearRing([(2, 1), (4, 3), (6, 5), (2, 1)]), - GI.LinearRing([(4, 3), (6, 5), (7, 6), (4, 3)])]) - @testset "Tables.jl support" begin # check to account for missing data missing_or_equal(x, y) = (ismissing(x) && ismissing(y)) || (x == y) @@ -73,9 +80,8 @@ end @testset "flatten" begin very_wrapped = [[GI.FeatureCollection([GI.Feature(poly; properties=(;))])]] - @test collect(GO.flatten(GI.PointTrait, very_wrapped)) == vcat(pv1, pv2) - @test collect(GO.flatten(GI.LinearRingTrait, [poly])) == [lr1, lr2] - @test collect(GO.flatten(GI.LinearRingTrait, [poly])) == [lr1, lr2] + @test GO._tuple_point.(GO.flatten(GI.PointTrait, very_wrapped)) == vcat(pv1, pv2) + @test collect(GO.flatten(GI.AbstractCurveTrait, [poly])) == [lr1, lr2] @test collect(GO.flatten(GI.x, GI.PointTrait, very_wrapped)) == first.(vcat(pv1, pv2)) end diff --git a/test/runtests.jl b/test/runtests.jl index 567ae3f5c..d4d25637a 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -1,48 +1,35 @@ using GeometryOps using Test +using SafeTestsets -using GeometryOps.GeoInterface -using GeometryOps.GeometryBasics -using GeoInterface.Extents: Extents -using ArchGDAL -using LibGEOS -using Random, Distributions -using Proj +include("helpers.jl") -const GI = GeoInterface -const AG = ArchGDAL -const LG = LibGEOS -const GO = GeometryOps - -@testset "GeometryOps.jl" begin - @testset "Primitives" begin include("primitives.jl") end - # Methods - @testset "Angles" begin include("methods/angles.jl") end - @testset "Area" begin include("methods/area.jl") end - @testset "Barycentric coordinate operations" begin include("methods/barycentric.jl") end - @testset "Orientation" begin include("methods/orientation.jl") end - @testset "Centroid" begin include("methods/centroid.jl") end - @testset "DE-9IM Geom Relations" begin include("methods/geom_relations.jl") end - @testset "Distance" begin include("methods/distance.jl") end - @testset "Equals" begin include("methods/equals.jl") end - # Clipping - @testset "Coverage" begin include("methods/clipping/coverage.jl") end - @testset "Cut" begin include("methods/clipping/cut.jl") end - @testset "Intersection Point" begin include("methods/clipping/intersection_points.jl") end - @testset "Polygon Clipping" begin include("methods/clipping/polygon_clipping.jl") end - # Transformations - @testset "Embed Extent" begin include("transformations/extent.jl") end - @testset "Reproject" begin include("transformations/reproject.jl") end - @testset "Flip" begin include("transformations/flip.jl") end - @testset "Simplify" begin include("transformations/simplify.jl") end - @testset "Segmentize" begin include("transformations/segmentize.jl") end - @testset "Transform" begin include("transformations/transform.jl") end - @testset "Geometry correction" begin - include("transformations/correction/geometry_correction.jl") - include("transformations/correction/closed_ring.jl") - include("transformations/correction/intersecting_polygons.jl") - end - # Extensions - @testset "FlexiJoins" begin include("extensions/flexijoins.jl") end - @testset "LibGEOS" begin include("extensions/libgeos.jl") end -end +@safetestset "Primitives" begin include("primitives.jl") end +# Methods +@safetestset "Angles" begin include("methods/angles.jl") end +@safetestset "Area" begin include("methods/area.jl") end +# @safetestset "Barycentric coordinate operations" begin include("methods/barycentric.jl") end +@safetestset "Orientation" begin include("methods/orientation.jl") end +@safetestset "Centroid" begin include("methods/centroid.jl") end +@safetestset "Convex Hull" begin include("methods/convex_hull.jl") end +@safetestset "DE-9IM Geom Relations" begin include("methods/geom_relations.jl") end +@safetestset "Distance" begin include("methods/distance.jl") end +@safetestset "Equals" begin include("methods/equals.jl") end +# Clipping +@safetestset "Coverage" begin include("methods/clipping/coverage.jl") end +@safetestset "Cut" begin include("methods/clipping/cut.jl") end +@safetestset "Intersection Point" begin include("methods/clipping/intersection_points.jl") end +@safetestset "Polygon Clipping" begin include("methods/clipping/polygon_clipping.jl") end +# Transformations +@safetestset "Embed Extent" begin include("transformations/extent.jl") end +@safetestset "Reproject" begin include("transformations/reproject.jl") end +@safetestset "Flip" begin include("transformations/flip.jl") end +@safetestset "Simplify" begin include("transformations/simplify.jl") end +@safetestset "Segmentize" begin include("transformations/segmentize.jl") end +@safetestset "Transform" begin include("transformations/transform.jl") end +@safetestset "Geometry Correction" begin include("transformations/correction/geometry_correction.jl") end +@safetestset "Closed Rings" begin include("transformations/correction/closed_ring.jl") end +@safetestset "Intersecting Polygons" begin include("transformations/correction/intersecting_polygons.jl") end +# Extensions +@safetestset "FlexiJoins" begin include("extensions/flexijoins.jl") end +@safetestset "LibGEOS" begin include("extensions/libgeos.jl") end diff --git a/test/transformations/correction/closed_ring.jl b/test/transformations/correction/closed_ring.jl index 939d88c0f..90819f40f 100644 --- a/test/transformations/correction/closed_ring.jl +++ b/test/transformations/correction/closed_ring.jl @@ -1,12 +1,16 @@ using Test +import ArchGDAL as AG +import GeoInterface as GI +import GeometryBasics as GB +import GeometryOps as GO +using ..TestHelpers -import GeoInterface as GI, GeometryOps as GO +open_rectangle = GI.Polygon([collect.([(0, 0), (10, 0), (10, 10), (0, 10)])]) -open_rectangle = GI.Wrappers.Polygon([collect.([(0, 0), (10, 0), (10, 10), (0, 10)])]) - -@testset "Closed Ring correction" begin - closed_rectangle = GO.ClosedRing()(open_rectangle) - @test GI.npoint(closed_rectangle) == GI.npoint(open_rectangle) + 1 # test that the rectangle is closed +# LibGEOS fails in GI.convert if the ring is not closed +@testset_implementations "Closed Ring correction" [GI, AG, GB] begin + closed_rectangle = GO.ClosedRing()($open_rectangle) + @test GI.npoint(closed_rectangle) == GI.npoint($open_rectangle) + 1 # test that the rectangle is closed @test GI.getpoint(closed_rectangle.geom[1], 1) == GI.getpoint(closed_rectangle.geom[1], GI.npoint(closed_rectangle)) @test all(GO.flatten(GI.PointTrait, closed_rectangle) .== GO.flatten(GI.PointTrait, GO.ClosedRing()(closed_rectangle))) -end \ No newline at end of file +end diff --git a/test/transformations/correction/geometry_correction.jl b/test/transformations/correction/geometry_correction.jl index e2a1a48ac..de8d51a09 100644 --- a/test/transformations/correction/geometry_correction.jl +++ b/test/transformations/correction/geometry_correction.jl @@ -1,12 +1,16 @@ using Test +import ArchGDAL as AG +import GeoInterface as GI +import GeometryBasics as GB +import GeometryOps as GO +using ..TestHelpers -import GeoInterface as GI, GeometryOps as GO +open_rectangle = GI.Polygon([collect.([(0, 0), (10, 0), (10, 10), (0, 10)])]) -open_rectangle = GI.Wrappers.Polygon([collect.([(0, 0), (10, 0), (10, 10), (0, 10)])]) - -@testset "Fix method" begin - closed_rectangle = GO.fix(open_rectangle) - @test GI.npoint(closed_rectangle) == GI.npoint(open_rectangle) + 1 # test that the rectangle is closed +# LibGEOS fails in GI.convert if the ring is not closed +@testset_implementations "fix open rectangle" [GI, AG, GB] begin + closed_rectangle = GO.fix($open_rectangle) + @test GI.npoint(closed_rectangle) == GI.npoint($open_rectangle) + 1 # test that the rectangle is closed @test GI.getpoint(closed_rectangle.geom[1], 1) == GI.getpoint(closed_rectangle.geom[1], GI.npoint(closed_rectangle)) @test all(GO.flatten(GI.PointTrait, closed_rectangle) .== GO.flatten(GI.PointTrait, GO.fix(closed_rectangle))) -end \ No newline at end of file +end diff --git a/test/transformations/correction/intersecting_polygons.jl b/test/transformations/correction/intersecting_polygons.jl index c55618b72..dfac723bf 100644 --- a/test/transformations/correction/intersecting_polygons.jl +++ b/test/transformations/correction/intersecting_polygons.jl @@ -1,14 +1,14 @@ using Test - import GeoInterface as GI import GeometryOps as GO +using ..TestHelpers p1 = GI.Polygon([[(0.0, 0.0), (3.0, 0.0), (3.0, 3.0), (0.0, 3.0), (0.0, 0.0)]]) # (p1, p2) -> one polygon inside of the other, sharing an edge p2 = GI.Polygon([[(1.0, 0.0), (2.0, 0.0), (2.0, 1.0), (1.0, 1.0), (1.0, 0.0)]]) # (p1, p3) -> polygons outside of one another, sharing an edge p3 = GI.Polygon([[(1.0, 0.0), (2.0, 0.0), (2.0, -1.0), (1.0, -1.0), (1.0, 0.0)]]) -# (p1, p4) -> polygons are completly disjoint (no holes) +# (p1, p4) -> polygons are completely disjoint (no holes) p4 = GI.Polygon([[(1.0, -1.0), (2.0, -1.0), (2.0, -2.0), (1.0, -2.0), (1.0, -1.0)]]) # (p1, p5) -> polygons cut through one another p5 = GI.Polygon([[(1.0, -1.0), (2.0, -1.0), (2.0, 4.0), (1.0, 4.0), (1.0, -1.0)]]) @@ -25,55 +25,61 @@ mp6 = GI.MultiPolygon([ # four interlocking polygons forming a hole [[(-5.0, -5.0), (-5.0, 15.0), (0.0, 15.0), (0.0, -5.0), (-5.0, -5.0)]], ]) -union_fixed_mp1 = GO.UnionIntersectingPolygons()(mp1) -@test GI.npolygon(union_fixed_mp1) == 1 -@test GO.equals(GI.getpolygon(union_fixed_mp1, 1), p1) +@testset_implementations begin + union_fixed_mp1 = GO.UnionIntersectingPolygons()($mp1) + @test GI.npolygon(union_fixed_mp1) == 1 + @test GO.equals(GI.getpolygon(union_fixed_mp1, 1), $p1) -diff_fixed_mp1 = GO.DiffIntersectingPolygons()(mp1) -@test GO.equals(diff_fixed_mp1, union_fixed_mp1) + diff_fixed_mp1 = GO.DiffIntersectingPolygons()($mp1) + @test GO.equals(diff_fixed_mp1, union_fixed_mp1) -union_fixed_mp2 = GO.UnionIntersectingPolygons()(mp2) -@test GI.npolygon(union_fixed_mp2) == 1 -@test GO.equals(GI.getpolygon(union_fixed_mp2, 1), p1) + union_fixed_mp2 = GO.UnionIntersectingPolygons()($mp2) + @test GI.npolygon(union_fixed_mp2) == 1 + @test GO.equals(GI.getpolygon(union_fixed_mp2, 1), $p1) -diff_fixed_mp2 = GO.DiffIntersectingPolygons()(mp2) -@test GO.equals(diff_fixed_mp2, union_fixed_mp2) + diff_fixed_mp2 = GO.DiffIntersectingPolygons()($mp2) + @test GO.equals(diff_fixed_mp2, union_fixed_mp2) -union_fixed_mp3 = GO.UnionIntersectingPolygons()(mp3) -@test GI.npolygon(union_fixed_mp3) == 1 -@test all((GO.coveredby(p, union_fixed_mp3) for p in GI.getpolygon(mp3))) -diff_polys = GO.difference(union_fixed_mp3, mp3; target = GI.PolygonTrait(), fix_multipoly = nothing) -@test sum(GO.area, diff_polys; init = 0.0) == 0 + union_fixed_mp3 = GO.UnionIntersectingPolygons()($mp3) + @test GI.npolygon(union_fixed_mp3) == 1 + @test all((GO.coveredby(p, union_fixed_mp3) for p in GI.getpolygon(mp3))) + diff_polys = GO.difference(union_fixed_mp3, $mp3; + target = GI.PolygonTrait(), fix_multipoly = nothing + ) + @test sum(GO.area, diff_polys; init = 0.0) == 0 -diff_fixed_mp3 = GO.DiffIntersectingPolygons()(mp3) -@test GI.npolygon(diff_fixed_mp3) == 3 -@test all((GO.coveredby(p, union_fixed_mp3) for p in GI.getpolygon(diff_fixed_mp3))) + diff_fixed_mp3 = GO.DiffIntersectingPolygons()($mp3) + @test GI.npolygon(diff_fixed_mp3) == 3 + @test all((GO.coveredby(p, union_fixed_mp3) for p in GI.getpolygon(diff_fixed_mp3))) -union_fixed_mp4 = GO.UnionIntersectingPolygons()(mp4) -@test GI.npolygon(union_fixed_mp4) == 2 -@test (GO.equals(GI.getpolygon(union_fixed_mp4, 1), p1) && GO.equals(GI.getpolygon(union_fixed_mp4, 2), p4)) || - (GO.equals(GI.getpolygon(union_fixed_mp4, 2), p1) && GO.equals(GI.getpolygon(union_fixed_mp4, 1), p4)) + union_fixed_mp4 = GO.UnionIntersectingPolygons()($mp4) + @test GI.npolygon(union_fixed_mp4) == 2 + @test (GO.equals(GI.getpolygon(union_fixed_mp4, 1), $p1) && GO.equals(GI.getpolygon(union_fixed_mp4, 2), $p4)) || + (GO.equals(GI.getpolygon(union_fixed_mp4, 2), $p1) && GO.equals(GI.getpolygon(union_fixed_mp4, 1), $p4)) -diff_fixed_mp4 = GO.DiffIntersectingPolygons()(mp4) -@test GO.equals(diff_fixed_mp4, union_fixed_mp4) + diff_fixed_mp4 = GO.DiffIntersectingPolygons()($mp4) + @test GO.equals(diff_fixed_mp4, union_fixed_mp4) -union_fixed_mp5 = GO.UnionIntersectingPolygons()(mp5) -@test GI.npolygon(union_fixed_mp5) == 1 -@test all((GO.coveredby(p, union_fixed_mp5) for p in GI.getpolygon(mp5))) -diff_polys = GO.difference(union_fixed_mp5, mp5; target = GI.PolygonTrait(), fix_multipoly = nothing) -@test sum(GO.area, diff_polys; init = 0.0) == 0 + union_fixed_mp5 = GO.UnionIntersectingPolygons()($mp5) + @test GI.npolygon(union_fixed_mp5) == 1 + @test all((GO.coveredby(p, union_fixed_mp5) for p in GI.getpolygon($mp5))) + diff_polys = GO.difference(union_fixed_mp5, $mp5; target = GI.PolygonTrait(), fix_multipoly = nothing) + @test sum(GO.area, diff_polys; init = 0.0) == 0 -diff_fixed_mp5 = GO.DiffIntersectingPolygons()(mp5) -@test GI.npolygon(diff_fixed_mp5) == 3 -@test all((GO.coveredby(p, union_fixed_mp5) for p in GI.getpolygon(diff_fixed_mp5))) + diff_fixed_mp5 = GO.DiffIntersectingPolygons()($mp5) + @test GI.npolygon(diff_fixed_mp5) == 3 + @test all((GO.coveredby(p, union_fixed_mp5) for p in GI.getpolygon(diff_fixed_mp5))) -union_fixed_mp6 = GO.UnionIntersectingPolygons()(mp6) -@test GI.npolygon(union_fixed_mp6) == 1 -@test GI.nhole(GI.getpolygon(union_fixed_mp6, 1)) == 1 -@test all((GO.coveredby(p, union_fixed_mp6) for p in GI.getpolygon(mp6))) -diff_polys = GO.difference(union_fixed_mp6, mp6; target = GI.PolygonTrait(), fix_multipoly = nothing) -@test sum(GO.area, diff_polys; init = 0.0) == 0 + union_fixed_mp6 = GO.UnionIntersectingPolygons()($mp6) + @test GI.npolygon(union_fixed_mp6) == 1 + @test GI.nhole(GI.getpolygon(union_fixed_mp6, 1)) == 1 + @test all((GO.coveredby(p, union_fixed_mp6) for p in GI.getpolygon(mp6))) + diff_polys = GO.difference(union_fixed_mp6, $mp6; + target = GI.PolygonTrait(), fix_multipoly = nothing + ) + @test sum(GO.area, diff_polys; init = 0.0) == 0 -diff_fixed_mp6 = GO.DiffIntersectingPolygons()(mp6) -@test GI.npolygon(diff_fixed_mp6) == 4 -@test all((GO.coveredby(p, union_fixed_mp6) for p in GI.getpolygon(diff_fixed_mp6))) + diff_fixed_mp6 = GO.DiffIntersectingPolygons()($mp6) + @test GI.npolygon(diff_fixed_mp6) == 4 + @test all((GO.coveredby(p, union_fixed_mp6) for p in GI.getpolygon(diff_fixed_mp6))) +end diff --git a/test/transformations/extent.jl b/test/transformations/extent.jl index e93390347..87026e9eb 100644 --- a/test/transformations/extent.jl +++ b/test/transformations/extent.jl @@ -1,15 +1,17 @@ using Test - -import GeoInterface as GI, GeometryOps as GO +import GeoInterface as GI +import GeometryOps as GO using GeoInterface: Extents +using ..TestHelpers -@testset "embed_extent" begin - poly = GI.Polygon([GI.LinearRing([(1, 2), (3, 4), (5, 6), (1, 2)]), - GI.LinearRing([(3, 4), (5, 6), (6, 7), (3, 4)])]) +poly = GI.Polygon([GI.LinearRing([(1, 2), (3, 4), (5, 6), (1, 2)]), + GI.LinearRing([(3, 4), (5, 6), (6, 7), (3, 4)])]) - ext_poly = GO.embed_extent(poly) +@testset_implementations "embed_extent" begin + ext_poly = GO.embed_extent($poly) lr1, lr2 = GI.getgeom(ext_poly) @test ext_poly.extent == Extents.Extent(X=(1, 6), Y=(2, 7)) @test lr1.extent == Extents.Extent(X=(1, 5), Y=(2, 6)) @test lr2.extent == Extents.Extent(X=(3, 6), Y=(4, 7)) end + diff --git a/test/transformations/flip.jl b/test/transformations/flip.jl index e977c7504..f02083ff5 100644 --- a/test/transformations/flip.jl +++ b/test/transformations/flip.jl @@ -1,12 +1,11 @@ using Test - -import GeoInterface as GI, GeometryOps as GO +import GeoInterface as GI +import GeometryOps as GO +using ..TestHelpers -@testset "flip" begin - geom = GI.Polygon([GI.LinearRing([(1, 2), (3, 4), (5, 6), (1, 2)]), - GI.LinearRing([(3, 4), (5, 6), (6, 7), (3, 4)])]) - - - @test GO.flip(geom) == GI.Polygon([GI.LinearRing([(2, 1), (4, 3), (6, 5), (2, 1)]), - GI.LinearRing([(4, 3), (6, 5), (7, 6), (4, 3)])]) +geom = GI.Polygon([GI.LinearRing([(1, 2), (3, 4), (5, 6), (1, 2)]), + GI.LinearRing([(3, 4), (5, 6), (6, 7), (3, 4)])]) +@testset_implementations "flip" begin + @test GO.flip($geom) == GI.Polygon([GI.LinearRing([(2, 1), (4, 3), (6, 5), (2, 1)]), + GI.LinearRing([(4, 3), (6, 5), (7, 6), (4, 3)])]) end diff --git a/test/transformations/reproject.jl b/test/transformations/reproject.jl index ac7d22faf..52044a6da 100644 --- a/test/transformations/reproject.jl +++ b/test/transformations/reproject.jl @@ -1,38 +1,43 @@ using Test - -import GeoInterface as GI, GeometryOps as GO using GeoFormatTypes +import GeoInterface as GI +import GeometryOps as GO import Proj +using ..TestHelpers -@testset "reproject" begin - ring1 = GI.LinearRing([(1, 2), (7, 4), (5, 6), (1, 2)]) - ring2 = GI.LinearRing([(11, 2), (20, 4), (15, 6), (11, 2)]) - hole2 = GI.LinearRing([(14, 4), (16, 4), (17, 5), (14, 4)]) +ring1 = GI.LinearRing([(1, 2), (7, 4), (5, 6), (1, 2)]) +ring2 = GI.LinearRing([(11, 2), (20, 4), (15, 6), (11, 2)]) +hole2 = GI.LinearRing([(14, 4), (16, 4), (17, 5), (14, 4)]) - # Set up a regular tranformation of the points for reference - source_crs = convert(Proj.CRS, EPSG(4326)) - target_crs = convert(Proj.CRS, EPSG(3857)) - trans = Proj.Transformation(source_crs, target_crs; always_xy=true) +# Set up a regular tranformation of the points for reference +source_crs = convert(Proj.CRS, EPSG(4326)) +target_crs = convert(Proj.CRS, EPSG(3857)) +trans = Proj.Transformation(source_crs, target_crs; always_xy=true) - polygon1 = GI.Polygon([ring1]) - polygon2 = GI.Polygon([ring2, hole2]) - multipolygon = GI.MultiPolygon([polygon1, polygon2]) +polygon1 = GI.Polygon([ring1]) +polygon2 = GI.Polygon([ring2, hole2]) +multipolygon = GI.MultiPolygon([polygon1, polygon2]) - ref_points3857 = map(GI.getpoint(multipolygon)) do p - trans([GI.x(p), GI.y(p)]) - end +ref_points3857 = map(GI.getpoint(multipolygon)) do p + trans([GI.x(p), GI.y(p)]) +end + +_xy(p) = GI.x(p), GI.y(p) - multipolygon3857 = GO.reproject(multipolygon, EPSG(4326), EPSG(3857)) - multipolygon4326 = GO.reproject(multipolygon3857; target_crs=EPSG(4326)) - points4326_1 = collect(GI.getpoint(multipolygon)) - points4326_2 = GI.getcoord.(GI.getpoint(multipolygon4326)) - points3857 = GI.getcoord.(GI.getpoint(multipolygon3857)) +@testset_implementations "reproject" begin + multipolygon3857 = GO.reproject($multipolygon, EPSG(4326), EPSG(3857)) + multipolygon4326 = GO.reproject($multipolygon; source_crs=EPSG(4326), target_crs=EPSG(4326)) + points4326_1 = collect(GI.getpoint($multipolygon)) + points4326_2 = collect.(GI.getcoord.(GI.getpoint(multipolygon4326))) + points3857 = collect.(GI.getcoord.(GI.getpoint(multipolygon3857))) # Comparison to regular `trans` on points - @test all(map((p1, p2) -> all(map(isapprox, p1, p2)), ref_points3857, points3857)) + @test map(ref_points3857, points3857) do p1, p2 + all(map(isapprox, _xy(p1), _xy(p2))) + end |> all # Round trip comparison - @test all(map((p1, p2) -> all(map(isapprox, p1, p2)), points4326_1, points4326_2)) + @test all(map((p1, p2) -> all(map(isapprox, _xy(p1), _xy(p2))), points4326_1, points4326_2)) # Embedded crs check @test GI.crs(multipolygon3857) == EPSG(3857) @@ -73,6 +78,5 @@ import Proj GO.reproject(multipolygon4326; target_crs=utm32_wkt) GO.reproject(multipolygon4326; target_crs=ProjString("+proj=moll")) - end diff --git a/test/transformations/segmentize.jl b/test/transformations/segmentize.jl index 3a5d18333..5095c4c84 100644 --- a/test/transformations/segmentize.jl +++ b/test/transformations/segmentize.jl @@ -1,6 +1,8 @@ using Test - -import GeometryOps as GO, GeoInterface as GI +using Proj +import GeometryOps as GO +import GeoInterface as GI +using ..TestHelpers @testset "Segmentation on multiple geometry levels" begin ls = GI.LineString([(0, 0), (1, 1), (2, 2), (3, 3)]) @@ -8,60 +10,60 @@ import GeometryOps as GO, GeoInterface as GI p = GI.Polygon([lr]) mp = GI.MultiPolygon([p, p, p]) mls = GI.MultiLineString([ls, ls, ls]) - - @testset "LinearSegments" begin - @test GO.segmentize(ls; max_distance = 0.5) isa GI.LineString - @test GO.segmentize(lr; max_distance = 0.5) isa GI.LinearRing + @testset_implementations "LinearSegments" begin + @test GO.segmentize($ls; max_distance = 0.5) isa GI.LineString + if GI.trait($lr) isa GI.LinearRingTrait + @test GO.segmentize($lr; max_distance = 0.5) isa GI.LinearRing + end # Test that linear rings are closed after segmentization - segmentized_ring = GO.segmentize(lr; max_distance = 0.5) - @test GI.getpoint(segmentized_ring, 1) == GI.getpoint(segmentized_ring, GI.npoint(segmentized_ring)) + segmentized_ring = GO.segmentize($lr; max_distance = 0.5) + @test GI.getpoint(segmentized_ring, 1) == + GI.getpoint(segmentized_ring, GI.npoint(segmentized_ring)) - @test GO.segmentize(p; max_distance = 0.5) isa GI.Polygon - @test GO.segmentize(mp; max_distance = 0.5) isa GI.MultiPolygon - @test GI.ngeom(GO.segmentize(mp; max_distance = 0.5)) == 3 + @test GO.segmentize($p; max_distance = 0.5) isa GI.Polygon + @test GO.segmentize($mp; max_distance = 0.5) isa GI.MultiPolygon + @test GI.ngeom(GO.segmentize($mp; max_distance = 0.5)) == 3 # Now test multilinestrings - @test GO.segmentize(mls; max_distance = 0.5) isa GI.MultiLineString - @test GI.ngeom(GO.segmentize(mls; max_distance = 0.5)) == 3 + @test GO.segmentize($mls; max_distance = 0.5) isa GI.MultiLineString + @test GI.ngeom(GO.segmentize($mls; max_distance = 0.5)) == 3 end - @testset "GeodesicSegments" begin - - @test GO.segmentize(GO.GeodesicSegments(; max_distance = 0.5*900), ls) isa GI.LineString - @test GO.segmentize(GO.GeodesicSegments(; max_distance = 0.5*900), lr) isa GI.LinearRing + @testset_implementations "GeodesicSegments" begin + @test GO.segmentize(GO.GeodesicSegments(; max_distance = 0.5*900), $ls) isa GI.LineString + if GI.trait($lr) isa GI.LinearRingTrait + @test GO.segmentize(GO.GeodesicSegments(; max_distance = 0.5*900), $lr) isa GI.LinearRing + end # Test that linear rings are closed after segmentization - segmentized_ring = GO.segmentize(GO.GeodesicSegments(; max_distance = 0.5*900), lr) + segmentized_ring = GO.segmentize(GO.GeodesicSegments(; max_distance = 0.5*900), $lr) @test GI.getpoint(segmentized_ring, 1) == GI.getpoint(segmentized_ring, GI.npoint(segmentized_ring)) - p = GI.Polygon([lr]) - mp = GI.MultiPolygon([p, p, p]) - @test GO.segmentize(GO.GeodesicSegments(; max_distance = 0.5*900), p) isa GI.Polygon - @test GO.segmentize(GO.GeodesicSegments(; max_distance = 0.5*900), mp) isa GI.MultiPolygon - @test GI.ngeom(GO.segmentize(GO.GeodesicSegments(; max_distance = 0.5*900), mp)) == 3 + @test GO.segmentize(GO.GeodesicSegments(; max_distance = 0.5*900), $p) isa GI.Polygon + @test GO.segmentize(GO.GeodesicSegments(; max_distance = 0.5*900), $mp) isa GI.MultiPolygon + @test GI.ngeom(GO.segmentize(GO.GeodesicSegments(; max_distance = 0.5*900), $mp)) == 3 # Now test multilinestrings - mls = GI.MultiLineString([ls, ls, ls]) - @test GO.segmentize(GO.GeodesicSegments(; max_distance = 0.5*900), mls) isa GI.MultiLineString - @test GI.ngeom(GO.segmentize(GO.GeodesicSegments(; max_distance = 0.5*900), mls)) == 3 + @test GO.segmentize(GO.GeodesicSegments(; max_distance = 0.5*900), $mls) isa GI.MultiLineString + @test GI.ngeom(GO.segmentize(GO.GeodesicSegments(; max_distance = 0.5*900), $mls)) == 3 end end -@testset "LinearSegments" begin - lr = GI.LinearRing([(0, 0), (1, 0), (1, 1), (0, 1), (0, 0)]) - ct = GO.centroid(lr) - ar = GO.area(lr) +lr = GI.LinearRing([(0, 0), (1, 0), (1, 1), (0, 1), (0, 0)]) +@testset_implementations "LinearSegments" begin + ct = GO.centroid($lr) + ar = GO.area($lr) for max_distance in exp10.(LinRange(log10(0.01), log10(1), 10)) - segmentized = GO.segmentize(GO.LinearSegments(; max_distance), lr) + segmentized = GO.segmentize(GO.LinearSegments(; max_distance), $lr) @test all(GO.centroid(segmentized) .≈ ct) @test GO.area(segmentized) ≈ ar end end -@testset "GeodesicSegments" begin - lr = GI.LinearRing([(0, 0), (1, 0), (1, 1), (0, 1), (0, 0)]) +lr = GI.LinearRing([(0, 0), (1, 0), (1, 1), (0, 1), (0, 0)]) +@testset_implementations "GeodesicSegments" begin for max_distance in exp10.(LinRange(log10(0.01), log10(1), 10)) .* 900 - @test_nowarn segmentized = GO.segmentize(GO.GeodesicSegments(; max_distance), lr) + @test_nowarn segmentized = GO.segmentize(GO.GeodesicSegments(; max_distance), $lr) end -end \ No newline at end of file +end diff --git a/test/transformations/simplify.jl b/test/transformations/simplify.jl index dd59ef56c..1661c566a 100644 --- a/test/transformations/simplify.jl +++ b/test/transformations/simplify.jl @@ -1,6 +1,9 @@ using Test import GeoJSON, JLD2 -import GeometryOps as GO, GeoInterface as GI, LibGEOS as LG +import GeometryOps as GO +import GeoInterface as GI +import LibGEOS as LG +using ..TestHelpers datadir = realpath(joinpath(dirname(pathof(GO)), "../test/data")) @testset "RadialDistance and VisvalingamWhyatt" begin @@ -8,6 +11,7 @@ datadir = realpath(joinpath(dirname(pathof(GO)), "../test/data")) fc2 = GeoJSON.read(joinpath(datadir, "simplify2.geojson")) fcs = [fc for i in 1:100] + # TODO: @testset_implementations doesn't handle feature collections yet for T in (GO.RadialDistance, GO.VisvalingamWhyatt) @test length(collect(GO.flatten(GI.PointTrait, GO.simplify(T(number=10), fc)))) == 10 @test length(collect(GO.flatten(GI.PointTrait, GO.simplify(T(ratio=0.5), fc)))) == 39 # Half of 78 @@ -15,20 +19,25 @@ datadir = realpath(joinpath(dirname(pathof(GO)), "../test/data")) GO.simplify(T(tol=0.001), fcs; threaded=true, calc_extent=true) end end + @testset "DouglasPeucker" begin poly_coords = JLD2.jldopen(joinpath(datadir, "complex_polygons.jld2"))["verts"][1:4] for c in poly_coords npoints = length(c[1]) poly = LG.Polygon(c) - lg_vals = GI.coordinates(LG.simplify(poly, 100.0))[1] - reduced_npoints = length(lg_vals) - @test all(GI.coordinates(GO.simplify(poly; tol = 100.0))[1] .== lg_vals) - @test all(GI.coordinates(GO.simplify(poly; number = reduced_npoints))[1] .== lg_vals) - @test all(GI.coordinates(GO.simplify(poly; ratio = (reduced_npoints/npoints)))[1] .== lg_vals) + @testset_implementations "Polygon coords match LibGEOS simplify" begin + lg_vals = GI.coordinates(LG.simplify($poly, 100.0))[1] + reduced_npoints = length(lg_vals) + @test all(GI.coordinates(GO.simplify($poly; tol = 100.0))[1] .== lg_vals) + @test all(GI.coordinates(GO.simplify($poly; number = reduced_npoints))[1] .== lg_vals) + @test all(GI.coordinates(GO.simplify($poly; ratio = (reduced_npoints/npoints)))[1] .== lg_vals) + end end # Ensure last point isn't removed with curve c = poly_coords[1] - line = LG.LineString(c[1]) - lg_vals = GI.coordinates(LG.simplify(line, 100.0)) - @test all(GI.coordinates(GO.simplify(line; tol = 100.0)) .== lg_vals) + linestring = LG.LineString(c[1]) + @testset_implementations "LineString coords match LibGEOS simplify" begin + lg_vals = GI.coordinates(LG.simplify($linestring, 100.0)) + @test all(GI.coordinates(GO.simplify($linestring; tol = 100.0)) .== lg_vals) + end end diff --git a/test/transformations/transform.jl b/test/transformations/transform.jl index f3d5b50ef..40a700e65 100644 --- a/test/transformations/transform.jl +++ b/test/transformations/transform.jl @@ -1,13 +1,15 @@ using Test - -import GeoInterface as GI, GeometryOps as GO using CoordinateTransformations +import GeoInterface as GI +import GeometryOps as GO +using ..TestHelpers -@testset "transform" begin - geom = GI.Polygon([GI.LinearRing([(1, 2), (3, 4), (5, 6), (1, 2)]), - GI.LinearRing([(3, 4), (5, 6), (6, 7), (3, 4)])]) +geom = GI.Polygon([GI.LinearRing([(1, 2), (3, 4), (5, 6), (1, 2)]), + GI.LinearRing([(3, 4), (5, 6), (6, 7), (3, 4)])]) + +@testset_implementations "transform" begin translated = GI.Polygon([GI.LinearRing([[4.5, 3.5], [6.5, 5.5], [8.5, 7.5], [4.5, 3.5]]), GI.LinearRing([[6.5, 5.5], [8.5, 7.5], [9.5, 8.5], [6.5, 5.5]])]) f = CoordinateTransformations.Translation(3.5, 1.5) - @test GO.transform(f, geom) == translated + @test GO.transform(f, $geom) == translated end