diff --git a/src/lookup.jl b/src/lookup.jl index 86751455..0adccfa6 100644 --- a/src/lookup.jl +++ b/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 """ diff --git a/src/methods/extract.jl b/src/methods/extract.jl index 4718069d..bf04ac94 100644 --- a/src/methods/extract.jl +++ b/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) diff --git a/src/methods/rasterize.jl b/src/methods/rasterize.jl index 1c1bd526..d459efeb 100644 --- a/src/methods/rasterize.jl +++ b/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`.")) diff --git a/src/openstack.jl b/src/openstack.jl index 5ba646e9..bb930b93 100644 --- a/src/openstack.jl +++ b/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] diff --git a/src/polygon_ops.jl b/src/polygon_ops.jl index cc71c598..b517191d 100644 --- a/src/polygon_ops.jl +++ b/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,6 +125,7 @@ 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) @@ -132,7 +133,7 @@ function _iyperm(dims::Tuple{<:Dimension,<:Dimension,<:Dimension}) 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")) diff --git a/test/methods.jl b/test/methods.jl index 0ea19bf7..996c1559 100644 --- a/test/methods.jl +++ b/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