Fix and test rasterize for feature collections (#651)

* fix rasterize points

* fix rasterize

* bugfix

* test

Author: rafaqz

src/methods/rasterize.jl

@@ -148,7 +148,7 @@ function Rasterizer(data::T; fill, geomcolumn=nothing, kw...) where T
         schema = Tables.schema(data)
         cols = Tables.columns(data)
         colnames = Tables.columnnames(Tables.columns(data))
-        fillitr = _iterable_fill(cols, fill)
+        fillitr = _iterable_fill(nothing, cols, fill)
         # If fill is a symbol or tuple of Symbol we need to allocate based on the column type
         geomcolname = if isnothing(geomcolumn)
             geomcols = GI.geometrycolumns(data)
@@ -171,27 +171,27 @@ function Rasterizer(data::T; fill, geomcolumn=nothing, kw...) where T
         Rasterizer(GeoInterface.trait(data), data; fill, kw...)
     end
 end
-function Rasterizer(::GI.AbstractFeatureCollectionTrait, fc; name, fill, kw...)
-    # TODO: how to handle when there are fillvals with different types
-    # fillval = _featurefillval(GI.getfeature(fc, 1), fill)
-    fillitr = _iterable_fill(fc, fill)
-    geometries = map(f -> GI.geometry(f), GI.getfeature(fc))
-    Rasterizer(geometries; kw..., fill, fillitr)
+function Rasterizer(trait::GI.AbstractFeatureCollectionTrait, fc; fill, kw...)
+    fillitr = _iterable_fill(trait, fc, fill)
+    geoms = map(f -> GI.geometry(f), GI.getfeature(fc))
+    Rasterizer(geoms, fill, fillitr; kw...)
 end
-function Rasterizer(::GI.AbstractFeatureTrait, feature; fill, kw...)
-    fillitr = _iterable_fill(feature, fill)
-    # fillval = _featurefillval(feature, fill)
-    Rasterizer(GI.geometry(feature), fill, fillitr; kw...)
+function Rasterizer(trait::GI.AbstractFeatureTrait, feature; fill, kw...)
+    fillitr = _iterable_fill(trait, feature, fill)
+    geom = GI.geometry(feature)
+    Rasterizer(geom, fill, fillitr; kw...)
 end
-function Rasterizer(::GI.GeometryCollectionTrait, collection; kw...)
-    Rasterizer(collect(GI.getgeom(collection)); kw...)
+function Rasterizer(trait::GI.GeometryCollectionTrait, collection; kw...)
+    geoms = collect(GI.getgeom(collection))
+    Rasterizer(geoms; kw...)
 end
-function Rasterizer(::Nothing, geoms; fill, kw...)
-    fillitr = _iterable_fill(geoms, fill)
+function Rasterizer(trait::Nothing, geoms; fill, kw...)
+    fillitr = _iterable_fill(trait, geoms, fill)
     Rasterizer(geoms, fill, fillitr; kw...)
 end
-function Rasterizer(::GI.AbstractGeometryTrait, geom; fill, kw...)
-    Rasterizer(geom, fill, _iterable_fill(geom, fill); kw...)
+function Rasterizer(trait::GI.AbstractGeometryTrait, geom; fill, kw...)
+    fillitr = _iterable_fill(trait, geom, fill)
+    Rasterizer(geom, fill, fillitr; kw...)
 end
 
 function get_eltype_missingval(eltype, missingval, fill, fillitr, init, filename, op, reducer)
@@ -274,51 +274,53 @@ function _filter_name(name, fill)
 end
 
 # A Tuple of `Symbol` is multiple keys to make a RasterStack
-_iterable_fill(data, keys::Tuple{Symbol,Vararg}) =
-    NamedTuple{keys}(map(k -> _iterable_fill(data, k), keys))
+_iterable_fill(trait, data, keys::Tuple{Symbol,Vararg}) =
+    NamedTuple{keys}(map(k -> _iterable_fill(trait, data, k), keys))
 # A Symbol is a Table or FeatureCollection key, it cant be used as fill itself
-function _iterable_fill(data, key::Symbol)
+function _iterable_fill(trait, data, key::Symbol)
     if GI.isfeature(data)
-        return get(x -> error("feature has no key $key"), GI.properties(data), key)
+        return get(() -> throw(ArgumentError("feature has no property `:$key`")), GI.properties(data), key)
+    elseif trait isa GI.FeatureCollectionTrait
+        return [get(() -> throw(ArgumentError("feature has no property `:$key`")), GI.properties(f), key) for f in GI.getfeature(data)]
     end
     cols = Tables.columns(data)
     # For column tables, get the column now
     names = Tables.columnnames(cols)
     key in names || _fill_key_error(names, key)
     return Tables.getcolumn(cols, key)
 end
-_iterable_fill(data, fill::Function) = fill
-_iterable_fill(data, fill::NamedTuple) = begin
-    map(f -> _iterable_fill(data, f), fill)
+_iterable_fill(trait, data, fill::Function) = fill
+_iterable_fill(trait, data, fill::NamedTuple) = begin
+    map(f -> _iterable_fill(trait, data, f), fill)
 end
-# Inspect our data and fill as much as possible to check they match
-# and cycle any fill of known length one
-function _iterable_fill(data, fill)
-    if GI.isgeometry(data) || GI.isfeature(data)
+function _iterable_fill(trait, data, fill)
+    # trait isa Union{GI.AbstractGeometryTrait,GI.FeatureTrait} && return fill
+    if trait isa GI.AbstractGeometryTrait || trait isa GI.FeatureTrait
         return fill
-    end
-    if fill isa Number
+    elseif fill isa Number 
         return Iterators.cycle(fill)
-    end
-    if Tables.istable(typeof(data))
+    elseif Tables.istable(typeof(data))
         # we don't need the keys, just the column length
         data = first(Tables.columns(data))
     end
-    if Base.IteratorSize(data) isa Union{Base.HasShape,Base.HasLength}
-        fillvec = collect(fill)
-        l = length(fillvec)
+
+    if trait isa GI.FeatureCollectionTrait
+        n = GI.nfeature(data)
+    elseif Base.IteratorSize(data) isa Union{Base.HasShape,Base.HasLength}
         n = length(data)
-        if l == 1
-            # Cycle all length one iterables to fill every row
-            return Iterators.cycle(fillvec[1])
-        elseif !(l == n)
-            throw(ArgumentError("Length of fill $l does not match length of iterator $n"))
-        else
-            return fillvec
-        end
     else
         return fill
     end
+    fillvec = collect(fill)
+    l = length(fillvec)
+    if l == 1
+        # Cycle all length one iterables to fill every row
+        return Iterators.cycle(fillvec[1])
+    elseif !(l == n)
+        throw(ArgumentError("Length of fill $l does not match length of iterator $n"))
+    else
+        return fillvec
+    end
 end
 
 _getfill(itrs::NamedTuple, i::Int) = map(itr -> _getfill(itr, i), itrs)
@@ -660,9 +662,10 @@ end
 _rasterize_points!(A, r::Rasterizer) = _rasterize_points!(A, r.geom, r.fillitr, r)
 _rasterize_points!(A, geom, fillitr, r::Rasterizer) =
     _rasterize_points!(A, GI.trait(geom), geom, fillitr, r)
-function _rasterize_points!(A, ::GI.AbstractGeometryTrait, geom, fill, r::Rasterizer)
+
+function _rasterize_points!(A, trait::GI.AbstractGeometryTrait, geom, fill, r::Rasterizer)
     points = GI.getpoint(geom)
-    fill1 =_iterable_fill(points, fill)
+    fill1 =_iterable_fill(nothing, points, fill)
     _rasterize_points!(A, nothing, points, fill1, r)
 end
 function _rasterize_points!(A, ::GI.GeometryCollectionTrait, collection, fill, r::Rasterizer)

test/rasterize.jl

@@ -1,5 +1,5 @@
 using Rasters, Test, ArchGDAL, ArchGDAL.GDAL, Dates, Statistics, DataFrames, Extents, Shapefile, GeometryBasics
-import GeoInterface
+import GeoInterface as GI
 using Rasters.Lookups, Rasters.Dimensions 
 using Rasters: bounds
 
@@ -25,11 +25,14 @@ multi_point = ArchGDAL.createmultipoint(pointvec)
 linestring = ArchGDAL.createlinestring(pointvec)
 multi_linestring = ArchGDAL.createmultilinestring([pointvec])
 linearring = ArchGDAL.createlinearring(pointvec)
-line_collection = GeoInterface.GeometryCollection([multi_linestring])
-poly_collection = GeoInterface.GeometryCollection([polygon])
-pointfc = map(GeoInterface.getpoint(polygon), vals) do geom, v
+line_collection = GI.GeometryCollection([multi_linestring])
+poly_collection = GI.GeometryCollection([polygon])
+pointtable = map(GI.getpoint(polygon), vals) do geom, v
     (geometry=geom, val1=v, val2=2.0f0v)
 end
+pointfc = map(GI.getpoint(polygon), vals) do geom, v
+    GI.Feature(geom; properties=(val1=v, val2=2.0f0v))
+end |> GI.FeatureCollection
 table = (X=first.(pointvec), Y=last.(pointvec), othercol=zero.(last.(pointvec)))
 
 test_shape_dir = realpath(joinpath(dirname(pathof(Shapefile)), "..", "test", "shapelib_testcases"))
@@ -48,10 +51,11 @@ st = RasterStack((A1, copy(A1)))
     geom = pointvec
     geom = line_collection
     geom = poly_collection
+    geom = pointfc
     threaded = true
 
     for A in (A1, A2), 
-        geom in (pointvec, pointfc, multi_point, linestring, multi_linestring, linearring, polygon, multi_polygon, table, line_collection, poly_collection),
+        geom in (pointvec, pointtable, pointfc, multi_point, linestring, multi_linestring, linearring, polygon, multi_polygon, table, line_collection, poly_collection),
         threaded in (true, false)
 
         fill!(A, 0)
@@ -64,15 +68,15 @@ st = RasterStack((A1, copy(A1)))
         @test sum(A) == 4.0
         @test sum(rasterize(last, geom; to=A, shape=:point, fill=1, missingval=0, threaded)) == 4.0
         fill!(A, 0)
-        if !Tables.istable(geom)
+        if !(Tables.istable(geom) || GI.isfeaturecollection(geom))
             rasterize!(count, A, [geom, geom]; shape=:point, threaded)
             @test sum(A) == 10.0
             fill!(A, 0)
         end
     end
     geom = multi_point
     for A in (A1, A2), 
-        geom in (table, pointvec, pointfc, multi_point, linestring, multi_linestring, linearring, polygon, multi_polygon, line_collection, poly_collection),
+        geom in (table, pointvec, pointtable, pointfc, multi_point, linestring, multi_linestring, linearring, polygon, multi_polygon, line_collection, poly_collection),
         threaded in (true, false)
 
         rasterize!(sum, st, geom; shape=:point, fill=(layer1=2, layer2=3), threaded)
@@ -94,7 +98,7 @@ st = RasterStack((A1, copy(A1)))
 
 
     for A in (A1, A2), 
-        geom in (table, pointvec, pointfc, multi_point, linestring, multi_linestring, linearring, polygon, multi_polygon),
+        geom in (table, pointvec, pointtable, pointfc, multi_point, linestring, multi_linestring, linearring, polygon, multi_polygon),
         threaded in (true, false)
 
         st[:layer1] .= 0; st[:layer2] .= 0
@@ -180,11 +184,11 @@ end
 
 @testset "from geometries, tables and features of points" begin
     A = A1
-    data = DataFrame(pointfc)
+    data = DataFrame(pointtable)
     data = multi_point
     data = pointfc
 
-    for data in (pointfc, DataFrame(pointfc), multi_point, pointvec, reverse(pointvec))
+    for data in (pointtable, pointfc, DataFrame(pointtable), multi_point, pointvec, reverse(pointvec))
         @test sum(skipmissing(rasterize(sum, data; to=A, fill=1))) == 5
         @testset "to and fill Keywords are required" begin
             @test_throws ArgumentError R = rasterize(data; fill=1) 
@@ -206,8 +210,7 @@ end
     end
 
     @testset "a single feature" begin
-        feature = pointfc[4]
-        GeoInterface.isfeature(feature)
+        feature = pointtable[4]
         @testset "NTuple of Symbol fill makes an stack" begin
             rst = rasterize(feature; to=A, fill=(:val1, :val2))
             @test parent(rst.val1) isa Array{Union{Missing,Int},2}
@@ -225,9 +228,9 @@ end
     end
 
     @testset "feature collection, table from fill of Symbol keys" begin
+        data = pointtable
         data = pointfc
-        Tables.istable(data)
-        for data in (pointfc, DataFrame(pointfc)), threaded in (true, false)
+        for data in (pointfc, pointtable, DataFrame(pointtable)), threaded in (true, false)
             @testset "NTuple of Symbol fill makes an stack" begin
                 rst = rasterize(sum, data; to=A, fill=(:val1, :val2), threaded)
                 @test parent(rst.val1) isa Array{Union{Missing,Int},2}
@@ -445,6 +448,17 @@ end
     end
 end
 
+@testset "Rasterizing feature collection of polygons" begin
+    p1 = GI.Polygon([[[-55965.680060140774, -31588.16072168928], [-55956.50771556479, -31478.09258677756], [-31577.548550575284, -6897.015828572996], [-15286.184961223798, -15386.952072224134], [-9074.387601621409, -27468.20712382156], [-8183.4538916097845, -31040.003969070774], [-27011.85123029944, -38229.02388009402], [-54954.72822634951, -32258.9734800704], [-55965.680060140774, -31588.16072168928]]])
+    p2 = GI.Polygon([[[-80000.0, -80000.0], [-80000.0, 80000.0], [-60000.0, 80000.0], [-50000.0, 40000.0], [-60000.0, -80000.0], [-80000.0, -80000.0]]])
+    fc = GI.FeatureCollection([GI.Feature(p1, properties = (;val=1)), GI.Feature(p2, properties = (;val=2))])
+
+    vecint = rasterize(sum, [p1, p2]; size=(100, 100), fill=[1, 2])
+    fcint = rasterize(sum, fc; size=(100, 100), fill=[1, 2])
+    fcval = rasterize(sum, fc; size=(100, 100), fill=:val)
+    @test all(vecint .=== fcval .=== fcint)
+end
+
 @testset "Rasterizing with extra dimensions" begin
     A3 = cat(A1, A1, A1; dims=Band)
     Rasters.rasterize!(last, A3, polygon; fill=true)