better primitives docs and comments, and some tweaks

src/primitives.jl

@@ -12,14 +12,66 @@ $CALC_EXTENT_KEYWORD
 
 # This file mainly defines the [`apply`](@ref) function.
 
+#=
+## What is `apply`?
+
+`apply` applies some function to every geometry matching the `Target`
+GeoInterface trait, in some arbitrarily nested object made up of:
+- `AbstractArray`s (we also try to iterate other non-GeoInteface compatible object)
+- `FeatureCollectionTrait` objects
+- `FeatureTrait` objects
+- `AbstractGeometryTrait` objects
+
+`apply` recursively calls itself through these nested
+layers until it reaches objects with the `Target` GeoInterface trait. When found `apply` applies the function `f`, and stops.
+
+The outer recursive functions then progressively rebuild the object
+using GeoInterface objects matching the original traits.
+
+If `PointTrait` is found  but it is not the `Target`, an error is thrown.
+This likely means the object contains a different geometry trait to 
+the target, such as `MultiPointTrait` when `LineStringTrait` was specified.
+
+To handle this possibility it may be necessary to make `Target` a
+`Union` of traits found at the same level of nesting, and define methods
+of `f` to handle all cases.
+
+Be careful making a union across "levels" of nesting, e.g. 
+`Union{FeatureTrait,PolygonTrait}`, as `_apply` will just never reach 
+`PolygonTrait` when all the polygons are wrapped in a `FeatureTrait` object.
+
+## Embedding:
+
+`extent` and `crs` can be embedded in all geometries, features, and
+feature collections as part of `apply`. Geometries deeper than `Target`
+will of course not have new `extent` or `crs` embedded.
+
+- `calc_extent` signals to recalculate an `Extent` and embed it. 
+- `crs` will be embedded as-is
+
+## Threading
+
+Threading is used at the outermost level possible - over
+an array, feature collection, or e.g. a MultiPolygonTrait where
+each `PolygonTrait` sub-geometry may be calculated on a different thread.
+=#
+
 """
     apply(f, target::Type{<:AbstractTrait}, obj; kw...)
 
-Reconstruct a geometry or feature using the function `f` on the `target` trait.
+Reconstruct a geometry, feature, feature collection, or nested vectors of
+either using the function `f` on the `target` trait.
 
-`f(target_geom) => x` where `x` also has the `target` trait, or an equivalent.
+`f(target_geom) => x` where `x` also has the `target` trait, or a trait that can 
+be substituted. For example, swapping `PolgonTrait` to `MultiPointTrait` will fail
+if the outer object has `MultiPolygonTrait`, but should work if it has `FeatureTrait`.
 
-The result is an functionally similar geometry with values depending on `f`
+Objects "shallower" than the target trait are always completely rebuilt, like
+a `Vector` of `FeatureCollectionTrait` of `FeatureTrait` when the target
+has `PolygonTrait` and is held in the features. But "deeper" objects may remain 
+unchanged - such as points and linear rings if the target is the same `PolygonTrait`.
+
+The result is a functionally similar geometry with values depending on `f`
 
 $APPLY_KEYWORDS
 
@@ -39,70 +91,92 @@ end
 """
 apply(f, ::Type{Target}, geom; kw...) where Target = _apply(f, Target, geom; kw...)
 
+# Call _apply again with the trait of `geom`
 _apply(f, ::Type{Target}, geom; kw...)  where Target =
     _apply(f, Target, GI.trait(geom), geom; kw...)
+# There is no trait and this is an AbstractArray - so just iterate over it calling _apply on the contents
 function _apply(f, ::Type{Target}, ::Nothing, A::AbstractArray; threaded=false, kw...) where Target
+    # For an Array there is nothing else to do but map `_apply` over all values
+    # _maptasks may run this level threaded if `threaded==true`, 
+    # but deeper `_apply` called in the closure will not be threaded 
     _maptasks(eachindex(A); threaded) do i
-        _apply(f, Target, A[i]; kw...)
+        _apply(f, Target, A[i]; threaded=false, kw...)
     end
 end
-# Try to _apply over iterables
+# There is no trait and this is not an AbstractArray.
+# Try to call _apply over it. We can't use threading 
+# as we don't know if we can can index into it. So just `map`.
+# (TODO: maybe `collect` first if `threaded=true` so we can thread?)
 _apply(f, ::Type{Target}, ::Nothing, iterable; kw...) where Target =
     map(x -> _apply(f, Target, x; kw...), iterable)
-# Rewrap feature collections
+# Rewrap all FeatureCollectionTrait feature collections as GI.FeatureCollection
 function _apply(f, ::Type{Target}, ::GI.FeatureCollectionTrait, fc; 
     crs=GI.crs(fc), calc_extent=false, threaded=false
 ) where Target
+    # Run _apply on all `features` in the feature collection, possibly threaded
     features = _maptasks(1:GI.nfeature(fc); threaded) do i
         feature = GI.getfeature(fc, i)
-        _apply(f, Target, feature; crs, calc_extent)::GI.Feature
+        _apply(f, Target, feature; crs, calc_extent, threaded=false)::GI.Feature
     end
     if calc_extent
-        extent = reduce(features; init=GI.extent(first(features))) do acc, f
-            Extents.union(acc, Extents.extent(f))
-        end
+        # Calculate the extent of the features
+        extent = mapreduce(GI.extent, Extents.union, features)
+        # Return a FeatureCollection with features, crs and caculated extent
         return GI.FeatureCollection(features; crs, extent)
     else
+        # Return a FeatureCollection with features and crs
         return GI.FeatureCollection(features; crs)
     end
 end
-# Rewrap features
+# Rewrap all FeatureTrait features as GI.Feature, keeping the properties
 function _apply(f, ::Type{Target}, ::GI.FeatureTrait, feature; 
     crs=GI.crs(feature), calc_extent=false, threaded=false
 ) where Target
+    # Run _apply on the contained geometry 
+    geometry = _apply(f, Target, GI.geometry(feature); crs, calc_extent, threaded)
+    # Get the feature properties
     properties = GI.properties(feature)
-    geometry = _apply(f, Target, GI.geometry(feature); crs, calc_extent)
     if calc_extent
+        # 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 GI.Feature(geometry; properties, crs, extent)
     else
+        # Return a new Feature with the new geometry, but the oroginal properties and crs 
         return GI.Feature(geometry; properties, crs)
     end
 end
 # Reconstruct nested geometries
 function _apply(f, ::Type{Target}, trait, geom; 
     crs=GI.crs(geom), calc_extent=false, threaded=false
 )::(GI.geointerface_geomtype(trait)) where Target
-    # TODO handle zero length...
+    # Map `_apply` over all sub geometries of `geom`
+    # to create a new vector of geometries
+    # TODO handle zero length
     geoms = _maptasks(1:GI.ngeom(geom); threaded) do i
-        _apply(f, Target, GI.getgeom(geom, i); crs, calc_extent)
+        _apply(f, Target, GI.getgeom(geom, i); crs, calc_extent, threaded=false)
     end
     if calc_extent
-        extent = GI.extent(first(geoms))
-        for g in geoms
-            extent = Extents.union(extent, GI.extent(g))
-        end
+        # 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
         return rebuild(geom, geoms; crs, extent)
     else
+        # Return a new geometryof the same trait as `geom`, holding the new `geoms` with `crs`
         return rebuild(geom, geoms; crs)
     end
 end
-# Apply f to the target geometry
-_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`
+# Fail loudly if we hit PointTrait without running `f`
+# (after PointTrait there is no further to dig with `_apply`)
 _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
+# Finally, these short methods are the main purpose of `apply`.
+# The `Trait` is a subtype of the `Target` (or identical to it)
+# So the `Target` is found. We apply `f` to geom and return it to previous 
+# _apply calls to be wrapped with the outer geometries/feature/featurecollection/array.
+_apply(f, ::Type{Target}, ::Trait, geom; crs=GI.crs(geom), kw...) where {Target,Trait<:Target} = f(geom)
+# Define some specific cases of this match to avoid method ambiguity
 _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)
@@ -111,7 +185,7 @@ _apply(f, ::Type{GI.FeatureCollectionTrait}, ::GI.FeatureCollectionTrait, fc; kw
     unwrap(target::Type{<:AbstractTrait}, obj)
     unwrap(f, target::Type{<:AbstractTrait}, obj)
 
-Unwrap the geometry to vectors, down to the target trait.
+Unwrap the object newst to vectors, down to the target trait.
 
 If `f` is passed in it will be applied to the target geometries
 as they are found.
@@ -139,10 +213,14 @@ unwrap(f, target::Type{GI.FeatureTrait}, ::GI.FeatureTrait, feature) = f(feature
 unwrap(f, target::Type{GI.FeatureCollectionTrait}, ::GI.FeatureCollectionTrait, fc) = f(fc)
 
 """
-    flatten(target::Type{<:GI.AbstractTrait}, geom)
+    flatten(target::Type{<:GI.AbstractTrait}, obj)
+    flatten(f, target::Type{<:GI.AbstractTrait}, obj)
 
-Lazily flatten any geometry, feature or iterator of geometries or features
-so that objects with the specified trait are returned by the iterator.
+Lazily flatten any `AbstractArray`, iterator, `FeatureCollectionTrait`, 
+`FeatureTrait` or `AbstractGeometryTrait` object `obj`, so that objects 
+with the `target` trait are returned by the iterator.
+
+If `f` is passed in it will be applied to the target geometries.
 """
 flatten(::Type{Target}, geom) where {Target<:GI.AbstractTrait} = flatten(identity, Target, geom) 
 flatten(f, ::Type{Target}, geom) where {Target<:GI.AbstractTrait} = _flatten(f, Target, geom) 
@@ -181,7 +259,10 @@ All objects in `components` must have the same `GeoInterface.trait`.
 
 Ususally used in combination with `flatten`.
 """
-reconstruct(geom, components) = first(_reconstruct(geom, components))
+function reconstruct(geom, components)
+    obj, iter = _reconstruct(geom, components)
+    return obj
+end
 
 _reconstruct(geom, components) = 
     _reconstruct(typeof(GI.trait(first(components))), geom, components, 1) 
@@ -190,6 +271,7 @@ _reconstruct(::Type{Target}, geom, components, iter) where Target =
 # Try to reconstruct over iterables
 function _reconstruct(::Type{Target}, ::Nothing, iterable, components, iter) where Target 
     vect = map(iterable) do x
+        # iter is updated by _reconstruct here
         obj, iter = _reconstruct(Target, x, components, iter)
         obj
     end
@@ -198,6 +280,7 @@ end
 # Reconstruct feature collections
 function _reconstruct(::Type{Target}, ::GI.FeatureCollectionTrait, fc, components, iter) where Target
     features = map(GI.getfeature(fc)) do feature
+        # iter is updated by _reconstruct here
         newfeature, iter = _reconstruct(Target, feature, components, iter)
         newfeature
     end
@@ -209,6 +292,7 @@ function _reconstruct(::Type{Target}, ::GI.FeatureTrait, feature, components, it
 end
 function _reconstruct(::Type{Target}, trait, geom, components, iter) where Target
     geoms = map(GI.getgeom(geom)) do subgeom
+        # iter is updated by _reconstruct here
         subgeom1, iter = _reconstruct(Target, GI.trait(subgeom), subgeom, components, iter)
         subgeom1
     end
@@ -234,7 +318,7 @@ const BasicsGeoms = Union{GB.AbstractGeometry,GB.AbstractFace,GB.AbstractPoint,G
 
 Rebuild a geometry from child geometries.
 
-By default geometries will be rebuilt as a GeoInterface.Wrappers 
+By default geometries will be rebuilt as a `GeoInterface.Wrappers`
 geometry, but `rebuild` can have methods added to it to dispatch
 on geometries from other packages and specify how to rebuild them.
 
@@ -286,7 +370,7 @@ function _maptasks(f, taskrange; threaded=false)
         end
 
         # Finally we join the results into a new vector
-        return reduce(vcat, map(fetch, tasks))
+        return mapreduce(fetch, vcat, tasks)
     else
         return map(f, taskrange)
     end