add continue_patch changes more cleanly: fixes Makie for new DD (#780)

* add continue_patch commits here

* fix refdims values in testts

* fix Vararg warnings

* bugfix concat

* fix concat with parent

* min bound DD at 0.28.2

* compare parents in reproject

* fix comparison

* bugfix comparisons

* bugfix tests

* fix allocs?

* fix allocs

* run code once before testing allocations

* allow more allocations for 1.11

Project.toml

@@ -52,15 +52,15 @@ CommonDataModel = "0.2.3, 0.3"
 ConstructionBase = "1"
 CoordinateTransformations = "0.6.2"
 DataFrames = "1"
-DimensionalData = "0.27.3"
+DimensionalData = "0.28.2"
 DiskArrays = "0.3, 0.4"
 Extents = "0.1"
 FillArrays = "0.12, 0.13, 1"
 Flatten = "0.4"
 GRIBDatasets = "0.2, 0.3"
 GeoFormatTypes = "0.4"
 GeoInterface = "1"
-Makie = "0.19, 0.20, 0.21"
+Makie = "0.20, 0.21"
 Missings = "0.4, 1"
 NCDatasets = "0.13, 0.14"
 OffsetArrays = "1"

docs/src/.vitepress/config.mts

@@ -29,7 +29,7 @@ export default defineConfig({
   themeConfig: {
     outline: 'deep',
     // https://vitepress.dev/reference/default-theme-config
-    logo: 'REPLACE_ME_DOCUMENTER_VITEPRESS',
+    logo: '/logo.png',
     search: {
       provider: 'local',
       options: {

docs/src/.vitepress/theme/style.css

@@ -279,4 +279,17 @@ kbd {
   .VPDoc {
     padding-left: 25px !important;  
   }
+}
+/* Component: Docstring Custom Block */
+
+.jldocstring.custom-block {
+  border: 1px solid var(--vp-c-gray-2);
+  color: var(--vp-c-text-1)
+}
+
+.jldocstring.custom-block summary {
+  font-weight: 700;
+  cursor: pointer;
+  user-select: none;
+  margin: 0 0 8px;
 }

docs/src/plot_makie.md

@@ -17,14 +17,10 @@ fig, ax, _ = plot(A)
 contour(fig[1, 2], A)
 ax = Axis(fig[2, 1]; aspect = DataAspect())
 contourf!(ax, A)
-# surface(fig[2, 2], A) # even a 3D plot works! # broken mutating method
+surface(fig[2, 2], A) # even a 3D plot work!
 fig
 ````
 
-even a 3D plot works!
-````@example makie
-surface(A)
-````
 ## 3-D rasters in Makie
 
 !!! warning
@@ -59,10 +55,9 @@ Makie.set_theme!()
 `Rasters.rplot` should support Observable input out of the box, but the dimensions of that input
 must remain the same - i.e., the element names of a RasterStack must remain the same.
 
-````julia
+````@example makie
 Makie.set_theme!(Rasters.theme_rasters())
 # `stack` is the WorldClim climate data for January
-# observables not working here, same error as for contourf: ERROR: MethodError: no method matching typemax(::Type{ColorTypes.RGBA{Float32}})
 stack_obs = Observable(stack)
 fig = Rasters.rplot(stack_obs;
     Colorbar=(; height=Relative(0.75), width=5)
@@ -72,7 +67,9 @@ record(fig, "rplot.mp4", 1:12; framerate = 3) do i
 end 
 ````
 
-<!-- <video src="./rplot.mp4" controls="controls" autoplay="autoplay"></video> -->
+```@raw html
+<video src="./rplot.mp4" controls="controls" autoplay="autoplay"></video>
+```
 
 ````@example makie
 Makie.set_theme!() # reset theme

docs/src/plotting.md

@@ -129,7 +129,7 @@ Then load raster data. We load some worldclim layers using `RasterDataSources` v
 
 ````@example plots
 using Rasters, RasterDataSources
-
+using Dates
 climate = RasterStack(WorldClim{Climate}, (:tmin, :tmax, :prec, :wind); month=July)
 ````
 

ext/RastersMakieExt/plotrecipes.jl

@@ -1,8 +1,6 @@
 
 const HIDE_DEC = (; label=true, grid=false, minorgrid=false, minorticks=false)
 
-const SurfaceLikeCompat = isdefined(Makie, :SurfaceLike) ? Makie.SurfaceLike : Union{Makie.VertexGrid,Makie.CellGrid,Makie.ImageLike}
-
 function Rasters.style_rasters()
     Makie.Attributes(
         Axis=(
@@ -288,22 +286,21 @@ function Makie.plottype(raster::AbstractRaster{<:Union{Missing,Real},3})
     end
 end
 
+# ## `convert_arguments`
+# We need to handle missing values properly here, which DimensionalData.jl can't.  
+# That's why we have the `_prepare_dimarray` function, and also why this extension is necessary.
+
+# 1d
 function Makie.convert_arguments(t::Makie.PointBased, A::AbstractRaster{<:Any,1})
     return Makie.convert_arguments(t, _prepare_dimarray(A))
 end
 function Makie.convert_arguments(t::Makie.PointBased, A::AbstractRaster{<:Number,2})
     return Makie.convert_arguments(t, _prepare_dimarray(A))
 end
-@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
+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))
 
 function Makie.convert_arguments(t::Makie.VolumeLike, A::AbstractRaster{<:Any,3}) 
     return Makie.convert_arguments(t, _prepare_dimarray(A))
@@ -319,10 +316,13 @@ function Makie.convert_arguments(x::Makie.ConversionTrait, raster::AbstractRaste
         return Makie.convert_arguments(x, view(raster, rebuild(D, 1)))
     end
 end
+
 function Makie.convert_arguments(x::Makie.ConversionTrait, series::AbstractRasterSeries)
     return Makie.convert_arguments(x, first(series))
 end
 
+# ## Utility / helper functions
+
 function _series_dim(A)
     spatialdims = (X(), Y(), Z())
     last((dims(A, spatialdims)..., otherdims(A, spatialdims)...))

src/array.jl

@@ -105,7 +105,7 @@ function DD.modify(f, A::AbstractRaster)
     return rebuild(A, newdata)
 end
 
-function DD.DimTable(As::Tuple{<:AbstractRaster,Vararg{<:AbstractRaster}}...)
+function DD.DimTable(As::Tuple{<:AbstractRaster,Vararg{AbstractRaster}}...)
     DD.DimTable(DimStack(map(read, As...)))
 end
 
@@ -372,7 +372,7 @@ end
 filekey(ds, name) = name
 filekey(filename::String) = Symbol(splitext(basename(filename))[1])
 
-DD.dimconstructor(::Tuple{<:Dimension{<:AbstractProjected},Vararg{<:Dimension}}) = Raster
+DD.dimconstructor(::Tuple{<:Dimension{<:AbstractProjected},Vararg{Dimension}}) = Raster
 
 
 function _drop_single_band(raster, lazy::Bool)

src/methods/aggregate.jl

@@ -1,6 +1,6 @@
 
-const DimOrDimTuple = Union{Dimension,Tuple{Vararg{<:Dimension}}}
-const IntOrIntTuple = Union{Int,Tuple{Vararg{<:Int}}}
+const DimOrDimTuple = Union{Dimension,Tuple{Vararg{Dimension}}}
+const IntOrIntTuple = Union{Int,Tuple{Vararg{Int}}}
 
 struct Ag end
 struct DisAg end
@@ -38,6 +38,7 @@ $PROGRESS_KEYWORD
 
 ```jldoctest
 using Rasters, RasterDataSources, Statistics, Plots
+import ArchGDAL
 using Rasters: Center
 st = read(RasterStack(WorldClim{Climate}; month=1))
 ag = aggregate(Center(), st, (Y(20), X(20)); skipmissingval=true, progress=false)

src/methods/burning/line.jl

@@ -178,7 +178,7 @@ function _burn_line!(A::AbstractRaster, line, fill)
     end
     return n_on_line
 end
-function _burn_line!(A::AbstractRaster, line, fill, order::Tuple{Vararg{<:Dimension}})
+function _burn_line!(A::AbstractRaster, line, fill, order::Tuple{Vararg{Dimension}})
     msg = """"
         Converting a `:line` geometry to raster is currently only implemented for 2d lines.
         Make a Rasters.jl github issue if you need this for more dimensions.

src/methods/mosaic.jl

@@ -175,7 +175,7 @@ function mosaic!(f::Function, st::AbstractRasterStack, regions::Tuple; kw...)
     end
 end
 
-_mosaic(alldims::Tuple{<:DimTuple,Vararg{<:DimTuple}}) = map(_mosaic, alldims...)
+_mosaic(alldims::Tuple{<:DimTuple,Vararg{DimTuple}}) = map(_mosaic, alldims...)
 function _mosaic(dims::Dimension...)
     map(dims) do d
         DD.comparedims(first(dims), d; val=false, length=false, valtype=true)

src/methods/slice_combine.jl

@@ -58,9 +58,9 @@ 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 = DA.ConcatDiskArray(map(parent, parent(ser_res)))
     data = lazy ? data : collect(data)
-    rebuild(ras1; data, dims = alldims, refdims = otherdims(dims(ras1, alldims)))
+    rebuild(ras1; data, dims=alldims, refdims=otherdims(dims(ras1, alldims)))
 end
 function _combine(ser::AbstractRasterSeries{<:AbstractRasterStack{K}}; kw...) where K
     r1 = first(ser)
@@ -69,4 +69,4 @@ function _combine(ser::AbstractRasterSeries{<:AbstractRasterStack{K}}; kw...) wh
     end |> NamedTuple{K}
     newdims = combinedims(new_layers...)
     DimensionalData.rebuild_from_arrays(r1, new_layers; refdims=otherdims(dims(r1), newdims))
-end
+end

src/stack.jl

@@ -83,7 +83,7 @@ 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
@@ -227,15 +227,15 @@ function RasterStack(
     return _postprocess_stack(st, dims; kw...)
 end
 # Convert Tuple/Array of array to NamedTuples using name/key
-function RasterStack(data::Tuple{Vararg{<:AbstractArray}}, dims::Tuple;
+function RasterStack(data::Tuple{Vararg{AbstractArray}}, dims::Tuple;
     name::Union{Tuple,AbstractArray,NamedTuple,Nothing}=nokw, 
     kw...
 )
     isnokw(name) && throw(ArgumentError("Pass a Tuple, Array or NamedTuple of names to the `name` keyword"))
     return RasterStack(NamedTuple{cleankeys(name)}(data), dims; kw...)
 end
 # Multi Raster stack from NamedTuple of AbstractArray
-function RasterStack(data::NamedTuple{<:Any,<:Tuple{Vararg{<:AbstractArray}}}, dims::Tuple; 
+function RasterStack(data::NamedTuple{<:Any,<:Tuple{Vararg{AbstractArray}}}, dims::Tuple; 
     layerdims=nokw,
     kw...
 )
@@ -254,15 +254,15 @@ end
 # Multi Raster stack from AbstractDimArray splat
 RasterStack(layers::AbstractDimArray...; kw...) = RasterStack(layers; kw...)
 # Multi Raster stack from tuple with `name` keyword
-function RasterStack(layers::Tuple{Vararg{<:AbstractDimArray}};
+function RasterStack(layers::Tuple{Vararg{AbstractDimArray}};
     name=map(name, layers),
     kw...
 )
     RasterStack(NamedTuple{cleankeys(name)}(layers); kw...)
 end
 # Multi RasterStack from NamedTuple
 # This method is called after most other RasterStack methods.
-function RasterStack(layers::NamedTuple{K,<:Tuple{Vararg{<:AbstractDimArray}}};
+function RasterStack(layers::NamedTuple{K,<:Tuple{Vararg{AbstractDimArray}}};
     resize::Union{Function,NoKW}=nokw,
     _layers=resize isa NoKW ? layers : resize(layers),
     dims::Tuple=DD.combinedims(_layers...),

src/table_ops.jl

@@ -8,7 +8,7 @@ function _auto_dim_columns(table, dims::Tuple)
 end
 
 _not_a_dimcol(table, dimcols::DimTuple) = _not_a_dimcol(table, map(name, dimcols))
-_not_a_dimcol(table, dimcols::Tuple{Vararg{<:Pair}}) = _not_a_dimcol(table, map(last, dimcols))
+_not_a_dimcol(table, dimcols::Tuple{Vararg{Pair}}) = _not_a_dimcol(table, map(last, dimcols))
 _not_a_dimcol(table, dimcols::Tuple{}) = ()
 function _not_a_dimcol(table, dimcols::Tuple{Vararg{Symbol}})
     dimcols = (dimcols..., :Band) # TODO deal with how annoying `Band` is 

test/reproject.jl

@@ -36,8 +36,8 @@ using Rasters: reproject, convertlookup
     reprojected_raster = reproject(raster; crs=projcea)
     # Dims have changed
     @test dims(reprojected_raster) == (x1, y1)
-    # Raster has not changed
-    @test reprojected_raster == raster
+    # Parent data has not changed
+    @test parent(reprojected_raster) == parent(raster)
 
     @test_throws ArgumentError reproject(cea, EPSG(32618), Y(), [-3.658789324855012e6])
     @test_throws ArgumentError reproject(cea, EPSG(32618), X(), [-3.658789324855012e6])

test/series.jl

@@ -76,6 +76,7 @@ end
     r1 = Raster(ones(4, 5, 10), (X(), Y(), Ti(10:10:100))) .* reshape(1:10, (1, 1, 10))
     r2 = r1 .* 2
     ser = slice(r1, Ti)
+
     @test size(ser) == (10,)
     combined = Rasters.combine(ser, Ti())
     @test combined == r1

test/sources/gdal.jl

@@ -422,7 +422,7 @@ gdalpath = maybedownload(url)
             @test (@allocations write(filename2, gdalarray; force=true)) < 1e4
             saved = Raster(filename2; crs=crs(gdalarray), mappedcrs=crs(gdalarray))
             @test size(saved) == size(gdalarray)
-            @test saved ≈ gdalarray
+            @test parent(saved) ≈ parent(gdalarray)
             clat, clon = DimensionalData.shiftlocus.(Ref(Center()), dims(gdalarray, (Y, X)))
             @test index(clat) ≈ index(saved, Y)
             @test index(clon) ≈ index(saved, X)
@@ -700,7 +700,7 @@ end
             # Test forcing
             write(filename, gdalstack; force=true);
             saved = RasterStack(filename);
-            @test all(read(saved[:a]) .== geoA)
+            @test all(parent(read(saved[:a])) .== parent(geoA))
             rm(filename)
         end
 

test/sources/grd.jl

@@ -224,7 +224,8 @@ grdpath = stem * ".gri"
             # 1 band is added again on save
             @test size(saved) == size(grdarray[Band(1)])
             @test parent(saved) == parent(grdarray[Band(1)])
-            @test (@allocations write(filename2, grdarray[Band(1)]; force = true)) < 1e3
+            write(filename2, grdarray; force = true)
+            @test (@allocations write(filename2, grdarray; force = true)) < 3e3
         end
 
         @testset "3d with subset" begin
@@ -250,7 +251,8 @@ grdpath = stem * ".gri"
             @test all(parent(saved) .=== parent(geoA))
             @test saved isa typeof(geoA)
             @test parent(saved) == parent(geoA)
-            @test (@allocations write(filename, GRDsource(), geoA; force = true)) < 1e3
+            write(filename, GRDsource(), geoA; force = true)
+            @test (@allocations write(filename, GRDsource(), geoA; force = true)) < 3e3
         end
 
         @testset "to netcdf" begin
@@ -259,20 +261,21 @@ grdpath = stem * ".gri"
             write(filename2, grdarray[Band(1)]; force = true)
             saved = Raster(filename2; crs=crs(grdarray))
             @test size(saved) == size(grdarray[Band(1)])
-            @test all(replace_missing(saved, missingval(grdarray)) .≈ grdarray[Band(1)])
+            @test all(parent(replace_missing(saved, missingval(grdarray))) .≈ parent(grdarray[Band(1)]))
             @test index(saved, X) ≈ index(grdarray, X) .+ 0.5
             @test index(saved, Y) ≈ index(grdarray, Y) .+ 0.5
             @test bounds(saved, Y) == bounds(grdarray, Y)
             @test bounds(saved, X) == bounds(grdarray, X)
-            @test (@allocations write(filename2, grdarray[Band(1)]; force = true)) < 1e3
+            write(filename2, grdarray; force = true)
+            @test (@allocations write(filename2, grdarray; force = true)) < 3e3
 
         end
 
         @testset "to gdal" begin
             # No Band
             gdalfilename = tempname() * ".tif"
             write(gdalfilename, GDALsource(), grdarray[Band(1)]; force = true)
-            @test (@allocations write(gdalfilename, GDALsource(), grdarray[Band(1)]; force = true)) < 1e4
+            @test (@allocations write(gdalfilename, GDALsource(), grdarray[Band(1)]; force = true)) < 3e4
             gdalarray = Raster(gdalfilename)
             # @test convert(ProjString, crs(gdalarray)) == convert(ProjString, EPSG(4326))
             @test val(dims(gdalarray, X)) ≈ val(dims(grdarray, X))

test/sources/ncdatasets.jl

@@ -331,8 +331,9 @@ end
             # Tiff locus = Start, Netcdf locus = Center
             @test index(gdalarray, Y) .+ 0.5 ≈ index(nccleaned, Y)
             @test index(gdalarray, X) .+ 1.0  ≈ index(nccleaned, X)
-            @test gdalarray ≈ nccleaned
+            @test parent(gdalarray) ≈ parent(nccleaned)
         end
+
         @testset "to grd" begin
             nccleaned = replace_missing(ncarray[Ti(1)], -9999.0)
             write("testgrd.gri", nccleaned; force=true)
@@ -342,9 +343,9 @@ end
             @test bounds(grdarray) == bounds(nccleaned)
             @test index(grdarray, Y) ≈ reverse(index(nccleaned, Y)) .- 0.5
             @test index(grdarray, X) ≈ index(nccleaned, X) .- 1.0
-            @test reverse(grdarray; dims=Y) ≈ nccleaned
-            # rm("testgrd.gri")
-            # rm("testgrd.grd")
+            @test parent(reverse(grdarray; dims=Y)) ≈ parent(nccleaned)
+            rm("testgrd.gri")
+            rm("testgrd.grd")
         end
 
         @testset "write points" begin