From c67ac1b15eb0563f1d85070356f449492f155560 Mon Sep 17 00:00:00 2001
From: Rafael Schouten <rafaelschouten@gmail.com>
Date: Sat, 12 Oct 2024 08:50:28 +0200
Subject: [PATCH 1/9] allow single dim constructor like DimArray (#798)

---
 src/array.jl  | 1 +
 test/array.jl | 5 +++++
 2 files changed, 6 insertions(+)

diff --git a/src/array.jl b/src/array.jl
index c40c57bf..7b0a02d7 100644
--- a/src/array.jl
+++ b/src/array.jl
@@ -269,6 +269,7 @@ function Raster(A::AbstractArray{T,1}, dims::Tuple{<:Dimension,<:Dimension,Varar
 )::Raster{T,length(dims)} where T
     Raster(reshape(A, map(length, dims)), dims; kw...)
 end
+Raster(A::AbstractArray{<:Any,1}, dim::Dimension; kw...) = Raster(A, (dim,); kw...)
 function Raster(table, dims::Tuple;
     name=nokw,
     kw...
diff --git a/test/array.jl b/test/array.jl
index f9c3c9d0..e7716709 100644
--- a/test/array.jl
+++ b/test/array.jl
@@ -22,6 +22,11 @@ ga1 = Raster(data1; dims=dims1, refdims=refdimz, name=nme, metadata=meta, missin
     @test_throws ArgumentError Raster("notafile", dims1)
 end
 
+@testset "Single dim vector constructor" begin
+    A = [1, 2, 3]
+    @test Raster(A, X()) === Raster(A, (X(),))
+end
+
 @testset "array properties" begin
     @test name(ga1) == :test
     @test missingval(ga1) == -9999.0

From b677b0551ce9d52c9eea817de7e330e936c32cd4 Mon Sep 17 00:00:00 2001
From: Anshul Singhvi <anshulsinghvi@gmail.com>
Date: Fri, 11 Oct 2024 23:52:08 -0700
Subject: [PATCH 2/9] forgot a space in showerror(::BackendException) (#795)

---
 src/sources/sources.jl | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/src/sources/sources.jl b/src/sources/sources.jl
index 1c2d0dc4..0c90cef6 100644
--- a/src/sources/sources.jl
+++ b/src/sources/sources.jl
@@ -61,7 +61,7 @@ function Base.showerror(io::IO, e::BackendException)
     printstyled(io, "Rasters.jl"; underline = true)
     printstyled(io, " requires backends to be loaded manually.  Run ")
     printstyled(io, "`import $(e.backend)`"; bold = true)
-    print(io, "to fix this error.")
+    print(io, " to fix this error.")
 end
 
 # Get the source backend for a file extension, falling back to GDALsource

From 8cf4942c4a8f80a95cb34821e6cc6d95373dbcb1 Mon Sep 17 00:00:00 2001
From: Tiem van der Deure <tiemvanderdeure@gmail.com>
Date: Sat, 12 Oct 2024 12:44:48 +0200
Subject: [PATCH 3/9] do not unnecessarily return a union in extract (#790)

* extract tweaks

* just write kw... in docs

Co-authored-by: Rafael Schouten <rafaelschouten@gmail.com>

* use Type{G}

Co-authored-by: Rafael Schouten <rafaelschouten@gmail.com>

* add a missing where G

* merge _rowtype changes from other branch

* optimization for extract points with skipmissing

---------

Co-authored-by: Rafael Schouten <rafaelschouten@gmail.com>
---
 src/methods/extract.jl | 106 +++++++++++------------------------------
 src/utils.jl           |  60 +++++++++++++++++++++++
 test/methods.jl        |   6 +--
 3 files changed, 92 insertions(+), 80 deletions(-)

diff --git a/src/methods/extract.jl b/src/methods/extract.jl
index 898bae16..5fd748ee 100644
--- a/src/methods/extract.jl
+++ b/src/methods/extract.jl
@@ -1,11 +1,5 @@
-using DimensionalData.Lookups: _True, _False
-
-_booltype(x) = x ? _True() : _False()
-istrue(::_True) = true
-istrue(::_False) = false
-
 """
-    extract(x, data; atol)
+    extract(x, data; kw...)
 
 Extracts the value of `Raster` or `RasterStack` at given points, returning
 an iterable of `NamedTuple` with properties for `:geometry` and raster or
@@ -78,7 +72,7 @@ function extract end
 end
 
 function _extract(A::RasterStackOrArray, geom::Missing, names, kw...)
-    T = _rowtype(A, geom; names, kw...)
+    T = _extractrowtype(A, geom; names, kw...)
     [_maybe_add_fields(T, map(_ -> missing, names), missing, missing)]
 end
 function _extract(A::RasterStackOrArray, geom; names, kw...)
@@ -89,9 +83,9 @@ function _extract(A::RasterStackOrArray, ::Nothing, data;
 )
     geoms = _get_geometries(data, geometrycolumn)
     T = if istrue(skipmissing)
-        _rowtype(A, nonmissingtype(eltype(geoms)); names, skipmissing, kw...)
+        _extractrowtype(A, nonmissingtype(eltype(geoms)); names, skipmissing, kw...)
     else
-        _rowtype(A, eltype(geoms); names, skipmissing, kw...)
+        _extractrowtype(A, eltype(geoms); names, skipmissing, kw...)
     end
     # Handle empty / all missing cases
     (length(geoms) > 0 && any(!ismissing, geoms)) || return T[]
@@ -101,13 +95,15 @@ function _extract(A::RasterStackOrArray, ::Nothing, data;
     # We need to split out points from other geoms
     # TODO this will fail with mixed point/geom vectors
     if trait1 isa GI.PointTrait
-        rows = Vector{T}(undef, length(geoms))
         if istrue(skipmissing)
+            T2 = _extractrowtype(A, eltype(geoms), Tuple{Int64, Int64}; 
+                names, skipmissing = _False(), kw...)
+            rows = Vector{T2}(undef, length(geoms))
             j = 1
             for i in eachindex(geoms)
                 g = geoms[i]
                 ismissing(g) && continue
-                e = _extract_point(T, A, g, skipmissing; names, kw...)
+                e = _extract_point(T2, A, g, skipmissing; names, kw...)
                 if !ismissing(e) 
                     rows[j] = e
                     j += 1
@@ -115,7 +111,9 @@ function _extract(A::RasterStackOrArray, ::Nothing, data;
                 nothing
             end
             deleteat!(rows, j:length(rows))
+            rows = T === T2 ? rows : T.(rows)
         else
+            rows = Vector{T}(undef, length(geoms))
             for i in eachindex(geoms)
                 g = geoms[i]
                 rows[i] = _extract_point(T, A, g, skipmissing; names, kw...)::T
@@ -143,14 +141,14 @@ end
 function _extract(A::RasterStackOrArray, ::GI.AbstractMultiPointTrait, geom; 
     skipmissing, kw...
 )
-    T = _rowtype(A, GI.getpoint(geom, 1); names, skipmissing, kw...)
+    T = _extractrowtype(A, GI.getpoint(geom, 1); names, skipmissing, kw...)
     rows = (_extract_point(T, A, p, skipmissing; kw...) for p in GI.getpoint(geom))
     return skipmissing isa _True ? collect(_skip_missing_rows(rows, _missingval_or_missing(A), names)) : collect(rows)
 end
 function _extract(A::RasterStackOrArray, ::GI.PointTrait, geom; 
     skipmissing, kw...
 )
-    T = _rowtype(A, geom; names, skipmissing, kw...)
+    T = _extractrowtype(A, geom; names, skipmissing, kw...)
     _extract_point(T, A, geom, skipmissing; kw...)
 end
 function _extract(A::RasterStackOrArray, t::GI.AbstractGeometryTrait, geom;
@@ -168,7 +166,7 @@ function _extract(A::RasterStackOrArray, t::GI.AbstractGeometryTrait, geom;
     else
         GI.x(p), GI.y(p)
     end
-    T = _rowtype(A, tuplepoint; names, skipmissing, kw...)
+    T = _extractrowtype(A, tuplepoint; names, skipmissing, kw...)
     B = boolmask(geom; to=template, kw...)
     offset = CartesianIndex(map(x -> first(x) - 1, parentindices(parent(template))))
     # Add a row for each pixel that is `true` in the mask
@@ -290,70 +288,24 @@ Base.@assume_effects :total function _maybe_add_fields(::Type{T}, props::NamedTu
     end
 end
 
-_names(A::AbstractRaster) = (Symbol(name(A)),)
-_names(A::AbstractRasterStack) = keys(A)
-
-@inline _nametypes(::Raster{T}, ::NamedTuple{Names}, skipmissing::_True) where {T,Names} = (nonmissingtype(T),)
-@inline _nametypes(::Raster{T}, ::NamedTuple{Names}, skipmissing::_False) where {T,Names} = (Union{Missing,T},)
-# This only compiles away when generated
-@generated function _nametypes(
-    ::RasterStack{<:Any,T}, ::NamedTuple{PropNames}, skipmissing::_True
-) where {T<:NamedTuple{StackNames,Types},PropNames} where {StackNames,Types}
-    nt = NamedTuple{StackNames}(map(nonmissingtype, Types.parameters))
-    return values(nt[PropNames])
-end
-@generated function _nametypes(
-    ::RasterStack{<:Any,T}, ::NamedTuple{PropNames}, skipmissing::_False
-) where {T<:NamedTuple{StackNames,Types},PropNames} where {StackNames,Types}
-    nt = NamedTuple{StackNames}(map(T -> Union{Missing,T}, Types.parameters))
-    return values(nt[PropNames])
-end
-
-# _rowtype returns the complete NamedTuple type for a point row
-# This code is entirely for types stability and performance.
-_rowtype(x, g; kw...) = _rowtype(x, typeof(g); kw...)
-_rowtype(x, g::Type; geometry, index, skipmissing, names, kw...) = 
-    _rowtype(x, g, geometry, index, skipmissing, names)
-function _rowtype(x, ::Type{G}, geometry::_True, index::_True, skipmissing::_True, names::NamedTuple{Names}) where {G,Names}
-    keys = (:geometry, :index, Names...,) 
-    types = Tuple{G,Tuple{Int,Int},_nametypes(x, names, skipmissing)...}
-    NamedTuple{keys,types}
-end
-function _rowtype(x, ::Type{G}, geometry::_True, index::_False, skipmissing::_True, names::NamedTuple{Names}) where {G,Names}
-    keys = (:geometry, Names...,)
-    types = Tuple{G,_nametypes(x, names, skipmissing)...}
-    NamedTuple{keys,types}
-end
-function _rowtype(x, ::Type{G}, geometry::_False, index::_True, skipmissing::_True, names::NamedTuple{Names}) where {G,Names}
-    keys = (:index, Names...,) 
-    types = Tuple{Tuple{Int,Int},_nametypes(x, names, skipmissing)...}
-    NamedTuple{keys,types}
-end
-function _rowtype(x, ::Type{G}, geometry::_False, index::_False, skipmissing::_True, names::NamedTuple{Names}) where {G,Names}
-    keys = Names
-    types = Tuple{_nametypes(x, names, skipmissing)...}
-    NamedTuple{keys,types}
-end
-function _rowtype(x, ::Type{G}, geometry::_True, index::_True, skipmissing::_False, names::NamedTuple{Names}) where {G,Names}
-    keys = (:geometry, :index, names...,) 
-    types = Tuple{Union{Missing,G},Union{Missing,Tuple{Int,Int}},_nametypes(x, names, skipmissing)...}
-    NamedTuple{keys,types}
-end
-function _rowtype(x, ::Type{G}, geometry::_True, index::_False, skipmissing::_False, names::NamedTuple{Names}) where {G,Names}
-    keys = (:geometry, Names...,)
-    types = Tuple{Union{Missing,G},_nametypes(x, names, skipmissing)...}
-    NamedTuple{keys,types}
-end
-function _rowtype(x, ::Type{G}, geometry::_False, index::_True, skipmissing::_False, names::NamedTuple{Names}) where {G,Names}
-    keys = (:index, Names...,) 
-    types = Tuple{Union{Missing,Tuple{Int,Int}},_nametypes(x, names, skipmissing)...}
-    NamedTuple{keys,types}
+function _extractrowtype(x, g; geometry, index, skipmissing, names, kw...)
+    I = if istrue(skipmissing)
+        Tuple{Int, Int}
+    else
+        Union{Missing, Tuple{Int, Int}}
+    end
+    _extractrowtype(x, g, I; geometry, index, skipmissing, names, kw...)
 end
-function _rowtype(x, ::Type{G}, geometry::_False, index::_False, skipmissing::_False, names::NamedTuple{Names}) where {G,Names}
-    keys = Names
-    types = Tuple{_nametypes(x, names, skipmissing)...}
-    NamedTuple{keys,types}
+function _extractrowtype(x, g, ::Type{I}; geometry, index, skipmissing, names, kw...) where I
+    G = if istrue(skipmissing)
+        nonmissingtype(typeof(g))
+    else
+        typeof(g)
+    end
+    _extractrowtype(x, G, I; geometry, index, skipmissing, names, kw...)
 end
+_extractrowtype(x, ::Type{G}, ::Type{I}; geometry, index, skipmissing, names, kw...) where {G, I} =
+    _rowtype(x, G, I; geometry, index, skipmissing, names)
 
 @inline _skip_missing_rows(rows, ::Missing, names) = 
     Iterators.filter(row -> !any(ismissing, row), rows)
diff --git a/src/utils.jl b/src/utils.jl
index 2588284a..67a8d2c4 100644
--- a/src/utils.jl
+++ b/src/utils.jl
@@ -359,3 +359,63 @@ function _no_memory_error(f, bytes)
     """
     return error(msg)
 end
+
+
+# _rowtype returns the complete NamedTuple type for a point row
+# This code is entirely for types stability and performance.
+# It is used in extract and Rasters.sample
+_names(A::AbstractRaster) = (Symbol(name(A)),)
+_names(A::AbstractRasterStack) = keys(A)
+
+using DimensionalData.Lookups: _True, _False
+_booltype(x) = x ? _True() : _False()
+istrue(::_True) = true
+istrue(::_False) = false
+
+function _rowtype(x, ::Type{G}, ::Type{I}; geometry, index, skipmissing, names) where {G, I}
+    keys = _rowkeys(geometry, index, names)
+    types = _rowtypes(x, G, I, geometry, index, skipmissing, names)
+    NamedTuple{keys,types}
+end
+
+function _rowtypes(
+    x, ::Type{G}, ::Type{I}, geometry::_True, index::_True, skipmissing, names::NamedTuple{Names}
+) where {G,I,Names}
+    Tuple{G,I,_nametypes(x, names, skipmissing)...}
+end
+function _rowtypes(
+    x, ::Type{G}, ::Type{I}, geometry::_True, index::_False, skipmissing, names::NamedTuple{Names}
+) where {G,I,Names}
+    Tuple{G,_nametypes(x, names, skipmissing)...}
+end
+function _rowtypes(
+    x, ::Type{G}, ::Type{I}, geometry::_False, index::_True, skipmissing, names::NamedTuple{Names}
+) where {G,I,Names}
+    Tuple{I,_nametypes(x, names, skipmissing)...}
+end
+function _rowtypes(
+    x, ::Type{G}, ::Type{I}, geometry::_False, index::_False, skipmissing, names::NamedTuple{Names}
+) where {G,I,Names}
+    Tuple{_nametypes(x, names, skipmissing)...}
+end
+
+@inline _nametypes(::Raster{T}, ::NamedTuple{Names}, skipmissing::_True) where {T,Names} = (nonmissingtype(T),)
+@inline _nametypes(::Raster{T}, ::NamedTuple{Names}, skipmissing::_False) where {T,Names} = (Union{Missing,T},)
+# This only compiles away when generated
+@generated function _nametypes(
+    ::RasterStack{<:Any,T}, ::NamedTuple{PropNames}, skipmissing::_True
+) where {T<:NamedTuple{StackNames,Types},PropNames} where {StackNames,Types}
+    nt = NamedTuple{StackNames}(map(nonmissingtype, Types.parameters))
+    return values(nt[PropNames])
+end
+@generated function _nametypes(
+    ::RasterStack{<:Any,T}, ::NamedTuple{PropNames}, skipmissing::_False
+) where {T<:NamedTuple{StackNames,Types},PropNames} where {StackNames,Types}
+    nt = NamedTuple{StackNames}(map(T -> Union{Missing,T}, Types.parameters))
+    return values(nt[PropNames])
+end
+
+_rowkeys(geometry::_False, index::_False, names::NamedTuple{Names}) where Names = Names
+_rowkeys(geometry::_True, index::_False, names::NamedTuple{Names}) where Names = (:geometry, Names...)
+_rowkeys(geometry::_True, index::_True, names::NamedTuple{Names}) where Names = (:geometry, :index, Names...)
+_rowkeys(geometry::_False, index::_True, names::NamedTuple{Names}) where Names = (:index, Names...)
diff --git a/test/methods.jl b/test/methods.jl
index cd07130d..df1a81e5 100644
--- a/test/methods.jl
+++ b/test/methods.jl
@@ -399,7 +399,7 @@ createpoint(args...) = ArchGDAL.createpoint(args...)
             (index = (1, 2), test = 2,)
         ])
         # NamedTuple (reversed) points - tests a Table that iterates over points
-        T = @NamedTuple{geometry::Union{Missing,@NamedTuple{Y::Float64,X::Float64}},test::Union{Missing,Int64}}
+        T = @NamedTuple{geometry::Union{@NamedTuple{Y::Float64,X::Float64}},test::Union{Missing,Int64}}
         @test all(extract(rast, [(Y=0.1, X=9.0), (Y=0.2, X=10.0), (Y=0.3, X=10.0)]) .=== T[
             (geometry = (Y = 0.1, X = 9.0), test = 1)
             (geometry = (Y = 0.2, X = 10.0), test = 4)
@@ -412,7 +412,7 @@ createpoint(args...) = ArchGDAL.createpoint(args...)
         ])
         # Extract a polygon
         p = ArchGDAL.createpolygon([[[8.0, 0.0], [11.0, 0.0], [11.0, 0.4], [8.0, 0.0]]])
-        T = @NamedTuple{geometry::Union{Missing,Tuple{Float64,Float64}},test::Union{Missing,Int64}}
+        T = @NamedTuple{geometry::Union{Tuple{Float64,Float64}},test::Union{Missing,Int64}}
         @test all(extract(rast_m, p) .=== T[
             (geometry = (9.0, 0.1), test = 1)
             (geometry = (10.0, 0.1), test = 3)
@@ -447,7 +447,7 @@ createpoint(args...) = ArchGDAL.createpoint(args...)
             (index = (2, 1), test = 3)
             (index = (2, 2), test = missing)
         ])
-        T = @NamedTuple{geometry::Union{Missing,Tuple{Float64,Float64}},index::Union{Missing,Tuple{Int,Int}},test::Union{Missing,Int64}}
+        T = @NamedTuple{geometry::Union{Tuple{Float64,Float64}},index::Union{Missing,Tuple{Int,Int}},test::Union{Missing,Int64}}
         @test all(extract(rast_m, p; index=true) .=== T[
              (geometry = (9.0, 0.1), index = (1, 1), test = 1)
              (geometry = (10.0, 0.1), index = (2, 1), test = 3)

From 320270763f29f4d2ce8fcc16a12ecc2e388b9831 Mon Sep 17 00:00:00 2001
From: Anshul Singhvi <anshulsinghvi@gmail.com>
Date: Tue, 15 Oct 2024 11:31:39 -0700
Subject: [PATCH 4/9] Make extracting `grid_mapping` from CDM datasets safer
 (#803)

* Update commondatamodel.jl

* Add test

* move resample tests to be last

* Include commondatamodel and zarr tests

* use Zarr from ZarrDatasets.jl

* Don't run Zarr tests

* Add the julia cache action to speed up CI
---
 .github/workflows/ci.yml        |  4 ++++
 src/sources/commondatamodel.jl  | 11 +++++++++--
 test/runtests.jl                |  5 ++++-
 test/sources/commondatamodel.jl | 23 ++++++++++++++++++++++-
 test/sources/zarr.jl            |  3 ++-
 5 files changed, 41 insertions(+), 5 deletions(-)

diff --git a/.github/workflows/ci.yml b/.github/workflows/ci.yml
index 6c8b91be..644eca74 100644
--- a/.github/workflows/ci.yml
+++ b/.github/workflows/ci.yml
@@ -8,6 +8,9 @@ defaults:
 jobs:
   test:
     name: Julia ${{ matrix.version }} - ${{ matrix.os }} - ${{ matrix.arch }}
+    permissions:
+      actions: write
+      contents: read
     runs-on: ${{ matrix.os }}
     strategy:
       fail-fast: true
@@ -26,6 +29,7 @@ jobs:
         with:
           version: ${{ matrix.version }}
           arch: ${{ matrix.arch }}
+      - uses: julia-actions/cache@v2
       - uses: julia-actions/julia-buildpkg@latest
       - uses: julia-actions/julia-runtest@latest
         env:
diff --git a/src/sources/commondatamodel.jl b/src/sources/commondatamodel.jl
index 4dd1b71a..5e9084da 100644
--- a/src/sources/commondatamodel.jl
+++ b/src/sources/commondatamodel.jl
@@ -241,7 +241,14 @@ function _layermetadata(ds::AbstractDataset; layers)
     map(layers.attrs) do attr
         md = _metadatadict(_sourcetrait(ds), attr)
         if haskey(attr, "grid_mapping")
-            md["grid_mapping"] = Dict(CDM.attribs(ds[attr["grid_mapping"]]))
+            if haskey(ds, attr["grid_mapping"])
+                md["grid_mapping"] = Dict(CDM.attribs(ds[attr["grid_mapping"]]))
+            else
+                global_attrs = CDM.attribs(ds)
+                if haskey(global_attrs, attr["grid_mapping"])
+                    md["grid_mapping"] = global_attrs["grid_mapping"]
+                end
+            end     
         end
         md
     end
@@ -571,4 +578,4 @@ function _def_dim_var!(ds::AbstractDataset, dim::Dimension)
     end
     CDM.defVar(ds, dimname, Vector(index(dim)), (dimname,); attrib=attrib)
     return nothing
-end
\ No newline at end of file
+end
diff --git a/test/runtests.jl b/test/runtests.jl
index 06cad227..589bf472 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -22,11 +22,12 @@ end
 @time @safetestset "adapt" begin include("adapt.jl") end
 @time @safetestset "reproject" begin include("reproject.jl") end
 @time @safetestset "warp" begin include("warp.jl") end
-@time @safetestset "resample" begin include("resample.jl") end
 @time @safetestset "cellarea" begin include("cellarea.jl") end
 
 # CommondataModel sources
+@time @safetestset "commondatamodel" begin include("sources/commondatamodel.jl") end
 @time @safetestset "ncdatasets" begin include("sources/ncdatasets.jl") end
+# @time @safetestset "zarr" begin include("sources/zarr.jl") end # TODO: FIXME
 if !Sys.iswindows()
     # GRIBDatasets doesn't work on Windows for now
     @time @safetestset "gribdatasets" begin include("sources/gribdatasets.jl") end
@@ -43,3 +44,5 @@ if !Sys.iswindows()
     @time @safetestset "grd" begin include("sources/grd.jl") end
 end
 @time @safetestset "plot recipes" begin include("plotrecipes.jl") end
+
+@time @safetestset "resample" begin include("resample.jl") end
diff --git a/test/sources/commondatamodel.jl b/test/sources/commondatamodel.jl
index 451d18b3..6e246009 100644
--- a/test/sources/commondatamodel.jl
+++ b/test/sources/commondatamodel.jl
@@ -21,4 +21,25 @@ import Rasters: ForwardOrdered, ReverseOrdered, Regular
     @test steps[1] == 0.05
     @test steps[2] == 1/3
 
-end
\ No newline at end of file
+end
+
+@testset "faulty / no grid mapping" begin
+    # Construct a NetCDF dataset
+    filepath = joinpath(tempdir(), "test.nc")
+    ds = NCDataset(filepath, "c")
+    data = [Float32(i + j) for i = 1:100, j = 1:110]
+    v = defVar(ds, "temperature", data, ("lon", "lat"))
+    # Give `v` a "grid_mapping" attribute
+    # that points to a non-existent variable
+    v.attrib["grid_mapping"] = "non_existent_variable"
+    close(ds)
+    
+    # try loading raster
+    @test_nowarn Raster(filepath)
+    # make sure we don't magically materialize a CRS
+    @test isnothing(Rasters.crs(Raster(filepath)))
+
+    # Once Rasters has better CF-CRS support,
+    # we should be able to load a splatted global CRS,
+    # or a CRS as a global attribute a la Zarr.
+end
diff --git a/test/sources/zarr.jl b/test/sources/zarr.jl
index ef4c2d29..8075d0f9 100644
--- a/test/sources/zarr.jl
+++ b/test/sources/zarr.jl
@@ -1,5 +1,6 @@
-using Rasters, Zarr
+using Rasters
 using ZarrDatasets
+using ZarrDatasets.Zarr
 using Rasters: FileArray, FileStack, Zarrsource, crs, bounds, name, trim
 
 path = "https://s3.bgc-jena.mpg.de:9000/esdl-esdc-v3.0.2/esdc-16d-2.5deg-46x72x1440-3.0.2.zarr"

From 1f00143fe722c27e2fb25899fc077212de7bc79f Mon Sep 17 00:00:00 2001
From: Tiem van der Deure <tiemvanderdeure@gmail.com>
Date: Wed, 16 Oct 2024 14:25:01 +0200
Subject: [PATCH 5/9] Add `Rasters.sample` (#791)

* extract tweaks

* add rasters.sample

* get rid of some toy code

* just write kw... in docs

Co-authored-by: Rafael Schouten <rafaelschouten@gmail.com>

* use Type{G}

Co-authored-by: Rafael Schouten <rafaelschouten@gmail.com>

* add a missing where G

* neater formatting

* better reuse _rowtypes code

* more code style

Co-authored-by: Rafael Schouten <rafaelschouten@gmail.com>

* extract tweaks

* add rasters.sample

* get rid of some toy code

* neater formatting

* better reuse _rowtypes code

* more code style

Co-authored-by: Rafael Schouten <rafaelschouten@gmail.com>

* merge _rowtype changes from other branch

* optimization for extract points with skipmissing

* streamline _rowtype

* correctly use _rowtype

* improve type stability

* switch to using dimindices internally

* let users specify geometry type

* just use a broadcast

* don't need to be compatible with julia < 1.9

Co-authored-by: Rafael Schouten <rafaelschouten@gmail.com>

* deal with type stability better

* allow not providing `n`

* tweak the docstring

---------

Co-authored-by: Rafael Schouten <rafaelschouten@gmail.com>
---
 Project.toml                                  |  7 +-
 .../RastersStatsBaseExt.jl                    | 13 +++
 ext/RastersStatsBaseExt/sample.jl             | 94 +++++++++++++++++++
 src/extensions.jl                             | 44 +++++++++
 src/methods/extract.jl                        | 43 ++-------
 src/utils.jl                                  | 12 ++-
 test/methods.jl                               | 80 ++++++++++++++++
 7 files changed, 255 insertions(+), 38 deletions(-)
 create mode 100644 ext/RastersStatsBaseExt/RastersStatsBaseExt.jl
 create mode 100644 ext/RastersStatsBaseExt/sample.jl

diff --git a/Project.toml b/Project.toml
index c6db47d0..a7ba6b71 100644
--- a/Project.toml
+++ b/Project.toml
@@ -34,6 +34,7 @@ Makie = "ee78f7c6-11fb-53f2-987a-cfe4a2b5a57a"
 NCDatasets = "85f8d34a-cbdd-5861-8df4-14fed0d494ab"
 Proj = "c94c279d-25a6-4763-9509-64d165bea63e"
 RasterDataSources = "3cb90ccd-e1b6-4867-9617-4276c8b2ca36"
+StatsBase = "2913bbd2-ae8a-5f71-8c99-4fb6c76f3a91"
 ZarrDatasets = "519a4cdf-1362-424a-9ea1-b1d782dbb24b"
 
 [extensions]
@@ -44,6 +45,7 @@ RastersMakieExt = "Makie"
 RastersNCDatasetsExt = "NCDatasets"
 RastersProjExt = "Proj"
 RastersRasterDataSourcesExt = "RasterDataSources"
+RastersStatsBaseExt = "StatsBase"
 RastersZarrDatasetsExt = "ZarrDatasets"
 
 [compat]
@@ -81,6 +83,7 @@ SafeTestsets = "0.1"
 Setfield = "0.6, 0.7, 0.8, 1"
 Shapefile = "0.10, 0.11"
 Statistics = "1"
+StatsBase = "0.34"
 Test = "1"
 ZarrDatasets = "0.1"
 julia = "1.10"
@@ -101,8 +104,10 @@ Proj = "c94c279d-25a6-4763-9509-64d165bea63e"
 RasterDataSources = "3cb90ccd-e1b6-4867-9617-4276c8b2ca36"
 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", "Statistics", "Test", "ZarrDatasets"]
+test = ["Aqua", "ArchGDAL", "CFTime", "CoordinateTransformations", "DataFrames", "GeoDataFrames", "GeometryBasics", "GRIBDatasets", "NCDatasets", "Plots", "Proj", "RasterDataSources", "SafeTestsets", "Shapefile", "StableRNGs", "Statistics", "StatsBase", "Test", "ZarrDatasets"]
diff --git a/ext/RastersStatsBaseExt/RastersStatsBaseExt.jl b/ext/RastersStatsBaseExt/RastersStatsBaseExt.jl
new file mode 100644
index 00000000..fd6f1e91
--- /dev/null
+++ b/ext/RastersStatsBaseExt/RastersStatsBaseExt.jl
@@ -0,0 +1,13 @@
+module RastersStatsBaseExt
+
+using Rasters, StatsBase
+using StatsBase.Random
+
+const RA = Rasters
+
+import Rasters: _True, _False, _booltype, istrue
+import Rasters.DimensionalData as DD
+
+include("sample.jl")
+
+end # Module
diff --git a/ext/RastersStatsBaseExt/sample.jl b/ext/RastersStatsBaseExt/sample.jl
new file mode 100644
index 00000000..f54f9d91
--- /dev/null
+++ b/ext/RastersStatsBaseExt/sample.jl
@@ -0,0 +1,94 @@
+Rasters.sample(x::RA.RasterStackOrArray, args...; kw...) = Rasters.sample(Random.GLOBAL_RNG, x, args...; kw...)
+@inline function Rasters.sample(
+    rng::Random.AbstractRNG, x::RA.RasterStackOrArray, args...; 
+    geometry=(X,Y),
+    index = false, 
+    names=RA._names(x), 
+    name=names, 
+    skipmissing=false,
+    weights=nothing, 
+    weightstype::Type{<:StatsBase.AbstractWeights}=StatsBase.Weights,
+    kw...
+)
+    na = DD._astuple(name)
+    geometry, geometrytype, dims = _geometrytype(x, geometry)
+
+    _sample(rng, x, args...;
+        dims,
+        names = NamedTuple{na}(na),
+        geometry,
+        geometrytype,  
+        # These keywords are converted to _True/_False for type stability later on
+        index = _booltype(index), 
+        skipmissing = _booltype(skipmissing), 
+        weights,
+        weightstype,
+        kw... # passed to StatsBase.sample, could be replace or ordered
+    )
+end
+function _sample(
+    rng, x, n::Integer;
+    dims, names::NamedTuple{K}, geometry, geometrytype::Type{G}, index, skipmissing, weights, weightstype, kw..., 
+) where {K, G}
+    indices = sample_indices(rng, x, skipmissing, weights, weightstype, n; kw...)
+    T = RA._rowtype(x, G; geometry, index, skipmissing, skipinvalid = _True(), names)
+    x2 = x isa AbstractRasterStack ? x[K] : RasterStack(NamedTuple{K}((x,)))
+    return _getindices(T, x2, dims, indices)
+end
+function _sample(
+    rng, x;
+    dims, names::NamedTuple{K}, geometry, geometrytype::Type{G}, index, skipmissing, weights, weightstype, kw..., 
+) where {K, G}
+    indices = sample_indices(rng, x, skipmissing, weights, weightstype)
+    T = RA._rowtype(x, G; geometry, index, skipmissing, skipinvalid = _True(), names)
+    x2 = x isa AbstractRasterStack ? x[K] : RasterStack(NamedTuple{K}((x,)))
+    return _getindex(T, x2, dims, indices)
+end
+
+_getindices(::Type{T}, x, dims, indices) where {T} = 
+    broadcast(I -> _getindex(T, x, dims, I), indices)
+
+function _getindex(::Type{T}, x::AbstractRasterStack{<:Any, NT}, dims, idx) where {T, NT}
+    RA._maybe_add_fields(
+        T, 
+        NT(x[RA.commondims(idx, x)]), 
+        DimPoints(dims)[RA.commondims(idx, dims)], 
+        val(idx)
+    )
+end
+
+# args may be an integer or nothing
+function sample_indices(rng, x, skipmissing, weights::Nothing, weightstype, args...; kw...)
+    if istrue(skipmissing)
+        wts = StatsBase.Weights(vec(boolmask(x)))
+        StatsBase.sample(rng, RA.DimIndices(x), wts, args...; kw...)
+    else
+        StatsBase.sample(rng, RA.DimIndices(x), args...; kw...)
+    end
+end
+function sample_indices(rng, x, skipmissing, weights::AbstractDimArray, ::Type{W}, args...; kw...) where W
+    wts = if istrue(skipmissing) 
+        @d boolmask(x) .* weights
+    elseif dims(weights) == dims(x)
+        weights
+    else
+        @d ones(eltype(weights), dims(x)) .* weights
+    end |> vec |> W
+    StatsBase.sample(rng, RA.DimIndices(x), wts, args...; kw...)
+end
+function _geometrytype(x, geometry::Bool)
+    if geometry
+        error("Specify a geometry type by setting `geometry` to a Tuple or NamedTuple of Dimensions. E.g. `geometry = (X, Y)`")
+    else
+        return _False(), Nothing, dims(x)
+    end
+end
+
+function _geometrytype(x, geometry::Tuple)
+    dims = DD.commondims(DD.dims(x), geometry)
+    return _True(), Tuple{map(eltype, dims)...}, dims
+end
+function _geometrytype(x, geometry::NamedTuple{K}) where K
+    dims = DD.commondims(DD.dims(x), values(geometry))
+    return _True(), NamedTuple{K, Tuple{map(eltype, dims)...}}, dims
+end
\ No newline at end of file
diff --git a/src/extensions.jl b/src/extensions.jl
index f56f5bd5..de579bde 100644
--- a/src/extensions.jl
+++ b/src/extensions.jl
@@ -228,6 +228,50 @@ Note that cellarea returns the area in square m, while cellsize still uses squar
     return cellarea(args...; kw..., radius = 6371.0088)
 end
 
+"""
+Rasters.sample([rng], x, [n::Integer]; kw...)
+
+Sample `n` random and optionally weighted points from from a `Raster` or `RasterStack`.
+Returns a `Vector` of `NamedTuple`, closely resembling the return type of [`extract`](@ref).
+
+Run `using StatsBase` to make this method available.
+Note that this function is not exported to avoid confusion with StatsBase.sample
+
+# Keywords
+
+- `geometry`: include `:geometry` in returned `NamedTuple`. Specify the type and dimensions of the returned geometry by
+    providing a `Tuple` or `NamedTuple` of dimensions. Defaults to `(X,Y)`
+- `index`: include `:index` of the `CartesianIndex` in returned `NamedTuple`, `false` by default.
+- `name`: a `Symbol` or `Tuple` of `Symbol` corresponding to layer/s of a `RasterStack` to extract. All layers by default.
+- `skipmissing`: skip missing points automatically.
+- `weights`: A DimArray that matches one or more of the dimensions of `x` with weights for sampling.
+- `weightstype`: a `StatsBase.AbstractWeights` specifying the type of weights. Defaults to `StatsBase.Weights`.
+- `replace`: sample with replacement, `true` by default. See `StatsBase.sample`
+- `ordered`: sample in order, `false` by default. See `StatsBase.sample`
+
+# Example
+This code draws 5 random points from a raster, weighted by cell area.
+```julia
+using Rasters, Rasters.Lookups, Proj, StatsBase
+xdim = X(Projected(90.0:10.0:120; sampling=Intervals(Start()), crs=EPSG(4326)))
+ydim = Y(Projected(0.0:10.0:50; sampling=Intervals(Start()), crs=EPSG(4326)))
+myraster = rand(xdim, ydim)
+Rasters.sample(myraster, 5; weights = cellarea(myraster))
+
+# output
+
+5-element Vector{@NamedTuple{geometry::Tuple{Float64, Float64}, ::Union{Missing, Float64}}}:
+ @NamedTuple{geometry::Tuple{Float64, Float64}, ::Union{Missing, Float64}}(((90.0, 10.0), 0.7360504790189618))
+ @NamedTuple{geometry::Tuple{Float64, Float64}, ::Union{Missing, Float64}}(((90.0, 30.0), 0.5447657183842469))
+ @NamedTuple{geometry::Tuple{Float64, Float64}, ::Union{Missing, Float64}}(((90.0, 30.0), 0.5447657183842469))
+ @NamedTuple{geometry::Tuple{Float64, Float64}, ::Union{Missing, Float64}}(((90.0, 10.0), 0.7360504790189618))
+ @NamedTuple{geometry::Tuple{Float64, Float64}, ::Union{Missing, Float64}}(((110.0, 10.0), 0.5291143028176258))
+```
+"""
+sample(args...; kw...) = throw_extension_error(sample, "StatsBase", :RastersStatsBaseExt, args)
+
+
+
 # Other shared stubs
 function layerkeys end
 function smapseries end
diff --git a/src/methods/extract.jl b/src/methods/extract.jl
index 5fd748ee..ef41b04e 100644
--- a/src/methods/extract.jl
+++ b/src/methods/extract.jl
@@ -72,7 +72,7 @@ function extract end
 end
 
 function _extract(A::RasterStackOrArray, geom::Missing, names, kw...)
-    T = _extractrowtype(A, geom; names, kw...)
+    T = _rowtype(A, geom; names, kw...)
     [_maybe_add_fields(T, map(_ -> missing, names), missing, missing)]
 end
 function _extract(A::RasterStackOrArray, geom; names, kw...)
@@ -82,11 +82,7 @@ function _extract(A::RasterStackOrArray, ::Nothing, data;
     names, skipmissing, geometrycolumn, kw...
 )
     geoms = _get_geometries(data, geometrycolumn)
-    T = if istrue(skipmissing)
-        _extractrowtype(A, nonmissingtype(eltype(geoms)); names, skipmissing, kw...)
-    else
-        _extractrowtype(A, eltype(geoms); names, skipmissing, kw...)
-    end
+    T = _rowtype(A, eltype(geoms); names, skipmissing, kw...)
     # Handle empty / all missing cases
     (length(geoms) > 0 && any(!ismissing, geoms)) || return T[]
     
@@ -95,15 +91,13 @@ function _extract(A::RasterStackOrArray, ::Nothing, data;
     # We need to split out points from other geoms
     # TODO this will fail with mixed point/geom vectors
     if trait1 isa GI.PointTrait
+    rows = Vector{T}(undef, length(geoms))
         if istrue(skipmissing)
-            T2 = _extractrowtype(A, eltype(geoms), Tuple{Int64, Int64}; 
-                names, skipmissing = _False(), kw...)
-            rows = Vector{T2}(undef, length(geoms))
             j = 1
             for i in eachindex(geoms)
                 g = geoms[i]
                 ismissing(g) && continue
-                e = _extract_point(T2, A, g, skipmissing; names, kw...)
+                e = _extract_point(T, A, g, skipmissing; names, kw...)
                 if !ismissing(e) 
                     rows[j] = e
                     j += 1
@@ -111,9 +105,7 @@ function _extract(A::RasterStackOrArray, ::Nothing, data;
                 nothing
             end
             deleteat!(rows, j:length(rows))
-            rows = T === T2 ? rows : T.(rows)
         else
-            rows = Vector{T}(undef, length(geoms))
             for i in eachindex(geoms)
                 g = geoms[i]
                 rows[i] = _extract_point(T, A, g, skipmissing; names, kw...)::T
@@ -141,14 +133,14 @@ end
 function _extract(A::RasterStackOrArray, ::GI.AbstractMultiPointTrait, geom; 
     skipmissing, kw...
 )
-    T = _extractrowtype(A, GI.getpoint(geom, 1); names, skipmissing, kw...)
+    T = _rowtype(A, GI.getpoint(geom, 1); names, skipmissing, kw...)
     rows = (_extract_point(T, A, p, skipmissing; kw...) for p in GI.getpoint(geom))
     return skipmissing isa _True ? collect(_skip_missing_rows(rows, _missingval_or_missing(A), names)) : collect(rows)
 end
 function _extract(A::RasterStackOrArray, ::GI.PointTrait, geom; 
     skipmissing, kw...
 )
-    T = _extractrowtype(A, geom; names, skipmissing, kw...)
+    T = _rowtype(A, geom; names, skipmissing, kw...)
     _extract_point(T, A, geom, skipmissing; kw...)
 end
 function _extract(A::RasterStackOrArray, t::GI.AbstractGeometryTrait, geom;
@@ -166,7 +158,7 @@ function _extract(A::RasterStackOrArray, t::GI.AbstractGeometryTrait, geom;
     else
         GI.x(p), GI.y(p)
     end
-    T = _extractrowtype(A, tuplepoint; names, skipmissing, kw...)
+    T = _rowtype(A, tuplepoint; names, skipmissing, kw...)
     B = boolmask(geom; to=template, kw...)
     offset = CartesianIndex(map(x -> first(x) - 1, parentindices(parent(template))))
     # Add a row for each pixel that is `true` in the mask
@@ -285,27 +277,8 @@ Base.@assume_effects :total function _maybe_add_fields(::Type{T}, props::NamedTu
         :index in K ? merge((; geometry=point, index=I), props) : merge((; geometry=point), props)
     else
         :index in K ? merge((; index=I), props) : props
-    end
-end
-
-function _extractrowtype(x, g; geometry, index, skipmissing, names, kw...)
-    I = if istrue(skipmissing)
-        Tuple{Int, Int}
-    else
-        Union{Missing, Tuple{Int, Int}}
-    end
-    _extractrowtype(x, g, I; geometry, index, skipmissing, names, kw...)
-end
-function _extractrowtype(x, g, ::Type{I}; geometry, index, skipmissing, names, kw...) where I
-    G = if istrue(skipmissing)
-        nonmissingtype(typeof(g))
-    else
-        typeof(g)
-    end
-    _extractrowtype(x, G, I; geometry, index, skipmissing, names, kw...)
+    end |> T
 end
-_extractrowtype(x, ::Type{G}, ::Type{I}; geometry, index, skipmissing, names, kw...) where {G, I} =
-    _rowtype(x, G, I; geometry, index, skipmissing, names)
 
 @inline _skip_missing_rows(rows, ::Missing, names) = 
     Iterators.filter(row -> !any(ismissing, row), rows)
diff --git a/src/utils.jl b/src/utils.jl
index 67a8d2c4..5fe150cc 100644
--- a/src/utils.jl
+++ b/src/utils.jl
@@ -372,12 +372,20 @@ _booltype(x) = x ? _True() : _False()
 istrue(::_True) = true
 istrue(::_False) = false
 
-function _rowtype(x, ::Type{G}, ::Type{I}; geometry, index, skipmissing, names) where {G, I}
+# skipinvalid: can G and I be missing. skipmissing: can nametypes be missing
+_rowtype(x, g, args...; kw...) = _rowtype(x, typeof(g), args...; kw...)
+function _rowtype(
+    x, ::Type{G}, i::Type{I} = typeof(size(x)); 
+    geometry, index, skipmissing, skipinvalid = skipmissing, names, kw...
+) where {G, I}
+    _G = istrue(skipinvalid) ? nonmissingtype(G) : G
+    _I = istrue(skipinvalid) ? I : Union{Missing, I}
     keys = _rowkeys(geometry, index, names)
-    types = _rowtypes(x, G, I, geometry, index, skipmissing, names)
+    types = _rowtypes(x, _G, _I, geometry, index, skipmissing, names)
     NamedTuple{keys,types}
 end
 
+
 function _rowtypes(
     x, ::Type{G}, ::Type{I}, geometry::_True, index::_True, skipmissing, names::NamedTuple{Names}
 ) where {G,I,Names}
diff --git a/test/methods.jl b/test/methods.jl
index df1a81e5..1b48f381 100644
--- a/test/methods.jl
+++ b/test/methods.jl
@@ -671,3 +671,83 @@ end
                    0.0     0.0     missing    
                    missing missing missing], dims=3))
 end
+
+using StableRNGs, StatsBase
+test = rebuild(ga; name = :test)
+@testset "sample" begin
+    # test that all keywords work and return the same thing as extract
+    @test all(Rasters.sample(StableRNG(123), test, 2) .=== extract(test, [(2.0,2.0), (1.0,2.0)]))
+    @test all(Rasters.sample(StableRNG(123), st2, 2) .=== extract(st2, [(2,2), (1,2)]))
+    @test all(Rasters.sample(StableRNG(123), test, 2; geometry = false) .=== extract(test, [(2.0,2.0), (1.0,2.0)]; geometry = false))
+    @test all(Rasters.sample(StableRNG(123), st2, 2; geometry = false) .=== extract(st2, [(2,2), (1,2)]; geometry = false))
+    
+    @test all(Rasters.sample(StableRNG(123), test, 2, skipmissing = true) .=== extract(test, [(2.0,1.0), (2.0,1.0)], skipmissing = true))
+    @test all(Rasters.sample(StableRNG(123), st2, 2, skipmissing = true) .=== extract(st2, [(2,1), (2,1)], skipmissing = true))
+
+    @test all(
+        Rasters.sample(StableRNG(123), test, 2, skipmissing = true, index = true) .=== 
+        extract(test, [(2.0,1.0), (2.0,1.0)], skipmissing = true, index = true))
+    @test all(
+        Rasters.sample(StableRNG(123), st2, 2, skipmissing = true, index = true) .=== 
+        extract(st2, [(2,1), (2,1)], skipmissing = true, index = true))
+
+    kws = [(:geometry => false,), (:index => true,), ()]
+    for kw in kws
+        @test all(        
+            Rasters.sample(StableRNG(123), test, 2; skipmissing = true, kw...) .=== 
+            extract(test, [(2.0,1.0), (2.0,1.0)]; skipmissing = true, kw...)
+            )
+            
+        @test all(
+            Rasters.sample(StableRNG(123), st2, 2; skipmissing = true, kw...) .== 
+            extract(st2, [(2,1), (2,1)]; skipmissing = true, kw...)
+        )
+    end
+    @test all(Rasters.sample(StableRNG(123), st2, 2, name = (:a,)) .=== extract(st2, [(2,2), (1,2)], name = (:a,)))
+
+    # in this case extract and sample always return different types
+    @test eltype(Rasters.sample(StableRNG(123), test, 2, index = true)) != eltype(extract(test, [(2.0,1.0), (2.0,1.0)], index = true))
+    @test eltype(Rasters.sample(StableRNG(123), st2, 2, index = true)) != eltype(extract(st2, [(2,1), (2,1)], index = true))
+
+    @test all(
+        Rasters.sample(StableRNG(123), test, 2, weights = DimArray([1,1000], X(1:2)), skipmissing = true) .===
+        [
+            (geometry = (2.0,1.0), test = 2.0f0)
+            (geometry = (2.0,1.0), test = 2.0f0)
+        ]
+    )
+
+    @test all(
+        Rasters.sample(StableRNG(123), test, 2, weights = DimArray([1,1000], X(1:2)), skipmissing = true, weightstype = StatsBase.FrequencyWeights) .===
+        [
+            (geometry = (2.0,1.0), test = 2.0f0)
+            (geometry = (2.0,1.0), test = 2.0f0)
+        ]
+    )
+
+    @test all(
+        Rasters.sample(StableRNG(123), test, 2, skipmissing = true, replace = false) .===
+        [
+            (geometry = (2.0,1.0), test = 2.0f0)
+            (geometry = (1.0,2.0), test = 7.0f0)
+        ]
+    )
+    @test all(
+        Rasters.sample(StableRNG(123), test, 2, skipmissing = true, replace = false, ordered = true) .===
+        [
+            (geometry = (2.0,1.0), test = 2.0f0)
+            (geometry = (1.0,2.0), test = 7.0f0)
+        ]
+    )
+
+    # test providing a geometry type works
+    @test typeof(first(
+        Rasters.sample(StableRNG(123), test, 2, geometry = (X = X, Y = Y))
+    ).geometry) <: NamedTuple{(:X, :Y)}
+
+    # test this works without providing n
+    @test Rasters.sample(test, geometry = (X = X, Y = Y)) isa NamedTuple{(:geometry, :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
\ No newline at end of file

From fc27e99953fbc2bc325a2e15de095eb917fa7c59 Mon Sep 17 00:00:00 2001
From: Tiem van der Deure <tiemvanderdeure@gmail.com>
Date: Thu, 24 Oct 2024 17:13:01 +0200
Subject: [PATCH 6/9] bugfix to _without_mapped_crs (#810)

* bugfix

* add tests
---
 src/utils.jl       |  2 +-
 test/methods.jl    | 87 ++++++++++++++++++++++++----------------------
 test/test_utils.jl |  8 ++++-
 3 files changed, 54 insertions(+), 43 deletions(-)

diff --git a/src/utils.jl b/src/utils.jl
index 5fe150cc..e069755c 100644
--- a/src/utils.jl
+++ b/src/utils.jl
@@ -132,7 +132,7 @@ function _without_mapped_crs(f, st::AbstractRasterStack, mappedcrs::GeoFormat)
     st1 = map(A -> setmappedcrs(A, nothing), st)
     x = f(st1)
     if x isa AbstractRasterStack
-        x = map(A -> setmappedcrs(A, mappedcrs(st)), x)
+        x = map(A -> setmappedcrs(A, mappedcrs), x)
     end
     return x
 end
diff --git a/test/methods.jl b/test/methods.jl
index 1b48f381..324edaa8 100644
--- a/test/methods.jl
+++ b/test/methods.jl
@@ -556,11 +556,16 @@ end
          missing 1.0     missing]
 
     r_fwd = Raster(A, (X(1.0:1.0:3.0), Y(1.0:1.0:3.0)); missingval=missing)
+    r_crs = Raster(A, (X(1.0:1.0:3.0), Y(1.0:1.0:3.0)); missingval=missing, crs = EPSG(4326))
+    r_mcrs = Raster(A, (X(1.0:1.0:3.0), Y(1.0:1.0:3.0)); missingval=missing, crs = EPSG(4326), mappedcrs = EPSG(4326))
     r_revX = reverse(r_fwd; dims=X)
     r_revY = reverse(r_fwd; dims=Y)
+    st1 = RasterStack((a = r_fwd, b = r_fwd))
+    st2 = RasterStack((a = r_fwd, b = r_fwd), crs = EPSG(4326))
+    st3 = RasterStack((a = r_fwd, b = r_fwd), crs = EPSG(4326), mappedcrs = EPSG(4326))
 
-    ga = r_revY
-    for ga in (r_fwd, r_revX, r_revY)
+    ga = r_fwd
+    for ga in (r_fwd, r_crs, r_mcrs, r_revX, r_revY, st1, st2, st3)
         # Test with missing on all sides
         ga_r = rot180(ga)
         trimmed = trim(ga)
@@ -574,55 +579,55 @@ end
         cropped1 = crop(crop(ga; to=dims(trimmed, X)); to=dims(trimmed, Y))
         _, cropped2 = crop(trimmed, ga)
         cropped_r = crop(ga_r; to=trimmed_r)
-        @test all(cropped .=== cropped1 .=== trimmed)
-        @test all(cropped_r .=== trimmed_r)
+        @test map(identicalelements, maybe_layers.((cropped, cropped1, trimmed))...) |> all
+        @test map(identicalelements, maybe_layers.((cropped_r, trimmed_r))...) |> all
         extended = extend(cropped, ga)[1]
         extended_r = extend(cropped_r; to=ga_r)
         extended1 = extend(extend(cropped; to=dims(ga, X)); to=dims(ga, Y))
         extended_d = extend(cropped; to=ga, filename="extended.tif")
-        @test all(extended .=== extended1 .=== replace_missing(extended_d) .=== ga) 
-        @test all(extended_r .=== ga_r)
+        @test map(identicalelements, maybe_layers.((extended, extended1, replace_missing(extended_d), ga))...) |> all
+        @test map(identicalelements, maybe_layers.((extended_r, ga_r))...) |> all
         @test all(map(==, lookup(extended_d), lookup(extended)))
+    end
 
-        @testset "unformatted dimension works in crop" begin
-            xdim, ydim = X(1.0:0.2:2.0), Y(1.0:1:2.0)
-            A = Raster(rand((X(1.0:0.2:4.0), ydim)))
-            @test lookup(crop(A; to=xdim), X) == 1.0:0.2:2.0
+    @testset "unformatted dimension works in crop" begin
+        xdim, ydim = X(1.0:0.2:2.0), Y(1.0:1:2.0)
+        A = Raster(rand((X(1.0:0.2:4.0), ydim)))
+        @test lookup(crop(A; to=xdim), X) == 1.0:0.2:2.0
 
-            A = Raster(rand((X(1.0:0.2:4.0), ydim)))
-            @test lookup(crop(A; to=xdim), X) == 1.0:0.2:2.0
+        A = Raster(rand((X(1.0:0.2:4.0), ydim)))
+        @test lookup(crop(A; to=xdim), X) == 1.0:0.2:2.0
 
-            A = Raster(ones((X(1.0:0.2:1.4), ydim)); missingval=0.0)
-            extend(A; to=xdim)
-            @test lookup(extend(A; to=xdim), X) == 1.0:0.2:2.0
-            @test extend(A; to=xdim) == [1.0 1.0; 1.0 1.0; 1.0 1.0; 0.0 0.0; 0.0 0.0; 0.0 0.0]  
-        end
+        A = Raster(ones((X(1.0:0.2:1.4), ydim)); missingval=0.0)
+        extend(A; to=xdim)
+        @test lookup(extend(A; to=xdim), X) == 1.0:0.2:2.0
+        @test extend(A; to=xdim) == [1.0 1.0; 1.0 1.0; 1.0 1.0; 0.0 0.0; 0.0 0.0; 0.0 0.0]  
+    end
 
-        @testset "to polygons" begin
-            A1 = Raster(zeros(X(-20:-5; sampling=Points()), Y(0:30; sampling=Points())))
-            A2 = Raster(ones(X(-20:-5; sampling=Points()), Y(0:30; sampling=Points())))
-            A1crop1 = crop(A1; to=polygon)
-            A1crop2, A2crop = crop(A1, A2; to=polygon)
-            size(A1crop1)
-            @test size(A1crop1) == size(A1crop2) == size(A2crop) == (16, 21)
-            @test bounds(A1crop1) == bounds(A1crop2) == bounds(A2crop) == ((-20, -5), (10, 30))
-            A1extend1 = extend(A1; to=polygon)
-            A1extend2, A2extend = extend(A1, A2; to=polygon)
-            @test all(A1extend1 .=== A1extend2)
-            @test size(A1extend1) == size(A1extend2) == size(A2extend) == (21, 31)
-            @test bounds(A1extend1) == bounds(A1extend2) == bounds(A2extend) == ((-20.0, 0.0), (0, 30))
-        end
-        @testset "to featurecollection and table" begin
-            A1 = Raster(zeros(X(-20:-5; sampling=Points()), Y(0:30; sampling=Points())))
-            featurecollection = map(GeoInterface.getpoint(polygon), vals) do geometry, v
-                (; geometry, val1=v, val2=2.0f0v)
-            end
-            fccrop = crop(A1; to=featurecollection)
-            table = DataFrame(featurecollection)
-            tablecrop = crop(A1; to=table)
-            @test size(fccrop) == size(tablecrop) == (16, 21)
-            @test bounds(fccrop) == bounds(tablecrop) == ((-20, -5), (10, 30))
+    @testset "to polygons" begin
+        A1 = Raster(zeros(X(-20:-5; sampling=Points()), Y(0:30; sampling=Points())))
+        A2 = Raster(ones(X(-20:-5; sampling=Points()), Y(0:30; sampling=Points())))
+        A1crop1 = crop(A1; to=polygon)
+        A1crop2, A2crop = crop(A1, A2; to=polygon)
+        size(A1crop1)
+        @test size(A1crop1) == size(A1crop2) == size(A2crop) == (16, 21)
+        @test bounds(A1crop1) == bounds(A1crop2) == bounds(A2crop) == ((-20, -5), (10, 30))
+        A1extend1 = extend(A1; to=polygon)
+        A1extend2, A2extend = extend(A1, A2; to=polygon)
+        @test all(A1extend1 .=== A1extend2)
+        @test size(A1extend1) == size(A1extend2) == size(A2extend) == (21, 31)
+        @test bounds(A1extend1) == bounds(A1extend2) == bounds(A2extend) == ((-20.0, 0.0), (0, 30))
+    end
+    @testset "to featurecollection and table" begin
+        A1 = Raster(zeros(X(-20:-5; sampling=Points()), Y(0:30; sampling=Points())))
+        featurecollection = map(GeoInterface.getpoint(polygon), vals) do geometry, v
+            (; geometry, val1=v, val2=2.0f0v)
         end
+        fccrop = crop(A1; to=featurecollection)
+        table = DataFrame(featurecollection)
+        tablecrop = crop(A1; to=table)
+        @test size(fccrop) == size(tablecrop) == (16, 21)
+        @test bounds(fccrop) == bounds(tablecrop) == ((-20, -5), (10, 30))
     end
 
 end
diff --git a/test/test_utils.jl b/test/test_utils.jl
index 1557dfc7..688f993d 100644
--- a/test/test_utils.jl
+++ b/test/test_utils.jl
@@ -20,4 +20,10 @@ function temporary_random_rasters(f, N, size, type=UInt8)
     finally
         rm.(filenames; force=true)
     end
-end
\ No newline at end of file
+end
+
+maybe_layers(r::AbstractRaster) = (r,)
+maybe_layers(r::AbstractRasterStack) = layers(r)
+
+identicalelements(x,y,args...) = all(x .=== y) && identicalelements(x, args...)
+identicalelements(x,y) = all(x .=== y)
\ No newline at end of file

From 93bf8eaa41a5d105b9548526179e7fb2814461f6 Mon Sep 17 00:00:00 2001
From: Rafael Schouten <rafaelschouten@gmail.com>
Date: Sat, 2 Nov 2024 14:02:28 +0100
Subject: [PATCH 7/9] Update issue templates

---
 .github/ISSUE_TEMPLATE/bug_report.md | 32 ++++++++++++++++++++++++++++
 1 file changed, 32 insertions(+)
 create mode 100644 .github/ISSUE_TEMPLATE/bug_report.md

diff --git a/.github/ISSUE_TEMPLATE/bug_report.md b/.github/ISSUE_TEMPLATE/bug_report.md
new file mode 100644
index 00000000..6ae6610e
--- /dev/null
+++ b/.github/ISSUE_TEMPLATE/bug_report.md
@@ -0,0 +1,32 @@
+---
+name: Bug report
+about: Describe your issue here, including a minimum working example all files.
+title: ''
+labels: bug
+assignees: ''
+
+---
+
+Raster data is complicated, there are a lot of formats and a lot of things that can go wrong. So we can't usually fix your problem unless you help us by filling out the code blocks below. 
+
+Its important that your example just runs, and we can copy and paste the code into Julia and it does that same thing for us that it does for you. That means either generating or downloading data inside the script, e.g. with `download`. 
+
+The faster we can get to your problem, the more likely it will be fixed quickly. If we have 15 minutes of free dev time for your issue, you want all 15 to be looking at code, not 10 trying to download some file from a website and making sure the file path matches in the script, leaving 5 minutes to look at the error and fix it.
+
+```julia
+All the julia code needed to make the error goes here including download or generation of all data used
+
+This is the most important thing in the issue, we cant fix bugs without it.
+```
+
+```julia
+The full complete error from start to finish goes here
+
+This is the second most important thing in the issue, we can't be sure we fixed the same problem unless we can see this.
+```
+
+```julia
+Paste the results of `versioninfo()` here.
+
+Paste the results of `] status Rasters` here.
+```

From 58a3a57d8b043b149188c9c24254023d1288e799 Mon Sep 17 00:00:00 2001
From: Rafael Schouten <rafaelschouten@gmail.com>
Date: Sat, 2 Nov 2024 14:03:30 +0100
Subject: [PATCH 8/9] Delete .github/ISSUE_TEMPLATE/bug-report.md

---
 .github/ISSUE_TEMPLATE/bug-report.md | 30 ----------------------------
 1 file changed, 30 deletions(-)
 delete mode 100644 .github/ISSUE_TEMPLATE/bug-report.md

diff --git a/.github/ISSUE_TEMPLATE/bug-report.md b/.github/ISSUE_TEMPLATE/bug-report.md
deleted file mode 100644
index 9551724b..00000000
--- a/.github/ISSUE_TEMPLATE/bug-report.md
+++ /dev/null
@@ -1,30 +0,0 @@
----
-name: Bug Report
-about: Describe your issue here, including a minimum working example and all file downloads without authentication.
-We cant fix your problem unless you help us by filling out the code blocks below.
-
-title: ''
-labels: ''
-assignees: ''
-
----
-
-Your problem. You don't have to write anything here for a bug, the code and error are what is important.
-
-```julia
-All the julia code needed to make the error goes here including download or generation of all data used
-
-This is the most important thing in the issue, we cant fix bugs without it.
-```
-
-```julia
-The full complete error from start to finish goes here
-
-This is the second most important thing in the issue, we can't be sure we fixed the same problem unless we can see this.
-```
-
-```julia
-Paste the results of `versioninfo()` here.
-
-Paste the results of `] status Rasters` here.
-```

From 251bbbdd7bfca54c22e582ee80d7b04d141de79e Mon Sep 17 00:00:00 2001
From: Rafael Schouten <rafaelschouten@gmail.com>
Date: Sat, 2 Nov 2024 14:07:49 +0100
Subject: [PATCH 9/9] Update issue templates

---
 .github/ISSUE_TEMPLATE/bug_report.md | 8 +++++---
 1 file changed, 5 insertions(+), 3 deletions(-)

diff --git a/.github/ISSUE_TEMPLATE/bug_report.md b/.github/ISSUE_TEMPLATE/bug_report.md
index 6ae6610e..dede49a5 100644
--- a/.github/ISSUE_TEMPLATE/bug_report.md
+++ b/.github/ISSUE_TEMPLATE/bug_report.md
@@ -7,11 +7,13 @@ assignees: ''
 
 ---
 
-Raster data is complicated, there are a lot of formats and a lot of things that can go wrong. So we can't usually fix your problem unless you help us by filling out the code blocks below. 
+########################
+Pleas fill all code blocks below, then read and remove this message.
 
-Its important that your example just runs, and we can copy and paste the code into Julia and it does that same thing for us that it does for you. That means either generating or downloading data inside the script, e.g. with `download`. 
+Its important that your example *just runs*, and we can copy and paste the code into Julia and it does that same thing for us that it does for you. 
 
-The faster we can get to your problem, the more likely it will be fixed quickly. If we have 15 minutes of free dev time for your issue, you want all 15 to be looking at code, not 10 trying to download some file from a website and making sure the file path matches in the script, leaving 5 minutes to look at the error and fix it.
+That means either generating or downloading data inside the script, e.g. with `download`. 
+########################
 
 ```julia
 All the julia code needed to make the error goes here including download or generation of all data used