fix rasterization for a single feature (#323)

* fix rasterization for a single feature

src/lookup.jl

@@ -167,7 +167,7 @@ LA.transformfunc(lookup::AffineProjected) = CoordinateTransformations.inv(lookup
 # DD.bounds(lookup::AffineProjected) = lookup.metadata
 
 function Dimensions.sliceunalligneddims(
- _, I::NTuple{<:Any,<:Union{Colon,AbstractArray}},
+ _, I::NTuple{2,<:Union{Colon,AbstractArray}},
  ud1::Dimension{<:AffineProjected}, ud2::Dimension{<:AffineProjected}
 )
  # swap colons for the dimension index, which is the same as the array axis
@@ -193,7 +193,7 @@ end
 
 function Base.reverse(lookup::AffineProjected)
  sp = reverse(span(lookup))
- rebuild(lookup; data=i, order=o, span=sp)
+ rebuild(lookup; order=o, span=sp)
 end
 
 """

src/methods/extract.jl

@@ -87,8 +87,9 @@ function _extract(x::RasterStackOrArray, ::GI.PointTrait, point; dims, names, at
 
  # Extract the values
  if any(map(ismissing, coords)) 
+ # TODO test this branch somehow
  geometry = map(_ -> missing, coords)
- layer_vals = map(_ -> missing, layer_keys)
+ layer_vals = map(_ -> missing, names)
  else
  selectors = map(dims, coords) do d, c
  _at_or_contains(d, c, atol)

src/methods/rasterize.jl

@@ -87,12 +87,15 @@ function _rasterize(to::DimTuple, data;
 )
  _rasterize(to, GeoInterface.trait(data), data; fill, name, kw...)
 end
-function _rasterize(to::DimTuple, ::GI.AbstractFeatureTrait, feature; fill, kw...)
- fill = _featurefillval(feature, fill)
- dest = _create_rasterize_dest(fill, to; kw...)
+function _rasterize(to::DimTuple, ::GI.AbstractFeatureTrait, feature; fill, name, kw...)
+ fillval = _featurefillval(feature, fill)
+ name = _filter_name(name, fill)
+ @show name fillval
+ dest = _create_rasterize_dest(fillval, to; name, kw...)
  return rasterize!(dest, feature; fill, kw...)
 end
 function _rasterize(to::DimTuple, ::GI.AbstractFeatureCollectionTrait, fc; name, fill, kw...)
+ # TODO: how to handle when there are fillvals with different types
  fillval = _featurefillval(GI.getfeature(fc, 1), fill)
  name = _filter_name(name, fill)
  dest = _create_rasterize_dest(fillval, to; name, kw...)
@@ -122,6 +125,7 @@ function _rasterize(to::DimTuple, ::Nothing, data; fill, name, kw...)
  return rasterize!(dest, data; fill, kw...)
 end
 
+# Create a destination raster to fill into using `rasterize!`
 function _create_rasterize_dest(fill, dims; name=nothing, kw...)
  _create_rasterize_dest(fill, name, dims; kw...)
 end
@@ -249,7 +253,8 @@ function _rasterize!(x, ::GI.AbstractFeatureCollectionTrait, fc; fill, kw...)
  return x
 end
 function _rasterize!(x, ::GI.AbstractFeatureTrait, feature; fill, kw...)
- rasterize!(x, geom; fill=_featurefillval(feature, fill), kw...)
+ # TODO test this branch
+ rasterize!(x, GI.geometry(feature); fill=_featurefillval(feature, fill), kw...)
 end
 function _rasterize!(x, ::GI.AbstractGeometryTrait, geom; fill, _buffer=nothing, kw...)
  # TODO fix DimensionalData selectors so this works without _pad
@@ -397,6 +402,7 @@ end
 _filter_name(name, fill::NamedTuple) = keys(fill)
 _filter_name(name::NamedTuple, fill::NamedTuple) = keys(fill)
 _filter_name(name::Nothing, fill::Nothing) = nothing
+_filter_name(name::DimensionalData.NoName, fill::Nothing) = nothing
 _filter_name(name::Union{NamedTuple,Tuple,Array}, fill::NTuple{<:Any,Symbol}) = fill
 function _filter_name(name::Union{NamedTuple,Tuple,Array}, fill::Union{Tuple,Array})
  length(name) == length(fill) || throw(ArgumentError("`name` keyword (possibly from `to` object) does not match length of fill. A fix is to use a `NamedTuple` for `fill`."))

src/openstack.jl

@@ -23,6 +23,7 @@ ConstructionBase.constructorof(::Type{<:OpenStack{X,K}}) where {X,K} = OpenStack
 
 dataset(os::OpenStack) = os.dataset
 Base.keys(os::OpenStack{<:Any,K}) where K = K
-Base.values(os::OpenStack{<:Any}) = (fs[k] for k in keys(fs))
+# TODO test this, and does it make sense to return an iterator here?
+Base.values(os::OpenStack{<:Any}) = (os[k] for k in keys(os))
 # Indexing OpenStack returns memory-backed Raster. 
 Base.getindex(os::OpenStack, key::Symbol) = dataset(os)[key]

src/polygon_ops.jl

@@ -9,8 +9,8 @@ const DEFAULT_TABLE_DIM_KEYS = (:X, :Y, :Z)
 function _fill_geometry!(B::AbstractRaster, geom; kw...)
  _fill_geometry!(B, GI.trait(geom), geom; kw...)
 end
-function _fill_geometry!(B::AbstractRaster, ::GI.AbstractFeatureTrait, geom; kw...)
- return _fill_geometry!(B, GI.geometry(geom), geom; kw...)
+function _fill_geometry!(B::AbstractRaster, ::GI.AbstractFeatureTrait, feature; kw...)
+ return _fill_geometry!(B, GI.geometry(feature); kw...)
 end
 function _fill_geometry!(B::AbstractRaster, ::GI.AbstractFeatureCollectionTrait, fc; kw...)
  for feature in GI.getfeature(fc)
@@ -25,7 +25,7 @@ function _fill_geometry!(B::AbstractRaster, ::GI.AbstractGeometryTrait, geom; sh
  _fill_linestring!(B, geom; shape, kw...)
  elseif shape === :polygon
  geomextent = _extent(geom)
- arrayextent = Extents.extent(B, DEFAULT_POINT_ORDER) 
+ arrayextent = Extents.extent(B, DEFAULT_POINT_ORDER)
  # Only fill if the gemoetry bounding box overlaps the array bounding box
  Extents.intersects(geomextent, arrayextent) || return B
  _fill_polygon!(B, geom; shape, geomextent, kw...)
@@ -125,14 +125,15 @@ function _iyperm(dims::Tuple{<:Dimension,<:Dimension})
  return iyperm
 end
 function _iyperm(dims::Tuple{<:Dimension,<:Dimension,<:Dimension})
+ # TODO: test this 3d case
  a1, a2, a3 = map(dims) do d
  l = parent(d)
  LA.ordered_firstindex(l):_order_step(l):LA.ordered_lastindex(l)
  end
  iyperm = Array{Int}(undef, length(a1) * length(a2) * length(a3))
  lis = (LinearIndices(size(dims))[i, j, k] for k in a3 for j in a2 for i in a1)
  for (i, li) in enumerate(lis)
- Iyperm[i] = li
+ iyperm[i] = li
  end
  return iyperm
 end
@@ -263,7 +264,7 @@ function _fill_line!(A::AbstractRaster, line, fill)
  # delta: How far to move along the ray to move 1 grid cell.
  cs = cos(angle)
  si = sin(angle)
- max_x, delta_x = if isapprox(cs, zero(cs); atol=1e-10) 
+ max_x, delta_x = if isapprox(cs, zero(cs); atol=1e-10)
  -Inf, Inf
  else
  1.0 / cs, xoffset / cs
@@ -280,7 +281,7 @@ function _fill_line!(A::AbstractRaster, line, fill)
  for t in 0:manhattan_distance
  D = map((d, o) -> d(o), dimconstructors, (x, y))
  if checkbounds(Bool, A, D...)
- if fill isa Function 
+ if fill isa Function
  @inbounds A[D...] = fill(A[D...])
  else
  @inbounds A[D...] = fill
@@ -333,7 +334,7 @@ function to_edges_and_nodes!(edges, nodes, lastnode, geom)
  for (n, point) in enumerate(GI.getpoint(geom))
  i = lastnode + n
  if n == npoints
- # The closing edge of a sub-polygon 
+ # The closing edge of a sub-polygon
  edges[i, 1] = i
  edges[i, 2] = lastnode + 1
  else
@@ -367,10 +368,11 @@ function _extent(::Nothing, data)
  Extents.union(ext, _extent(geom))
  end
  else
+ # TODO: test this branch
  # Table of points with dimension columns
  reduce(DEFAULT_TABLE_DIM_KEYS; init=(;)) do acc, key
- if key in Tables.columnnames(cols) 
- merge(acc, (; key=extrema(columns)))
+ if key in Tables.columnnames(cols)
+ merge(acc, (; key=extrema(cols[key])))
  else
  acc
  end
@@ -393,7 +395,7 @@ function _extent(::GI.AbstractTrait, geom)
  end
 end
 _extent(::GI.AbstractFeatureTrait, feature) = _extent(GI.geometry(feature))
-function _extent(::GI.AbstractFeatureCollectionTrait, features) 
+function _extent(::GI.AbstractFeatureCollectionTrait, features)
  features = GI.getfeature(features)
  reduce(features; init=_extent(first(features))) do acc, f
  Extents.union(acc, _extent(f))
@@ -403,12 +405,13 @@ end
 # extent_may_intersect
 # Check if there is an extent for the geometry
 function extent_may_intersect(x, geom)
+ # TODO: this is not actually used.
  rasterextent = Extents.extent(x, DEFAULT_POINT_ORDER)
  geomextent = GI.extent(geom)
- if ext isa Nothing 
+ if isnothing(rasterextent) || isnothing(geomextent)
  return true
  else
- return Extents.intersects(geomextent, rasterext)
+ return Extents.intersects(geomextent, rasterextent)
  end
 end
 
@@ -422,7 +425,7 @@ _dimcoord(::ZDim, point) = GI.z(point)
 # Get the shape category for a geometry
 @inline _geom_shape(geom) = _geom_shape(GI.geomtrait(geom))
 @inline _geom_shape(geom::Union{<:GI.PointTrait,<:GI.MultiPointTrait}) = :point
-@inline _geom_shape(geom::Union{<:GI.LineStringTrait,<:GI.MultiLineStringTrait}) = :line 
+@inline _geom_shape(geom::Union{<:GI.LineStringTrait,<:GI.MultiLineStringTrait}) = :line
 @inline _geom_shape(geom::Union{<:GI.LinearRingTrait,<:GI.PolygonTrait,<:GI.MultiPolygonTrait}) = :polygon
 # _geom_shape(trait, geom) = throw(ArgumentError("Geometry trait $trait not handled by Rasters.jl"))
 

test/methods.jl

@@ -12,7 +12,6 @@ gaNaN = replace_missing(ga, NaN32)
 gaMi = replace_missing(ga)
 st = RasterStack((a=A, b=B), (X, Y); missingval=(a=missing,b=missing))
 
-
 pointvec = [(-20.0, 30.0),
  (-20.0, 10.0),
  (0.0, 10.0),
@@ -403,7 +402,22 @@ end
  @test sum(ra) == 12
  end
  end
-
+ @testset "a single feature" begin
+ feature = pointfc[4]
+ GeoInterface.isfeature(feature)
+ @testset "NTuple of Symbol fill makes an stack" begin
+ rst = rasterize(feature; to=A, fill=(:val1, :val2))
+ rst = rasterize(feature; to=A, fill=(:val1, :val2))
+ @test keys(rst) == (:val1, :val2)
+ @test dims(rst) == dims(A)
+ @test map(sum ∘ skipmissing, rst) === (val1=4, val2=8.0f0)
+ end
+ @testset "Symbol fill makes an array" begin
+ ra = rasterize(feature; to=A, fill=:val1)
+ @test ra isa Raster{Union{Missing,Int64}}
+ @test name(ra) == :val1
+ end
+ end
  @testset "feature collection, table from fill of Symbol keys" begin
  for data in (pointfc, DataFrame(pointfc))
  @testset "NTuple of Symbol fill makes an stack" begin