Breaking: Ambiguities and DD 0.29

Project.toml

@@ -25,6 +25,7 @@ ProgressMeter = "92933f4c-e287-5a05-a399-4b506db050ca"
 RecipesBase = "3cdcf5f2-1ef4-517c-9805-6587b60abb01"
 Reexport = "189a3867-3050-52da-a836-e630ba90ab69"
 Setfield = "efcf1570-3423-57d1-acb7-fd33fddbac46"
+Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"
 
 [weakdeps]
 ArchGDAL = "c9ce4bd3-c3d5-55b8-8973-c0e20141b8c3"
@@ -58,7 +59,7 @@ CommonDataModel = "0.2.3, 0.3"
 ConstructionBase = "1"
 CoordinateTransformations = "0.6.2"
 DataFrames = "1"
-DimensionalData = "0.28.2"
+DimensionalData = "0.28.2, 0.29"
 DiskArrays = "0.3, 0.4"
 Extents = "0.1"
 FillArrays = "0.12, 0.13, 1"
@@ -106,8 +107,7 @@ SafeTestsets = "1bc83da4-3b8d-516f-aca4-4fe02f6d838f"
 Shapefile = "8e980c4a-a4fe-5da2-b3a7-4b4b0353a2f4"
 StableRNGs = "860ef19b-820b-49d6-a774-d7a799459cd3"
 StatsBase = "2913bbd2-ae8a-5f71-8c99-4fb6c76f3a91"
-Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 
 [targets]
-test = ["Aqua", "ArchGDAL", "CFTime", "CoordinateTransformations", "DataFrames", "GeoDataFrames", "GeometryBasics", "GRIBDatasets", "NCDatasets", "Plots", "Proj", "RasterDataSources", "SafeTestsets", "Shapefile", "StableRNGs", "Statistics", "StatsBase", "Test", "ZarrDatasets"]
+test = ["Aqua", "ArchGDAL", "CFTime", "CoordinateTransformations", "DataFrames", "GeoDataFrames", "GeometryBasics", "GRIBDatasets", "NCDatasets", "Plots", "Proj", "RasterDataSources", "SafeTestsets", "Shapefile", "StableRNGs", "StatsBase", "Test", "ZarrDatasets"]

src/Rasters.jl

@@ -28,19 +28,21 @@ import Adapt,
        Mmap,
        RecipesBase,
        Reexport,
-       Setfield
+       Setfield,
+       Statistics
 
 Reexport.@reexport using DimensionalData, GeoFormatTypes
 
 using DimensionalData.Tables,
       DimensionalData.Lookups,
-      DimensionalData.Dimensions
+      DimensionalData.Dimensions,
       DimensionalData.Lookups.IntervalSets
 
 using DimensionalData: Name, NoName
 using .Dimensions: StandardIndices, DimTuple
 using .Lookups: LookupTuple
 
+using Statistics: mean
 using RecipesBase: @recipe, @series
 using Base: tail, @propagate_inbounds
 

src/array.jl

@@ -105,20 +105,17 @@ function DD.modify(f, A::AbstractRaster)
     return rebuild(A, newdata)
 end
 
-function DD.DimTable(As::Tuple{<:AbstractRaster,Vararg{AbstractRaster}}...)
-    DD.DimTable(DimStack(map(read, As...)))
-end
+DD.DimTable(As::Tuple{<:AbstractRaster,Vararg{AbstractRaster}}) =
+    DD.DimTable(DimStack(map(read, As)))
 
 # DiskArrays methods
 
 DiskArrays.eachchunk(A::AbstractRaster) = DiskArrays.eachchunk(parent(A))
 DiskArrays.haschunks(A::AbstractRaster) = DiskArrays.haschunks(parent(A))
-function DA.readblock!(A::AbstractRaster, dst, r::AbstractUnitRange...)
+DA.readblock!(A::AbstractRaster, dst, r::AbstractUnitRange...) =
     DA.readblock!(parent(A), dst, r...)
-end
-function DA.writeblock!(A::AbstractRaster, src, r::AbstractUnitRange...)
+DA.writeblock!(A::AbstractRaster, src, r::AbstractUnitRange...) =
     DA.writeblock!(parent(A), src, r...)
-end
 
 # Base methods
 

src/lookup.jl

@@ -89,18 +89,30 @@ GeoInterface.crs(lookup::Projected) = lookup.crs
 mappedcrs(lookup::Projected) = lookup.mappedcrs
 dim(lookup::Projected) = lookup.dim
 
-function LA.selectindices(l::Projected, sel::Contains)
+@inline function LA.selectindices(l::Projected, sel::LA.Selector; kw...)
     selval = reproject(mappedcrs(l), crs(l), dim(l), val(sel))
-    LA.contains(l, rebuild(sel; val=selval))
+    LA._selectindices(l, rebuild(sel; val=selval); kw...)
 end
-function LA.selectindices(l::Projected, sel::At)
+@inline function LA.selectindices(l::Projected, sel::LA.Selector{<:AbstractVector}; kw...)
     selval = reproject(mappedcrs(l), crs(l), dim(l), val(sel))
-    LA.at(l, rebuild(sel; val=selval))
+    _selectvec(l, rebuild(sel; val=selval)sel; kw...)
 end
-function LA.selectindices(l::Projected, sel::Between)
+@inline function LA.selectindices(l::Projected, sel::LA.IntSelector{<:Tuple}; kw...)
+    selval = reproject(mappedcrs(l), crs(l), dim(l), val(sel))
+    _selecttuple(l, rebuild(sel; val=selval); kw...)
+end
+@inline LA.selectindices(l::Projected{<:Tuple}, sel::LA.IntSelector{<:Tuple}; kw...) = LA._selectindices(l, sel; kw...)
+@inline LA.selectindices(l::Projected{<:Tuple}, sel::LA.IntSelector{<:Tuple{<:Tuple,<:Tuple}}; kw...) = 
+    _selecttuple(l, sel; kw...)
+
+function LA.selectindices(l::Projected, sel::Between{<:Tuple})
     selval = map(v -> reproject(mappedcrs(l), crs(l), dim(l), v), val(sel))
     LA.between(l, rebuild(sel; val=selval))
 end
+function LA.selectindices(l::Projected, sel::T) where T<:DD.IntervalSets.Interval
+    left, right = map(v -> reproject(mappedcrs(l), crs(l), dim(l), v), (sel.left, sel.right))
+    LA.between(l, T(left, right))
+end
 
 """
     Mapped <: AbstractProjected

src/methods/burning/utils.jl

@@ -4,16 +4,16 @@ _geomindices(::GI.FeatureCollectionTrait, geoms) = 1:GI.nfeature(geoms)
 _geomindices(::GI.FeatureTrait, geoms) = _geomindices(GI.geometry(geoms))
 _geomindices(::GI.AbstractGeometryTrait, geoms) = 1:GI.ngeom(geoms)
 
-_getgeom(geoms, i::Integer) = _getgeom(GI.trait(geoms), geoms, i)
-_getgeom(::GI.FeatureCollectionTrait, geoms, i::Integer) = GI.geometry(GI.getfeature(geoms, i))
-_getgeom(::GI.FeatureTrait, geoms, i::Integer) = GI.getgeom(GI.geometry(geoms), i)
-_getgeom(::GI.AbstractGeometryTrait, geoms, i::Integer) = GI.getgeom(geom, i)
-_getgeom(::GI.PointTrait, geom, i::Integer) = error("PointTrait should not be reached")
-_getgeom(::Nothing, geoms, i::Integer) = geoms[i] # Otherwise we can probably just index?
+_getgeom(geoms, i::Integer) = __getgeom(GI.trait(geoms), geoms, i)
+__getgeom(::GI.FeatureCollectionTrait, geoms, i::Integer) = GI.geometry(GI.getfeature(geoms, i))
+__getgeom(::GI.FeatureTrait, geoms, i::Integer) = GI.getgeom(GI.geometry(geoms), i)
+__getgeom(::GI.AbstractGeometryTrait, geoms, i::Integer) = GI.getgeom(geom, i)
+__getgeom(::GI.PointTrait, geom, i::Integer) = error("PointTrait should not be reached")
+__getgeom(::Nothing, geoms, i::Integer) = geoms[i] # Otherwise we can probably just index?
 
-_getgeom(geoms) = _getgeom(GI.trait(geoms), geoms)
-_getgeom(::GI.FeatureCollectionTrait, geoms) = (GI.geometry(f) for f in GI.getfeature(geoms))
-_getgeom(::GI.FeatureTrait, geoms) = GI.getgeom(GI.geometry(geoms))
-_getgeom(::GI.AbstractGeometryTrait, geoms) = GI.getgeom(geoms)
-_getgeom(::GI.PointTrait, geom) = error("PointTrait should not be reached")
-_getgeom(::Nothing, geoms) = geoms
+_getgeom(geoms) = __getgeom(GI.trait(geoms), geoms)
+__getgeom(::GI.FeatureCollectionTrait, geoms) = (GI.geometry(f) for f in GI.getfeature(geoms))
+__getgeom(::GI.FeatureTrait, geoms) = GI.getgeom(GI.geometry(geoms))
+__getgeom(::GI.AbstractGeometryTrait, geoms) = GI.getgeom(geoms)
+__getgeom(::GI.PointTrait, geom) = error("PointTrait should not be reached")
+__getgeom(::Nothing, geoms) = geoms

src/methods/mask.jl

@@ -170,7 +170,7 @@ function _mask!(x::RasterStackOrArray, geom; kw...)
 end
 # Array mask
 function _mask!(st::RasterStack, with::AbstractRaster; kw...)
-    map(A -> mask!(A; with, kw...), st)
+    maplayers(A -> mask!(A; with, kw...), st)
     return st
 end
 
@@ -431,7 +431,9 @@ end
 
 _false_to_missing(b::Bool) = (b ? true : missing)::Union{Missing,Bool}
 
-function _mask_multilayer(layers::Union{<:AbstractRasterStack,<:AbstractRasterSeries}, to; 
+_mask_multilayer(st::AbstractRasterStack, to; kw...) =
+    _mask_multilayer(layers(st), to; kw...)
+function _mask_multilayer(layers::Union{<:NamedTuple,<:AbstractRasterSeries}, to; 
     _dest_presentval,
     _dest_missingval, 
     missingval=nothing, 

src/methods/mosaic.jl

@@ -135,8 +135,9 @@ nothing
 
 $EXPERIMENTAL
 """
-mosaic!(f::Function, x::RasterStackOrArray, regions::RasterStackOrArray...; kw...) = mosaic!(f, x, regions; kw...)
-function mosaic!(f::Function, A::AbstractRaster{T}, regions;
+mosaic!(f::Function, x::RasterStackOrArray, regions::RasterStackOrArray...; kw...) =
+    _mosaic!(f, x, regions; kw...)
+function _mosaic!(f::Function, A::AbstractRaster{T}, regions::Tuple;
     missingval=missingval(A), atol=maybe_eps(T)
 ) where T
     _without_mapped_crs(A) do A1
@@ -169,7 +170,7 @@ function mosaic!(f::Function, A::AbstractRaster{T}, regions;
     end
     return A
 end
-function mosaic!(f::Function, st::AbstractRasterStack, regions::Tuple; kw...)
+function _mosaic!(f::Function, st::AbstractRasterStack, regions::Tuple; kw...)
     map(st, regions...) do A, r...
         mosaic!(f, A, r; kw...)
     end

src/methods/rasterize.jl

@@ -420,7 +420,7 @@ function rasterize(reducer::typeof(count), data; fill=nothing, init=nothing, kw.
 end
 # `mean` is sum ./ count. This is actually optimal with threading, 
 # as its means order is irrelevant so its threadsafe.
-function rasterize(reducer::typeof(DD.Statistics.mean), data; fill, kw...)
+function rasterize(reducer::typeof(Statistics.mean), data; fill, kw...)
     sums = rasterize(sum, data; kw..., fill)
     counts = rasterize(count, data; kw..., fill=nothing)
     rebuild(sums ./ counts; name=:mean)

src/methods/trim.jl

@@ -131,5 +131,5 @@ end
 
 # Broadcast over the array and tracker to mark axis indices as being missing or not
 _update!(tr::AxisTrackers, A::AbstractRaster) = tr .= A .!== missingval(A)
-_update!(tr::AxisTrackers, st::AbstractRasterStack) = map(A -> tr .= A .!== missingval(A), st)
+_update!(tr::AxisTrackers, st::AbstractRasterStack) = maplayers(A -> tr .= A .!== missingval(A), st)
 

src/methods/zonal.jl

@@ -84,7 +84,7 @@ _zonal(f, x::Raster, of::Extents.Extent; skipmissing=true) =
 function _zonal(f, x::RasterStack, ext::Extents.Extent; skipmissing=true)
     cropped = crop(x; to=ext, touches=true)
     prod(size(cropped)) > 0 || return missing
-    return map(cropped) do A
+    return maplayers(cropped) do A
         _maybe_skipmissing_call(f, A, skipmissing)
     end
 end
@@ -107,9 +107,9 @@ function _zonal(f, st::AbstractRasterStack, ::GI.AbstractGeometryTrait, geom;
     skipmissing=true, kw...
 )
     cropped = crop(st; to=geom, touches=true)
-    prod(size(cropped)) > 0 || return map(_ -> missing, st)
+    prod(size(cropped)) > 0 || return map(_ -> missing, layerdims(st))
     masked = mask(cropped; with=geom, kw...)
-    return map(masked) do A
+    return maplayers(masked) do A
         prod(size(A)) > 0 || return missing
         _maybe_skipmissing_call(f, A, skipmissing)
     end

src/read.jl

@@ -50,6 +50,10 @@ function Base.read!(src::AbstractRasterSeries, dst::AbstractRasterSeries)
     map(read!, src, dst)
     return dst
 end
+Base.read!(::AbstractRaster{<:Union{AbstractString,NamedTuple},1}, ::AbstractRasterSeries) =
+    error("Cant read Raster to RasterSeries")
+Base.read!(::AbstractRaster, ::AbstractRasterSeries) =
+    error("Cant read Raster to RasterSeries")
 
 # Filename methods
 function Base.read!(filename::AbstractString, dst::AbstractRaster)

src/series.jl

@@ -111,10 +111,21 @@ When loading a series from a single `String` path:
 Others:
 - `refdims`: existing reference dimension/s, normally not required.
 """
-struct RasterSeries{T,N,D,R,A<:AbstractArray{T,N}} <: AbstractRasterSeries{T,N,D,A}
+struct RasterSeries{T<:Union{AbstractRaster,AbstractRasterStack},N,D<:Tuple,R<:Tuple,A<:AbstractArray{T,N}} <: AbstractRasterSeries{T,N,D,A}
     data::A
     dims::D
     refdims::R
+    # Needed for ambiguity with DD
+    function RasterSeries{T,N,D,R,A}(
+        data::AbstractArray, dims::Tuple, refdims::Tuple
+    ) where {T,N,D,R,A}
+        new{T,N,D,R,A}(data, dims, refdims)
+    end
+end
+function RasterSeries(
+    data::A, dims::D, refdims::R
+) where {A<:AbstractArray{T,N},D<:Tuple,R<:Tuple} where {T,N}
+    RasterSeries{T,N,D,R,A}(data, dims, refdims)
 end
 function RasterSeries(data::AbstractArray{<:Union{AbstractRasterStack,AbstractRaster}}, dims;
     refdims=()
@@ -208,11 +219,21 @@ end
 )
     RasterSeries(data, dims, refdims)
 end
+@inline function DD.rebuild(
+    A::RasterSeries, 
+    data::AbstractArray{<:Union{AbstractRaster,AbstractRasterStack}}, 
+    dims::Tuple,
+    refdims=(),
+    name=nothing,
+    metadata=nothing,
+)
+    RasterSeries(data, dims, refdims)
+end
 @inline function DD.rebuild(
     A::RasterSeries;
     data=parent(A), dims=dims(A), refdims=refdims(A), name=nothing, metadata=nothing,
 )
-    RasterSeries(data, dims, refdims)
+    rebuild(A, data, dims, refdims)
 end
 
 function Base.map(f, series::RasterSeries)

src/stack.jl

@@ -31,8 +31,10 @@ abstract type AbstractRasterStack{K,T,N,L} <: AbstractDimStack{K,T,N,L} end
 missingval(stack::AbstractRasterStack) = getfield(stack, :missingval)
 missingval(s::AbstractRasterStack, name::Symbol) = _singlemissingval(missingval(s), name)
 
-filename(stack::AbstractRasterStack{<:Any,<:Any,<:Any,<:NamedTuple}) = map(s -> filename(s), stack)
-filename(stack::AbstractRasterStack{<:Any,<:Any,<:Any,<:Union{FileStack,OpenStack}}) = filename(parent(stack))
+filename(stack::AbstractRasterStack{<:Any,<:Any,<:Any,<:NamedTuple}) = 
+    map(s -> filename(s), layers(stack))
+filename(stack::AbstractRasterStack{<:Any,<:Any,<:Any,<:Union{FileStack,OpenStack}}) = 
+    filename(parent(stack))
 
 isdisk(st::AbstractRasterStack) = any(isdisk, layers(st))
 
@@ -83,12 +85,15 @@ function DD.rebuild(s::AbstractRasterStack;
 end
 
 function DD.rebuild_from_arrays(
-    s::AbstractRasterStack{<:Union{FileStack{<:Any,Keys},OpenStack{<:Any,Keys}}}, das::Tuple{Vararg{AbstractDimArray}}; kw...
+    s::AbstractRasterStack{<:Union{FileStack{<:Any,Keys},OpenStack{<:Any,Keys}}}, 
+    das::Tuple{Vararg{AbstractDimArray}}; 
+    kw...
 ) where Keys
     DD.rebuild_from_arrays(s, NamedTuple{Keys}(das); kw...)
 end
 function DD.rebuild_from_arrays(
-    s::AbstractRasterStack, das::NamedTuple{<:Any,<:Tuple{Vararg{AbstractDimArray}}};
+    s::AbstractRasterStack, 
+    das::NamedTuple{<:Any,<:Tuple{Vararg{AbstractDimArray}}};
     refdims=refdims(s),
     metadata=DD.metadata(s),
     dims=nothing,
@@ -112,7 +117,7 @@ end
 
 DD.name(s::AbstractRasterStack) = keys(s)
 Base.names(s::AbstractRasterStack) = keys(s)
-Base.copy(stack::AbstractRasterStack) = map(copy, stack)
+Base.copy(stack::AbstractRasterStack) = maplayers(copy, stack)
 
 #### Stack getindex ####
 # Different to DimensionalData as we construct a Raster
@@ -205,17 +210,17 @@ function RasterStack(
     data::Union{FileStack,OpenStack,NamedTuple};
     dims::Tuple,
     refdims::Tuple=(),
-    layerdims,
+    layerdims::NamedTuple,
     metadata=nokw,
     layermetadata=nokw,
     missingval=nokw,
     kw...
 )
-    # Handle values that musbe be `NamedTuple`
+    # Handle values that must be be `NamedTuple`
     layermetadata = if layermetadata isa NamedTuple
         layermetadata
-    elseif layermetadata isa Union{Nothing,NoMetadata}
-        map(_ -> NoMetadata(), layers)
+    elseif layermetadata isa Union{NoKW,Nothing,NoMetadata}
+        map(_ -> NoMetadata(), layerdims)
     else
         throw(ArgumentError("$layermetadata is not a valid input for `layermetadata`. Try a `NamedTuple` of `Dict`, `MetaData` or `NoMetadata`"))
     end
@@ -392,10 +397,10 @@ function RasterStack(filename::AbstractString;
             l_st = _layer_stack(filename; source, name, lazy, group, replace_missing, kw...)
 
             # Maybe split the stack into separate arrays to remove extra dims.
-            if !isnokw(name)
-                map(identity, l_st)
-            else
+            if isnokw(name)
                 l_st
+            else
+                maplayers(identity, l_st)
             end
         else
             # With bands actings as layers
@@ -415,12 +420,14 @@ function RasterStack(filename::AbstractString;
     end
 end
 
+# TODO test this properly
 function DD.modify(f, s::AbstractRasterStack{<:FileStack{<:Any,K}}) where K
-    open(s) do o
+    data = open(s) do o
         map(K) do k
-            Array(parent(ost)[k])
+            f(parent(ost)[k])
         end
     end
+    rebuild(s; data)
 end
 
 # Open a single file stack

test/methods.jl

@@ -59,7 +59,7 @@ gaMi = replace_missing(ga)
     @test all(isequal.(dNaN, [NaN32 7.0f0; 2.0f0 NaN32]))
     rm("test.tif")
     stNaN = replace_missing(st, NaN32; filename="teststack.tif")
-    @test all(map(stNaN[Band(1)], (a=[NaN32 7.0f0; 2.0f0 NaN32], b=[1.0 0.4; 2.0 NaN])) do x, y
+    @test all(maplayers(stNaN[Band(1)], (a=[NaN32 7.0f0; 2.0f0 NaN32], b=[1.0 0.4; 2.0 NaN])) do x, y
         all(x .=== y)
     end)
     rm("teststack_a.tif")
@@ -755,4 +755,4 @@ test = rebuild(ga; name = :test)
 
     @test_throws "strictly positive" Rasters.sample(StableRNG(123), test, 3, skipmissing = true, replace = false)
     @test_throws "Cannot draw" Rasters.sample(StableRNG(123), test, 5, replace = false)
-end
+end