Merge branch 'main' into sg/area_dist

src/methods/equals.jl

@@ -39,9 +39,13 @@ implementation, since it's a lot less work!
 Note that while we need the same set of points and edges, they don't need to be
 provided in the same order for polygons. For for example, we need the same set
 points for two multipoints to be equal, but they don't have to be saved in the
-same order. This requires checking every point against every other point in the
-two geometries we are comparing. Additionally, geometries and multi-geometries
-can be equal if the multi-geometry only includes that single geometry.
+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.
+Additionally, geometries and multi-geometries can be equal if the multi-geometry
+only includes that single geometry.
 =#
 
 """
@@ -96,11 +100,23 @@ function equals(::GI.PointTrait, p1, ::GI.PointTrait, p2)
     return true
 end
 
+"""
+    equals(::GI.PointTrait, p1, ::GI.MultiPointTrait, mp2)::Bool
+
+A point and a multipoint are equal if the multipoint is composed of a single
+point that is equivalent to the given point.
+"""
 function equals(::GI.PointTrait, p1, ::GI.MultiPointTrait, mp2)
     GI.npoint(mp2) == 1 || return false
     return equals(p1, GI.getpoint(mp2, 1))
 end
 
+"""
+    equals(::GI.MultiPointTrait, mp1, ::GI.PointTrait, p2)::Bool
+
+A point and a multipoint are equal if the multipoint is composed of a single
+point that is equivalent to the given point.
+"""
 equals(trait1::GI.MultiPointTrait, mp1, trait2::GI.PointTrait, p2) =
     equals(trait2, p2, trait1, mp1)
 
@@ -125,58 +141,147 @@ function equals(::GI.MultiPointTrait, mp1, ::GI.MultiPointTrait, mp2)
 end
 
 """
-    equals(::T, l1, ::T, l2) where {T<:GI.AbstractCurveTrait} ::Bool
-
-Two curves are equal if they share the same set of points going around the
-curve. 
+    _equals_curves(c1, c2, closed_type1, closed_type2)::Bool
+
+Two curves are equal if they share the same set of point, representing the same
+geometry. Both curves must must be composed of the same set of points, however,
+they do not have to wind in the same direction, or start on the same point to be
+equivalent.
+Inputs:
+    c1 first geometry
+    c2 second geometry
+    closed_type1::Bool true if c1 is closed by definition (polygon, linear ring)
+    closed_type2::Bool true if c2 is closed by definition (polygon, linear ring)
 """
-function equals(
-    ::Union{GI.LineTrait, GI.LineStringTrait, GI.LinearRingTrait}, l1,
-    ::Union{GI.LineTrait, GI.LineStringTrait, GI.LinearRingTrait}, l2,
-)
-    # Check line lengths match
-    n1 = GI.npoint(l1)
-    n2 = GI.npoint(l2)
-    # TODO: do we need to account for repeated last point??
+function _equals_curves(c1, c2, closed_type1, closed_type2)
+    # Check if both curves are closed or not
+    n1 = GI.npoint(c1)
+    n2 = GI.npoint(c2)
+    c1_repeat_point = GI.getpoint(c1, 1) == GI.getpoint(c1, n1)
+    n2 = GI.npoint(c2)
+    c2_repeat_point = GI.getpoint(c2, 1) == GI.getpoint(c2, n2)
+    closed1 = closed_type1 || c1_repeat_point
+    closed2 = closed_type2 || c2_repeat_point
+    closed1 == closed2 || return false
+    # How many points in each curve
+    n1 -= c1_repeat_point ? 1 : 0
+    n2 -= c2_repeat_point ? 1 : 0
     n1 == n2 || return false
-
-    # Find first matching point if it exists
-    p1 = GI.getpoint(l1, 1)
-    offset = nothing
+    n1 == 0 && return true
+    # Find offset between curves
+    jstart = nothing
+    p1 = GI.getpoint(c1, 1)
     for i in 1:n2
-        if equals(p1, GI.getpoint(l2, i))
-            offset = i - 1
+        if equals(p1, GI.getpoint(c2, i))
+            jstart = i
             break
         end
     end
-    isnothing(offset) && return false
-
-    # Then check all points are the same wrapping around line
-    for i in 1:n1
-        pi = GI.getpoint(l1, i)
-        j = i + offset
-        j = j <= n1 ? j : (j - n1)
-        pj = GI.getpoint(l2, j)
-        equals(pi, pj) || return false
+    # no point matches the first point
+    isnothing(jstart) && return false
+    # found match for only point
+    n1 == 1 && return true
+    # if isn't closed and first or last point don't match, not same curve
+    !closed_type1 && (jstart != 1 && jstart != n1) && return false
+    # Check if curves are going in same direction
+    i = 2
+    j = jstart + 1
+    j -= j > n2 ? n2 : 0
+    same_direction = equals(GI.getpoint(c1, i), GI.getpoint(c2, j))
+    # if only 2 points, we have already compared both
+    n1 == 2 && return same_direction
+    # Check all remaining points are the same wrapping around line
+    jstep = same_direction ? 1 : -1
+    for i in 2:n1
+        ip = GI.getpoint(c1, i)
+        j = jstart + (i - 1) * jstep
+        j += (0 < j <= n2) ? 0 : (n2 * -jstep)
+        jp = GI.getpoint(c2, j)
+        equals(ip, jp) || return false
     end
     return true
 end
 
+"""
+    equals(
+        ::Union{GI.LineTrait, GI.LineStringTrait}, l1,
+        ::Union{GI.LineTrait, GI.LineStringTrait}, l2,
+    )::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.
+"""
+equals(
+    ::Union{GI.LineTrait, GI.LineStringTrait}, l1,
+    ::Union{GI.LineTrait, GI.LineStringTrait}, l2,
+) = _equals_curves(l1, l2, false, false)
+
+"""
+    equals(
+        ::Union{GI.LineTrait, GI.LineStringTrait}, l1,
+        ::GI.LinearRingTrait, l2,
+    )::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
+rings are, so the line must have a repeated last point to be equal
+"""
+equals(
+    ::Union{GI.LineTrait, GI.LineStringTrait}, l1,
+    ::GI.LinearRingTrait, l2,
+) = _equals_curves(l1, l2, false, true)
+
+"""
+    equals(
+        ::GI.LinearRingTrait, l1,
+        ::Union{GI.LineTrait, GI.LineStringTrait}, l2,
+    )::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
+rings are, so the line must have a repeated last point to be equal
+"""
+equals(
+    ::GI.LinearRingTrait, l1,
+    ::Union{GI.LineTrait, GI.LineStringTrait}, l2,
+) = _equals_curves(l1, l2, true, false)
+
+"""
+    equals(
+        ::GI.LinearRingTrait, l1,
+        ::GI.LinearRingTrait, l2,
+    )::Bool
+
+Two linear rings are equal if they share the same set of points going along the
+curve. Note that rings are closed by definition, so they can have, but don't
+need, a repeated last point to be equal.
+"""
+equals(
+    ::GI.LinearRingTrait, l1,
+    ::GI.LinearRingTrait, l2,
+) = _equals_curves(l1, l2, true, true)
+
 """
     equals(::GI.PolygonTrait, geom_a, ::GI.PolygonTrait, geom_b)::Bool
 
 Two polygons are equal if they share the same exterior edge and holes.
 """
 function equals(::GI.PolygonTrait, geom_a, ::GI.PolygonTrait, geom_b)
     # Check if exterior is equal
-    equals(GI.getexterior(geom_a), GI.getexterior(geom_b)) || return false
+    _equals_curves(
+        GI.getexterior(geom_a), GI.getexterior(geom_b),
+        true, true,  # linear rings are closed by definition
+    ) || return false
     # Check if number of holes are equal
     GI.nhole(geom_a) == GI.nhole(geom_b) || return false
     # Check if holes are equal
     for ihole in GI.gethole(geom_a)
         has_match = false
         for jhole in GI.gethole(geom_b)
-            if equals(ihole, jhole)
+            if _equals_curves(
+                ihole, jhole,
+                true, true,  # linear rings are closed by definition
+            )
                 has_match = true
                 break
             end
@@ -186,11 +291,23 @@ function equals(::GI.PolygonTrait, geom_a, ::GI.PolygonTrait, geom_b)
     return true
 end
 
+"""
+    equals(::GI.PolygonTrait, geom_a, ::GI.MultiPolygonTrait, geom_b)::Bool
+
+A polygon and a multipolygon are equal if the multipolygon is composed of a
+single polygon that is equivalent to the given polygon.
+"""
 function equals(::GI.PolygonTrait, geom_a, ::MultiPolygonTrait, geom_b)
     GI.npolygon(geom_b) == 1 || return false
     return equals(geom_a, GI.getpolygon(geom_b, 1))
 end
 
+"""
+    equals(::GI.MultiPolygonTrait, geom_a, ::GI.PolygonTrait, geom_b)::Bool
+
+A polygon and a multipolygon are equal if the multipolygon is composed of a
+single polygon that is equivalent to the given polygon.
+"""
 equals(trait_a::GI.MultiPolygonTrait, geom_a, trait_b::PolygonTrait, geom_b) = 
     equals(trait_b, geom_b, trait_a, geom_a)
 

src/primitives.jl

@@ -1,3 +1,4 @@
+
 # # Primitive functions
 
 # This file mainly defines the [`apply`](@ref) function.
@@ -27,15 +28,24 @@ apply(f, ::Type{Target}, geom; kw...) where Target = _apply(f, Target, geom; kw.
 
 _apply(f, ::Type{Target}, geom; kw...)  where Target =
     _apply(f, Target, GI.trait(geom), geom; kw...)
+function _apply(f, ::Type{Target}, ::Nothing, A::AbstractArray; threaded=false, kw...) where Target
+    _maptasks(eachindex(A); threaded) do i
+        _apply(f, Target, A[i]; kw...)
+    end
+end
 # Try to _apply over iterables
 _apply(f, ::Type{Target}, ::Nothing, iterable; kw...) where Target =
     map(x -> _apply(f, Target, x; kw...), iterable)
 # Rewrap feature collections
-function _apply(f, ::Type{Target}, ::GI.FeatureCollectionTrait, fc; crs=GI.crs(fc), calc_extent=false) where Target
-    applicator(feature) = _apply(f, Target, feature; crs, calc_extent)::GI.Feature
-    features = map(applicator, GI.getfeature(fc))
+function _apply(f, ::Type{Target}, ::GI.FeatureCollectionTrait, fc; 
+    crs=GI.crs(fc), calc_extent=false, threaded=false
+) where Target
+    features = _maptasks(1:GI.nfeature(fc); threaded) do i
+        feature = GI.getfeature(fc, i)
+        _apply(f, Target, feature; crs, calc_extent)::GI.Feature
+    end
     if calc_extent
-        extent = rebuce(features; init=GI.extent(first(features))) do (acc, f)
+        extent = reduce(features; init=GI.extent(first(features))) do acc, f
             Extents.union(acc, Extents.extent(f))
         end
         return GI.FeatureCollection(features; crs, extent)
@@ -44,7 +54,9 @@ function _apply(f, ::Type{Target}, ::GI.FeatureCollectionTrait, fc; crs=GI.crs(f
     end
 end
 # Rewrap features
-function _apply(f, ::Type{Target}, ::GI.FeatureTrait, feature; crs=GI.crs(feature), calc_extent=false) where Target
+function _apply(f, ::Type{Target}, ::GI.FeatureTrait, feature; 
+    crs=GI.crs(feature), calc_extent=false, threaded=false
+) where Target
     properties = GI.properties(feature)
     geometry = _apply(f, Target, GI.geometry(feature); crs, calc_extent)
     if calc_extent
@@ -56,11 +68,12 @@ function _apply(f, ::Type{Target}, ::GI.FeatureTrait, feature; crs=GI.crs(featur
 end
 # Reconstruct nested geometries
 function _apply(f, ::Type{Target}, trait, geom; 
-    crs=GI.crs(geom), calc_extent=false
+    crs=GI.crs(geom), calc_extent=false, threaded=false
 )::(GI.geointerface_geomtype(trait)) where Target
     # TODO handle zero length...
-    applicator(g) = _apply(f, Target, g; crs, calc_extent)
-    geoms = map(applicator, GI.getgeom(geom))
+    geoms = _maptasks(1:GI.ngeom(geom); threaded) do i
+        _apply(f, Target, GI.getgeom(geom, i); crs, calc_extent)
+    end
     if calc_extent
         extent = GI.extent(first(geoms))
         for g in geoms
@@ -72,14 +85,14 @@ function _apply(f, ::Type{Target}, trait, geom;
     end
 end
 # Apply f to the target geometry
-_apply(f, ::Type{Target}, ::Trait, geom; crs=GI.crs(geom), calc_extent=false) where {Target,Trait<:Target} = f(geom)
+_apply(f, ::Type{Target}, ::Trait, geom; crs=GI.crs(geom), kw...) where {Target,Trait<:Target} = f(geom)
 # Fail if we hit PointTrait without running `f`
-_apply(f, ::Type{Target}, trait::GI.PointTrait, geom; crs=nothing, calc_extent=false) where Target =
+_apply(f, ::Type{Target}, trait::GI.PointTrait, geom; crs=nothing, kw...) where Target =
     throw(ArgumentError("target $Target not found, but reached a `PointTrait` leaf"))
 # Specific cases to avoid method ambiguity
-_apply(f, ::Type{GI.PointTrait}, trait::GI.PointTrait, geom; crs=nothing, calc_extent=false) = f(geom)
-_apply(f, ::Type{GI.FeatureTrait}, ::GI.FeatureTrait, feature; crs=GI.crs(feature), calc_extent=false) = f(feature)
-_apply(f, ::Type{GI.FeatureCollectionTrait}, ::GI.FeatureCollectionTrait, fc; crs=GI.crs(fc)) = f(fc)
+_apply(f, ::Type{GI.PointTrait}, trait::GI.PointTrait, geom; kw...) = f(geom)
+_apply(f, ::Type{GI.FeatureTrait}, ::GI.FeatureTrait, feature; kw...) = f(feature)
+_apply(f, ::Type{GI.FeatureCollectionTrait}, ::GI.FeatureCollectionTrait, fc; kw...) = f(fc)
 
 """
     unwrap(target::Type{<:AbstractTrait}, obj)
@@ -234,3 +247,34 @@ end
 function rebuild(trait::GI.PolygonTrait, geom::GB.Polygon, child_geoms; crs=nothing)
     Polygon(child_geoms[1], child_geoms[2:end])
 end
+
+using Base.Threads: nthreads, @threads, @spawn
+
+
+# Threading utility, modified Mason Protters threading PSA
+# run `f` over ntasks, where f recieves an AbstractArray/range
+# of linear indices
+function _maptasks(f, taskrange; threaded=false)
+    if threaded
+        ntasks = length(taskrange)
+        # Customize this as needed. 
+        # More tasks have more overhead, but better load balancing
+        tasks_per_thread = 2 
+        chunk_size = max(1, ntasks ÷ (tasks_per_thread * nthreads()))
+        # partition the range into chunks
+        task_chunks = Iterators.partition(taskrange, chunk_size) 
+        # Map over the chunks
+        tasks = map(task_chunks) do chunk
+            # Spawn a task to process this chunk
+            @spawn begin
+                # Where we map `f` over the chunk indices 
+                map(f, chunk)
+            end
+        end
+
+        # Finally we join the results into a new vector
+        return reduce(vcat, map(fetch, tasks))
+    else
+        return map(f, taskrange)
+    end
+end

src/transformations/extent.jl

@@ -6,7 +6,7 @@ calculating and adding an `Extents.Extent` to all objects.
 
 This can improve performance when extents need to be checked multiple times.
 """
-embed_extent(x) = apply(extent_applicator, AbstractTrait, x)
+embed_extent(x; kw...) = apply(AbstractTrait, x; kw...)
 
 extent_applicator(x) = extent_applicator(trait(x), x)
 extent_applicator(::Nothing, xs::AbstractArray) = embed_extent.(xs)