merge changes

Project.toml

@@ -1,7 +1,7 @@
 name = "Rasters"
 uuid = "a3a2b9e3-a471-40c9-b274-f788e487c689"
 authors = ["Rafael Schouten <[email protected]>"]
-version = "0.11.4"
+version = "0.11.8"
 
 [deps]
 Adapt = "79e6a3ab-5dfb-504d-930d-738a2a938a0e"
@@ -82,6 +82,7 @@ ArchGDAL = "c9ce4bd3-c3d5-55b8-8973-c0e20141b8c3"
 CFTime = "179af706-886a-5703-950a-314cd64e0468"
 CoordinateTransformations = "150eb455-5306-5404-9cee-2592286d6298"
 DataFrames = "a93c6f00-e57d-5684-b7b6-d8193f3e46c0"
+GeoDataFrames = "62cb38b5-d8d2-4862-a48e-6a340996859f"
 GRIBDatasets = "82be9cdb-ee19-4151-bdb3-b400788d9abc"
 GeometryBasics = "5c1252a2-5f33-56bf-86c9-59e7332b4326"
 Makie = "ee78f7c6-11fb-53f2-987a-cfe4a2b5a57a"
@@ -94,4 +95,4 @@ Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 
 [targets]
-test = ["Aqua", "ArchGDAL", "CFTime", "CoordinateTransformations", "DataFrames", "GeometryBasics", "GRIBDatasets", "NCDatasets", "Plots", "RasterDataSources", "SafeTestsets", "Shapefile", "Statistics", "Test", "ZarrDatasets"]
+test = ["Aqua", "ArchGDAL", "CFTime", "CoordinateTransformations", "DataFrames", "GeoDataFrames", "GeometryBasics", "GRIBDatasets", "NCDatasets", "Plots", "RasterDataSources", "SafeTestsets", "Shapefile", "Statistics", "Test", "ZarrDatasets"]

docs/src/index.md

@@ -36,17 +36,11 @@ features:
     Rasters provides a standardised interface that allows many source data types to
     be used with identical syntax.
 
-    - Scripts and packages building on Rasters.jl can treat `Raster`,
-    `RasterStack`, and `RasterSeries` as black boxes.
-    - The data could hold GeoTiff or NetCDF files, `Array`s in memory or
-        `CuArray`s on the GPU - they will all behave in the same way.
-    - `RasterStack` can be backed by a Netcdf or HDF5 file, or a `NamedTuple` of
-        `Raster` holding `.tif` files, or all `Raster` in memory.
+    - Scripts and packages building on Rasters.jl can treat `Raster`, `RasterStack`, and `RasterSeries` as black boxes.
+    - The data could hold GeoTiff or NetCDF files, `Array`s in memory or `CuArray`s on the GPU - they will all behave in the same way.
+    - `RasterStack` can be backed by a Netcdf or HDF5 file, or a `NamedTuple` of `Raster` holding `.tif` files, or all `Raster` in memory.
     - Users do not have to deal with the specifics of spatial file types.
-    - `Projected` lookups with Cylindrical projections can by indexed using other Cylindrical projections
-    by setting the `mappedcrs` keyword on construction. You don't need to know the underlying
-    projection, the conversion is handled automatically. This means lat/lon
-    `EPSG(4326)` can be used seamlessly if you need that.
+    - `Projected` lookups with Cylindrical projections can by indexed using other Cylindrical projections by setting the `mappedcrs` keyword on construction. You don't need to know the underlying projection, the conversion is handled automatically. This means lat/lon `EPSG(4326)` can be used seamlessly if you need that.
 
 ## Installation
 
@@ -104,4 +98,4 @@ To make an issue we can fix quickly (or at all) there are three key steps:
 
 Good issues are really appreciated, but they do take just a little extra effort with Rasters.jl because of this need for files.
 
-:::
+:::

ext/RastersArchGDALExt/cellsize.jl

@@ -65,7 +65,7 @@ function cellsize(dims::Tuple{<:XDim, <:YDim})
             ]))
             for xb in xbnds, yb in ybnds]
     else 
-        areas = ArchGDAL.crs2transform(crs(dims), EPSG(4326)) do transform
+        areas = ArchGDAL.crs2transform(crs(dims), EPSG(4326), order = :trad) do transform
             [_area_from_coords(
                 transform,         
                 GI.LinearRing([

ext/RastersArchGDALExt/resample.jl

@@ -55,7 +55,16 @@ function resample(A::RasterStackOrArray;
             # get crs from `to` or `x` if none was passed in
             isnothing(Rasters.crs(to)) ? Rasters.crs(A) : Rasters.crs(to)
         end
-    else
+    else # issomething(crs)
+        if crs isa String
+            error("""
+                Strings as CRS aren't yet supported.  
+                Please pass a `GeoFormatTypes.jl` CRS format, like `ESRIWellKnownText`, `ProjString`, or similar.
+
+                You can find out more about GeoFormatTypes.jl at https://juliageo.org/GeoFormatTypes.jl/stable/.
+                """
+            )
+        end
         crs
     end
     if !isnothing(crs)

ext/RastersMakieExt/plotrecipes.jl

@@ -294,9 +294,17 @@ end
 function Makie.convert_arguments(t::Makie.PointBased, A::AbstractRaster{<:Number,2})
     return Makie.convert_arguments(t, _prepare_dimarray(A))
 end
-function Makie.convert_arguments(t::SurfaceLikeCompat, A::AbstractRaster{<:Any,2})
-    return Makie.convert_arguments(t, _prepare_dimarray(A))
+@static if isdefined(Makie, :SurfaceLike)
+
+    function Makie.convert_arguments(t::SurfaceLike, A::AbstractRaster{var"#s115", 2, D}) where {var"#s115", D<:Tuple}
+        return Makie.convert_arguments(t, _prepare_dimarray(A))
+    end
+else # surfacelike is not a thing
+    Makie.convert_arguments(t::Makie.VertexGrid, A::AbstractRaster{<: Any, 2}) = Makie.convert_arguments(t, _prepare_dimarray(A))
+    Makie.convert_arguments(t::Makie.CellGrid, A::AbstractRaster{<: Any, 2}) = Makie.convert_arguments(t, _prepare_dimarray(A))
+    Makie.convert_arguments(t::Makie.ImageLike, A::AbstractRaster{<: Any, 2}) = Makie.convert_arguments(t, _prepare_dimarray(A))
 end
+
 function Makie.convert_arguments(t::Makie.DiscreteSurface, A::AbstractRaster{<:Any,2})
     return Makie.convert_arguments(t, _prepare_dimarray(A))
 end

src/Rasters.jl

@@ -145,17 +145,4 @@ include("sources/grd.jl")
 include("sources/commondatamodel.jl")
 include("extensions.jl")
 
-# Compatibility with pre-1.9 julia
-function __init__()
-    @static if !isdefined(Base, :get_extension)
-        @require ArchGDAL = "c9ce4bd3-c3d5-55b8-8973-c0e20141b8c3" include("../ext/RastersArchGDALExt/RastersArchGDALExt.jl")
-        @require CoordinateTransformations = "150eb455-5306-5404-9cee-2592286d6298" include("../ext/RastersCoordinateTransformationsExt/RastersCoordinateTransformationsExt.jl")
-        @require HDF5 = "f67ccb44-e63f-5c2f-98bd-6dc0ccc4ba2f" include("../ext/RastersHDF5Ext/RastersHDF5Ext.jl")
-        @require Makie = "ee78f7c6-11fb-53f2-987a-cfe4a2b5a57a" include("../ext/RastersMakieExt/RastersMakieExt.jl")
-        @require NCDatasets = "85f8d34a-cbdd-5861-8df4-14fed0d494ab" include("../ext/RastersNCDatasetsExt/RastersNCDatasetsExt.jl")
-        @require GRIBDatasets = "82be9cdb-ee19-4151-bdb3-b400788d9abc" include("../ext/RastersGRIBDatasetsExt/RastersGRIBDatasetsExt.jl")
-        @require RasterDataSources = "3cb90ccd-e1b6-4867-9617-4276c8b2ca36" include("../ext/RastersRasterDataSourcesExt/RastersRasterDataSourcesExt.jl")
-    end
-end
-
 end

src/methods/crop_extend.jl

@@ -67,7 +67,6 @@ savefig("build/argentina_crop_example.png"); nothing
 
 $EXPERIMENTAL
 """
-
 function crop end
 function crop(l1::RasterStackOrArray, l2::RasterStackOrArray, ls::RasterStackOrArray...; kw...)
     crop((l1, l2, ls...); kw...)
@@ -80,10 +79,10 @@ function crop(xs; to=nothing, kw...)
         map(l -> crop(l; to, kw...), xs)
     end
 end
-crop(x::RasterStackOrArray; to, geometrycolumn=nothing, kw...) = _crop_to(x, to; geometrycolumn, kw...)
+crop(x::RasterStackOrArray; to, kw...) = _crop_to(x, to; kw...)
 
 # crop `A` to values of dims of `to`
-function _crop_to(x, to; geometrycolumn, kw...)
+function _crop_to(x, to; geometrycolumn=nothing, kw...)
     ext = _extent(to; geometrycolumn)
     if isnothing(ext)
         if isnothing(dims(to))
@@ -99,7 +98,7 @@ _crop_to(A, to::RasterStackOrArray; dims=DD.dims(to), kw...) = _crop_to(A, DD.di
 _crop_to(x, to::Dimension; kw...) = _crop_to(x, (to,); kw...)
 function _crop_to(x, to::DimTuple; touches=false, kw...)
     # We can only crop to sampled dims (e.g. not categorical dims like Band)
-    sampled = reduce(to; init=()) do acc, d
+    sampled = reduce(format(to); init=()) do acc, d
         lookup(d) isa AbstractSampled ? (acc..., d) : acc
     end
     return _crop_to(x, Extents.extent(sampled); touches)

src/methods/fractional_resample.jl

@@ -0,0 +1,33 @@
+using DimensionalData
+
+fractional_resample(x, res; kw...) = fractional_resample(x; res, kw...)
+fractional_resample(xs::RasterStackOrArray...; kw...) = fractional_resample(xs; kw...)
+function fractional_resample(ser::AbstractRasterSeries, args...; kw...)
+    map(x -> resample(x, args...; kw...), ser)
+end
+function fractional_resample(A::AbstractRaster; 
+    categories::NamedTuple, to=nothing, res=nothing, size=nothing
+)
+    ds = _extent2dims(to; size, res, crs(A))
+    shared_dims = commondims(A, ds)
+    isintervals(shared_dims) || throw(ArgumentError("Fractional resampling only supported for `Intervals`."))
+    outputdims = setdims(dims(A), ds)
+    outputdata = Array{typeof(categories)}(undef, outputdims)
+    output = rebuild(A; data=outputdata, dims=outputdims, refdims=())
+    map(DimIndices(otherdims(A, ds))) do D
+        slice = view(A, D...)
+        fractions = read_fractions(slice, categories, weights)
+    end
+end
+
+function read_fractions(slice, categories, weights, indices::CartesianIndices)
+    vals = view(slice, indices)
+    map(categories) do c
+        sum(zip(vals, weights)) do (v, w)
+            _compare(c, v) * w
+        end / sum(weights)
+    end
+end
+
+_compare(c::Interval, v) = v in c
+_compare(c, v) = v == c

src/methods/mosaic.jl

@@ -166,7 +166,7 @@ mosaic!(f::Function, x::RasterStackOrArray, regions::RasterStackOrArray...; kw..
     mosaic!(f, x, regions; kw...)
 function mosaic!(f::Function, A::AbstractRaster{T}, regions;
     missingval=Rasters.missingval(A),
-    atol=_default_atol(T)
+    atol=maybe_eps(T),
 ) where T
     isnokwornothing(missingval) && throw(ArgumentError("destination array must have a `missingval`"))
     _without_mapped_crs(A) do A1
@@ -262,9 +262,3 @@ function _mosaic(span::Explicit, lookup::AbstractSampled, lookups::LookupTuple)
     newbounds = vcat(permutedims(newlower), permutedims(newupper))
     return rebuild(lookup; data=newindex, span=Explicit(newbounds))
 end
-
-# These are pretty random default, but seem to work
-_default_atol(T::Type{<:Float32}) = 100eps(T)
-_default_atol(T::Type{<:Float64}) = 1000eps(T)
-_default_atol(T::Type{<:Integer}) = T(1)
-_default_atol(::Type) = nothing

src/methods/rasterize.jl

@@ -429,7 +429,7 @@ function rasterize(reducer::typeof(DD.Statistics.mean), data; fill, kw...)
     rebuild(sums ./ counts; name=:mean)
 end
 function rasterize(data; to=nothing, fill, threaded=false, geometrycolumn=nothing, kw...)
-    r = Rasterizer(data; fill, threaded, kw...)
+    r = Rasterizer(data; fill, threaded, geometrycolumn, kw...)
     rc = RasterCreator(to, data; geometrycolumn, kw..., eltype=r.eltype, fill, missingval=r.missingval)
     allocs = r.shape == :points ? nothing : _burning_allocs(rc.to; threaded)
     return create_rasterize_dest(rc) do dest

src/methods/slice_combine.jl

@@ -30,69 +30,43 @@ slice(ser::AbstractRasterSeries, dims) = cat(map(x -> slice(x, dims), ser)...; d
     throw(ArgumentError("Dimensions $(map(name, targets)) were not found in $(map(name, dims))"))
 
 """
-    combine(A::Union{AbstractRaster,AbstractRasterStack,AbstracRasterSeries}, [dims]) => Raster
+    combine(A::AbstracRasterSeries; [dims], [lazy]) => Raster
 
 Combine a `RasterSeries` along some dimension/s, creating a new `Raster` or `RasterStack`,
 depending on the contents of the series.
 
 If `dims` are passed, only the specified dimensions will be combined
 with a `RasterSeries` returned, unless `dims` is all the dims in the series.
 
+If `lazy`, concatenate lazily. The default is to concatenate lazily for lazy `Raster`s and eagerly otherwise.
+
 $EXPERIMENTAL
 """
-function combine(ser::AbstractRasterSeries, dims)
+function combine(ser::AbstractRasterSeries; dims=nokw, kw...)
+    isnokw(dims) ? _combine(ser; kw...) : _combine(ser, dims; kw...)
+end
+combine(ser::AbstractRasterSeries, dims; kw...) = _combine(ser, dims; kw...)
+function _combine(ser::AbstractRasterSeries, dims; lazy=isdisk(ser))
     ods = otherdims(ser, dims)
     if length(ods) > 0
-        map(DimIndices(ods)) do D
-            combine(view(ser, D...))
-        end |> RasterSeries
-    else
-        combine(ser)
-    end
-end
-function combine(ser::AbstractRasterSeries{<:Any,N}) where N
-    r1 = first(ser)
-    dest = _alloc_combine_dest(ser)
-    for sD in DimIndices(ser)
-        rD = map(d -> rebuild(d, :), DD.dims(r1))
-        source = ser[sD...]
-        if dest isa RasterStack
-            foreach(layers(source), layers(dest)) do source_r, dest_r
-                view(dest_r, rD..., sD...) .= source_r
-            end
-        else
-            # TODO we shouldn't need a view here??
-            view(dest, rD..., sD...) .= source
-        end
-    end
-    return dest
-end
-
-_alloc_combine_dest(s::AbstractRasterSeries) = _alloc_combine_dest(first(s), s)
-_alloc_combine_dest(r::AbstractRaster, s) = similar(r, (dims(r)..., dims(s)...))
-_alloc_combine_dest(r::AbstractRasterStack, s) = map(r -> _alloc_combine_dest(r, s), first(s))
-
-function _maybereshape(A::AbstractRaster{<:Any,N}, acc, dim) where N
-    if ndims(acc) != ndims(A)
-        newdata = reshape(parent(A), Val{N+1}())
-        d = if hasdim(refdims(A), dim)
-            dims(refdims(A), dim)
-        else
-            DD.basetypeof(dim)(1:1; lookup=NoLookup())
-        end
-        newdims = (DD.dims(A)..., d)
-        return rebuild(A; data=newdata, dims=newdims)
+        map(x -> combine(x; lazy), eachslice(ser; dims = ods))
     else
-        return A
+        combine(ser; lazy)
     end
 end
-function _maybereshape(st::AbstractRasterStack, acc, dim)
-    map((s, a) -> _maybereshape(s, a, dim), st, acc)
-end
-
-# See iterate(::GridChunks) in Diskarrays.jl
-function _chunk_inds(g, ichunk)
-    outinds = map(ichunk.I, g.chunksize, g.parentsize, g.offset) do ic, cs, ps, of
-        max((ic - 1) * cs + 1 -of, 1):min(ic * cs - of, ps)
-    end
+function _combine(ser::AbstractRasterSeries; lazy = isdisk(ser))
+    ras1 = first(ser)
+    alldims = (dims(ras1)..., dims(ser)...)
+    ser_res = DD._insert_length_one_dims(ser, alldims)
+    data = DA.ConcatDiskArray(ser_res)
+    data = lazy ? data : collect(data)
+    rebuild(ras1; data, dims = alldims, refdims = otherdims(dims(ras1, alldims)))
 end
+function _combine(ser::AbstractRasterSeries{<:AbstractRasterStack{K}}; kw...) where K
+    r1 = first(ser)
+    new_layers = map(K) do k
+        _combine(map(s -> s[k], ser); kw...)
+    end |> NamedTuple{K}
+    newdims = combinedims(new_layers...)
+    DimensionalData.rebuild_from_arrays(r1, new_layers; refdims=otherdims(dims(r1), newdims))
+end

src/methods/zonal.jl

@@ -19,7 +19,10 @@ These can be used when `of` is or contains (a) GeoInterface.jl compatible object
     where the line `:touches` the pixel, or that are completely `:inside` inside the polygon.
     The default is `:center`.
 - `progress`: show a progress bar, `true` by default, `false` to hide..
-
+- `skipmissing`: wether to apply `f` to the result of `skipmissing(A)` or not. If `true`
+    `f` will be passed an iterator over the values, which loses all spatial information.
+    if `false` `f` will be passes a masked `Raster` or `RasterStack`, and will be responsible
+    for handling missing values itself. The default value is `true`.
 
 # Example
 
@@ -71,15 +74,18 @@ insertcols!(january_stats, 1, :country => first.(split.(countries.ADMIN, r"[^A-Z
 """
 zonal(f, x::RasterStackOrArray; of, kw...) = _zonal(f, x, of; kw...)
 
-_zonal(f, x::RasterStackOrArray, of::RasterStackOrArray) = _zonal(f, x, Extents.extent(of))
-_zonal(f, x::RasterStackOrArray, of::DimTuple) = _zonal(f, x, Extents.extent(of))
+_zonal(f, x::RasterStackOrArray, of::RasterStackOrArray; kw...) = 
+    _zonal(f, x, Extents.extent(of); kw...)
+_zonal(f, x::RasterStackOrArray, of::DimTuple; kw...) = 
+    _zonal(f, x, Extents.extent(of); kw...)
 # We don't need to `mask` with an extent, it's square so `crop` will do enough.
-_zonal(f, x::Raster, of::Extents.Extent) = f(skipmissing(crop(x; to=of, touches=true)))
-function _zonal(f, x::RasterStack, ext::Extents.Extent)
+_zonal(f, x::Raster, of::Extents.Extent; skipmissing=true) =
+    _maybe_skipmissing_call(f, crop(x; to=of, touches=true), skipmissing)
+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
-        f(skipmissing(A))
+        _maybe_skipmissing_call(f, A, skipmissing)
     end
 end
 # Otherwise of is a geom, table or vector
@@ -89,22 +95,28 @@ _zonal(f, x, ::GI.AbstractFeatureCollectionTrait, fc; kw...) =
     _zonal(f, x, nothing, fc; kw...)
 _zonal(f, x::RasterStackOrArray, ::GI.AbstractFeatureTrait, feature; kw...) =
     _zonal(f, x, GI.geometry(feature); kw...)
-function _zonal(f, x::AbstractRaster, ::GI.AbstractGeometryTrait, geom; kw...)
+function _zonal(f, x::AbstractRaster, ::GI.AbstractGeometryTrait, geom; 
+    skipmissing=true, kw...
+)
     cropped = crop(x; to=geom, touches=true)
     prod(size(cropped)) > 0 || return missing
     masked = mask(cropped; with=geom, kw...)
-    return f(skipmissing(masked))
+    return _maybe_skipmissing_call(f, masked, skipmissing)
 end
-function _zonal(f, st::AbstractRasterStack, ::GI.AbstractGeometryTrait, geom; kw...)
+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)
     masked = mask(cropped; with=geom, kw...)
     return map(masked) do A
         prod(size(A)) > 0 || return missing
-        f(skipmissing(A))
+        _maybe_skipmissing_call(f, A, skipmissing)
     end
 end
-function _zonal(f, x::RasterStackOrArray, ::Nothing, data; progress=true, threaded=true, geometrycolumn=nothing, kw...)
+function _zonal(f, x::RasterStackOrArray, ::Nothing, data; 
+    progress=true, threaded=true, geometrycolumn=nothing, kw...
+)
     geoms = _get_geometries(data, geometrycolumn)
     n = length(geoms)
     n == 0 && return []
@@ -136,3 +148,5 @@ function _alloc_zonal(f, x, geoms, n; kw...)
     zs[n_missing + 1] = z1
     return zs, n_missing + 1
 end
+
+_maybe_skipmissing_call(f, A, sm) = sm ? f(skipmissing(A)) : f(A)

src/series.jl

@@ -24,6 +24,7 @@ abstract type AbstractRasterSeries{T,N,D,A} <: AbstractDimArray{T,N,D,A} end
 DD.metadata(A::AbstractRasterSeries) = NoMetadata()
 DD.name(A::AbstractRasterSeries) = NoName()
 DD.label(A::AbstractRasterSeries) = ""
+isdisk(A::AbstractRasterSeries) = any(isdisk, A)
 
 """
     modify(f, series::AbstractRasterSeries)

src/stack.jl

@@ -34,7 +34,7 @@ missingval(s::AbstractRasterStack, name::Symbol) = _singlemissingval(missingval(
 filename(stack::AbstractRasterStack{<:Any,<:Any,<:Any,<:NamedTuple}) = map(s -> filename(s), stack)
 filename(stack::AbstractRasterStack{<:Any,<:Any,<:Any,<:Union{FileStack,OpenStack}}) = filename(parent(stack))
 
-isdisk(st::AbstractRasterStack) = isdisk(layers(st, 1))
+isdisk(st::AbstractRasterStack) = any(isdisk, layers(st))
 
 setcrs(x::AbstractRasterStack, crs) = set(x, setcrs(dims(x), crs)...)
 setmappedcrs(x::AbstractRasterStack, mappedcrs) = set(x, setmappedcrs(dims(x), mappedcrs)...)