[WIP] Vancouver polys demo

docs/Project.toml

@@ -17,6 +17,7 @@ Downloads = "f43a241f-c20a-4ad4-852c-f6b1247861c6"
 ExactPredicates = "429591f6-91af-11e9-00e2-59fbe8cec110"
 FlexiJoins = "e37f2e79-19fa-4eb7-8510-b63b51fe0a37"
 GADM = "a8dd9ffe-31dc-4cf5-a379-ea69100a8233"
+GZip = "92fee26a-97fe-5a0c-ad85-20a5f3185b63"
 GeoDataFrames = "62cb38b5-d8d2-4862-a48e-6a340996859f"
 GeoDatasets = "ddc7317b-88db-5cb5-a849-8449e5df04f9"
 GeoFormatTypes = "68eda718-8dee-11e9-39e7-89f7f65f511f"
@@ -41,3 +42,4 @@ Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
 Shapefile = "8e980c4a-a4fe-5da2-b3a7-4b4b0353a2f4"
 Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"
 Weave = "44d3d7a6-8a23-5bf8-98c5-b353f8df5ec9"
+WellKnownGeometry = "0f680547-7be7-4555-8820-bb198eeb646b"

docs/make.jl

@@ -93,16 +93,19 @@ makedocs(;
     ),
     pages=[
         "Introduction" => "introduction.md",
-        "API Reference" => "api.md",
         "Tutorials" => [
             "Creating Geometry" => "tutorials/creating_geometry.md",
             "Spatial Joins" => "tutorials/spatial_joins.md",
         ],
+        "How-To" => [
+            "Unary union" => "howto/unary_union.md",
+        ],
         "Explanations" => [
             "Paradigms" => "explanations/paradigms.md",
             "Peculiarities" => "explanations/peculiarities.md",
         ],
         "Source code" => literate_pages,
+        "API Reference" => "api.md",
     ],
     warnonly = true,
 )

docs/src/howto/unary_union.md

@@ -0,0 +1,45 @@
+# Unary union
+
+```@example unary
+using Makie, GeoMakie
+import GeometryOps as GO, GeoInterface as GI
+```
+
+Furst, we get the data.  These are watersheds in Vancouver Island, Canada - the very same used in the GEOS benchmarks.
+
+Since the file is zipped, we have to unzip it, which is why this code is a bit longer than it otherwise might be.
+
+```@example unary
+wkt_gz_file = download("https://rawcdn.githack.com/pramsey/geos-performance/b54d92a678e2174059d1b0ff233e275e4bd02084/data/watersheds.wkt.gz", joinpath(tempdir(), "watersheds.wkt.gz"))
+using GZip
+handle = GZip.open(wkt_gz_file)
+
+using WellKnownGeometry, GeoFormatTypes
+wkt = GeoFormatTypes.WellKnownText.((GeoFormatTypes.Geom(),), eachline(handle))
+close(handle)
+
+geoms = GO.tuples(wkt)
+
+plot(geoms; color = 1:length(geoms), axis = (; aspect = DataAspect()))
+```
+
+Now that we have the geometries, we reduce over the vector, performing unions along the way.
+
+```@example unary
+fixed_geoms = GO.buffer(geoms, 0)
+@time GO.fix(GI.MultiPolygon(fixed_geoms); corrections = (GO.UnionIntersectingPolygons(),))
+```
+
+```@example unary
+@time final_multipoly = reduce(
+    (x, y) -> GO.union(x, y; target = GI.MultiPolygonTrait, fix = GO.UnionIntersectingPolygons()), 
+    GO.fix(geoms)
+)
+```
+
+```@example unary
+@time final_multipoly = reduce(
+    (x, y) -> GO.union(x, y; target = GI.MultiPolygonTrait, fix = GO.UnionIntersectingPolygons()), 
+    fixed_geoms
+)
+```

src/GeometryOps.jl

@@ -60,6 +60,7 @@ include("transformations/transform.jl")
 include("transformations/correction/geometry_correction.jl")
 include("transformations/correction/closed_ring.jl")
 include("transformations/correction/intersecting_polygons.jl")
+include("transformations/correction/polygon_contents.jl")
 
 # Import all names from GeoInterface and Extents, so users can do `GO.extent` or `GO.trait`.
 for name in names(GeoInterface)

src/methods/clipping/union.jl

@@ -88,6 +88,20 @@ function _union(
     return polys
 end
 
+function _union(
+    ::TraitTarget{GI.MultiPolygonTrait}, ::Type{T},
+    ::GI.PolygonTrait, poly_a,
+    ::GI.PolygonTrait, poly_b;
+    exact, kwargs...,
+) where T
+    return GI.MultiPolygon(_union(
+        TraitTarget(GI.PolygonTrait()), T, 
+        GI.trait(poly_a), poly_a, 
+        GI.trait(poly_b), poly_b;
+        exact, kwargs...,
+    ))
+end
+
 # # Helper functions for Unions with Greiner and Hormann Polygon Clipping
 
 #= When marking the crossing status of a delayed crossing, the chain start point is crossing
@@ -217,7 +231,7 @@ function _union(
     if !isnothing(fix_multipoly) # Fix multipoly_b to prevent repeated regions in the output
         multipoly_b = fix_multipoly(multipoly_b)
     end
-    polys = [tuples(poly_a, T)]
+    polys = [fix(tuples(poly_a, T); corrections = (PolygonContents(), ))]
     for poly_b in GI.getpolygon(multipoly_b)
         if intersects(polys[1], poly_b)
             # If polygons intersect and form a new polygon, swap out polygon
@@ -229,11 +243,24 @@ function _union(
             end
         else
             # If they don't intersect, poly_b is now a part of the union as its own polygon
-            push!(polys, tuples(poly_b, T))
+            push!(polys, fix(tuples(poly_b, T); corrections = (PolygonContents(), )))
         end
     end
     return polys
 end
+function _union(
+    target::TraitTarget{GI.MultiPolygonTrait}, ::Type{T},
+    ::GI.PolygonTrait, poly_a,
+    ::GI.MultiPolygonTrait, multipoly_b;
+    fix_multipoly = UnionIntersectingPolygons(), kwargs...,
+) where T
+    return GI.MultiPolygon(_union(
+        TraitTarget(GI.PolygonTrait()), T, 
+        GI.trait(poly_a), poly_a, 
+        GI.trait(multipoly_b), multipoly_b;
+        fix_multipoly, kwargs...,
+    ))
+end
 
 #= Multipolygon with polygon union is equivalent to taking the union of the poylgon with the
 multipolygon and thus simply switches the order of operations and calls the above method. =#
@@ -244,6 +271,13 @@ _union(
     kwargs...,
 ) where T = union(poly_b, multipoly_a; target, kwargs...)
 
+_union(
+    target::TraitTarget{GI.MultiPolygonTrait}, ::Type{T},
+    ::GI.MultiPolygonTrait, multipoly_a,
+    ::GI.PolygonTrait, poly_b;
+    kwargs...,
+) where T = union(poly_b, multipoly_a; target, kwargs...)
+
 #= Multipolygon with multipolygon union - note that all of the sub-polygons of `multipoly_a`
 and the sub-polygons of `multipoly_b` are included and combined together where there are
 intersections. Unless specified with `fix_multipoly = nothing`, `multipolygon_b` will be
@@ -267,6 +301,20 @@ function _union(
     return polys
 end
 
+function _union(
+    target::TraitTarget{GI.MultiPolygonTrait}, ::Type{T},
+    ::GI.MultiPolygonTrait, multipoly_a,
+    ::GI.MultiPolygonTrait, multipoly_b;
+    fix_multipoly = UnionIntersectingPolygons(), kwargs...,
+) where T
+    return GI.MultiPolygon(_union( # this is the method directly above, and returns a vector of polygons
+        TraitTarget{GI.PolygonTrait}(), T, 
+        GI.MultiPolygonTrait(), multipoly_a, 
+        GI.MultiPolygonTrait(), multipoly_b; 
+        fix_multipoly, kwargs...
+    ))
+end
+
 # Many type and target combos aren't implemented
 function _union(
     ::TraitTarget{Target}, ::Type{T},

src/transformations/correction/geometry_correction.jl

@@ -44,7 +44,18 @@ application_level(gc::GeometryCorrection) = error("Not implemented yet for $(gc)
 
 (gc::GeometryCorrection)(trait::GI.AbstractGeometryTrait, geometry) = error("Not implemented yet for $(gc) and $(trait).")
 
-function fix(geometry; corrections = GeometryCorrection[ClosedRing(),], kwargs...)
+"""
+    fix(x; corrections = GeometryCorrection[], kwargs...)
+
+Apply the given corrections to `x`, and return the corrected version.
+
+`x` may be a geometry, vector of geometries, feature collection, or table - 
+anything which [`apply`](@ref) will accept!
+
+Some available corrections are: [`ClosedRing`](@ref), [`PolygonContents`](@ref), [`UnionIntersectingPolygons`](@ref), [`DiffIntersectingPolygons`](@ref).
+
+"""
+function fix(geometry; corrections = GeometryCorrection[PolygonContents(), ClosedRing(),], kwargs...)
     traits = application_level.(corrections)
     final_geometry = geometry
     for Trait in (GI.PointTrait, GI.MultiPointTrait, GI.LineStringTrait, GI.LinearRingTrait, GI.MultiLineStringTrait, GI.PolygonTrait, GI.MultiPolygonTrait)

src/transformations/correction/polygon_contents.jl

@@ -0,0 +1,33 @@
+# # PolygonContents
+
+export PolygonContents
+
+#=
+Polygons should only contain linear rings.  This fix checks
+whether the contents of the polygon are linear rings or linestrings,
+and converts linestrings to linear rings.
+
+It does **NOT** check whether the linear rings are valid - it only checks
+the types of the polygon's constituent geometries.  You can use the [`ClosedRing`](@ref)
+geometry fix to check for validity after applying this fix.
+=#
+
+struct PolygonContents <: GeometryCorrection end
+
+application_level(::PolygonContents) = GI.PolygonTrait
+
+function (::PolygonContents)(::GI.PolygonTrait, polygon)
+    exterior = GI.getexterior(polygon)
+    fixed_exterior = _ls2lr(exterior)
+    holes = GI.gethole(polygon)
+    if isempty(holes)
+        return GI.Polygon([fixed_exterior])
+    end
+    fixed_holes = _ls2lr.(holes)
+    return GI.Polygon([fixed_exterior, fixed_holes...])
+end
+
+_ls2lr(x) = _ls2lr(GI.geomtrait(x), x)
+
+_ls2lr(::GI.LineStringTrait, x) = GI.LinearRing(GI.getpoint(x))
+_ls2lr(::GI.LinearRingTrait, x) = x