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.29.4"
 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/classify.jl

@@ -152,7 +152,11 @@ function classify!(A::AbstractRaster, pairs;
     end
     return rebuild(out; missingval=missingval)
 end
-function classify!(xs::RasterSeriesOrStack, pairs...; kw...)
+function classify!(xs::RasterStack, pairs...; kw...)
+    maplayers(x -> classify!(x, pairs...; kw...),  xs)
+    return xs
+end
+function classify!(xs::RasterSeries, pairs...; kw...)
     map(x -> classify!(x, pairs...; kw...),  xs)
     return xs
 end

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

@@ -74,7 +74,7 @@ function _mosaic(f::Function, ::AbstractRaster, regions;
     l1 = first(regions)
     A = create(filename, T, dims; name=name(l1), missingval, metadata=metadata(l1))
     open(A; write=true) do a
-        mosaic!(f, a, regions; missingval, kw...)
+        _mosaic!(f, a, regions; missingval, kw...)
     end
     return A
 end
@@ -135,8 +135,10 @@ 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, dest::RasterStackOrArray, regions::RasterStackOrArray...; kw...) =
+    _mosaic!(f, dest, regions; kw...)
+
+function _mosaic!(f::Function, A::AbstractRaster{T}, regions::Union{Tuple,AbstractArray};
     missingval=missingval(A), atol=maybe_eps(T)
 ) where T
     _without_mapped_crs(A) do A1
@@ -169,10 +171,11 @@ function mosaic!(f::Function, A::AbstractRaster{T}, regions;
     end
     return A
 end
-function mosaic!(f::Function, st::AbstractRasterStack, regions::Tuple; kw...)
-    map(st, regions...) do A, r...
-        mosaic!(f, A, r; kw...)
+function _mosaic!(f::Function, st::AbstractRasterStack, regions::Union{Tuple,AbstractArray}; kw...)
+    map(values(st), map(values, regions)...) do A, r...
+        mosaic!(f, A, r...; kw...)
     end
+    return st
 end
 
 _mosaic(alldims::Tuple{<:DimTuple,Vararg{DimTuple}}) = map(_mosaic, alldims...)

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)
@@ -831,6 +831,7 @@ end
 # We get 64 Bool values to a regular `Int` meaning this doesn't scale too
 # badly for large tables of geometries. 64k geometries and a 1000 * 1000
 # raster needs 1GB of memory just for the `BitArray`.
+# TODO combine these theyre nearly the same
 function _reduce_bitarray!(f, st::AbstractRasterStack, geoms, fill::NamedTuple, r::Rasterizer, allocs)
     (; lock, shape, boundary, verbose, progress, threaded) = r
     # Define mask dimensions, the same size as the spatial dims of x
@@ -840,7 +841,7 @@ function _reduce_bitarray!(f, st::AbstractRasterStack, geoms, fill::NamedTuple,
     # Use a generator over the array axis in case the iterator has no length
     geom_axis = axes(masks, Dim{:geometry}())
     fill = map(itr -> [v for (_, v) in zip(geom_axis, itr)], fill)
-    T = NamedTuple{keys(st),Tuple{map(eltype, st)...}}
+    T = eltype(st)
     range = axes(st, Y())
     _run(range, threaded, progress, "Reducing...") do y
         _reduce_bitarray_loop(f, st, T, fill, masks, y)

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=()
@@ -205,14 +216,25 @@ end
 
 @inline function DD.rebuild(
     A::RasterSeries, data, dims::Tuple, refdims=(), name=nothing, metadata=nothing,
+)
+    # if `data` is not an AbstractArray of Raster or RasterStacks, return a Raster
+    Raster(data, dims; refdims, name, metadata)
+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)