From 2adfc0988e67986787ad300795e5cc64958a602c Mon Sep 17 00:00:00 2001
From: Rafael Schouten <rafaelschouten@gmail.com>
Date: Sun, 18 Aug 2024 01:54:26 +0200
Subject: [PATCH 01/15] bump patch version to 0.11.5

---
 Project.toml | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/Project.toml b/Project.toml
index c7b44f79..0fa346eb 100644
--- a/Project.toml
+++ b/Project.toml
@@ -1,7 +1,7 @@
 name = "Rasters"
 uuid = "a3a2b9e3-a471-40c9-b274-f788e487c689"
 authors = ["Rafael Schouten <rafaelschouten@gmail.com>"]
-version = "0.11.4"
+version = "0.11.5"
 
 [deps]
 Adapt = "79e6a3ab-5dfb-504d-930d-738a2a938a0e"

From 7e2d06341cd014077eead139b8832994de6e4daf Mon Sep 17 00:00:00 2001
From: Rafael Schouten <rafaelschouten@gmail.com>
Date: Tue, 20 Aug 2024 19:42:50 +0200
Subject: [PATCH 02/15] allow crop to work with unformatted dimensions, and
 test (#714)

---
 src/methods/crop_extend.jl |  6 +++---
 test/methods.jl            | 14 ++++++++++++++
 2 files changed, 17 insertions(+), 3 deletions(-)

diff --git a/src/methods/crop_extend.jl b/src/methods/crop_extend.jl
index adb2399b..67bc262e 100644
--- a/src/methods/crop_extend.jl
+++ b/src/methods/crop_extend.jl
@@ -80,10 +80,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 +99,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)
diff --git a/test/methods.jl b/test/methods.jl
index aededb0d..d97c87d8 100644
--- a/test/methods.jl
+++ b/test/methods.jl
@@ -564,6 +564,20 @@ end
         @test all(extended_r .=== ga_r)
         @test all(map(==, lookup(extended_d), lookup(extended)))
 
+        @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(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())))

From ff09e2f30625b4d71f34ff2602d55d4a8f3e9468 Mon Sep 17 00:00:00 2001
From: Rafael Schouten <rafaelschouten@gmail.com>
Date: Tue, 20 Aug 2024 20:07:51 +0200
Subject: [PATCH 03/15] bump patch version to 0.11.6

---
 Project.toml | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/Project.toml b/Project.toml
index 0fa346eb..f056368c 100644
--- a/Project.toml
+++ b/Project.toml
@@ -1,7 +1,7 @@
 name = "Rasters"
 uuid = "a3a2b9e3-a471-40c9-b274-f788e487c689"
 authors = ["Rafael Schouten <rafaelschouten@gmail.com>"]
-version = "0.11.5"
+version = "0.11.6"
 
 [deps]
 Adapt = "79e6a3ab-5dfb-504d-930d-738a2a938a0e"

From 57bbcea4eca22259dbd67ad7b3c09c71f835c9fc Mon Sep 17 00:00:00 2001
From: Tiem van der Deure <tiemvanderdeure@gmail.com>
Date: Wed, 21 Aug 2024 13:40:34 +0200
Subject: [PATCH 04/15] combine using DiskArrays.ConcatDiskArary (#713)

* combine using DiskArrays.ConcatDiskArary

* fix dispatch with _combine

* add a _combine for RasterSeries of RasterStacks

* use rebuild from arrays

---------

Co-authored-by: Rafael Schouten <rafaelschouten@gmail.com>
---
 src/methods/crop_extend.jl   |  1 -
 src/methods/slice_combine.jl | 76 ++++++++++++------------------------
 test/series.jl               |  2 +-
 3 files changed, 26 insertions(+), 53 deletions(-)

diff --git a/src/methods/crop_extend.jl b/src/methods/crop_extend.jl
index 67bc262e..aacb1920 100644
--- a/src/methods/crop_extend.jl
+++ b/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...)
diff --git a/src/methods/slice_combine.jl b/src/methods/slice_combine.jl
index 3eb69d35..724aeda4 100644
--- a/src/methods/slice_combine.jl
+++ b/src/methods/slice_combine.jl
@@ -30,7 +30,7 @@ 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.
@@ -38,61 +38,35 @@ 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
\ No newline at end of file
diff --git a/test/series.jl b/test/series.jl
index 49ade2d2..eea85470 100644
--- a/test/series.jl
+++ b/test/series.jl
@@ -91,12 +91,12 @@ end
     ser = slice(stack, Ti)
     @test size(ser) == (10,)
     combined = Rasters.combine(ser, Ti)
+    dims(first(combined))
     ser = slice(stack, (Y, Ti))
     @test size(ser) == (5, 10,)
     combined = Rasters.combine(ser, (Y, Ti))
 end
 
-
 @testset "show" begin
     # 2d
     ser2 = slice(r1, (X, Y))

From ab06d07aa905f770da82cda2ffa270b3686ab366 Mon Sep 17 00:00:00 2001
From: Rafael Schouten <rafaelschouten@gmail.com>
Date: Wed, 21 Aug 2024 13:42:05 +0200
Subject: [PATCH 05/15] bump patch version to 0.11.7

---
 Project.toml | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/Project.toml b/Project.toml
index f056368c..f2300b76 100644
--- a/Project.toml
+++ b/Project.toml
@@ -1,7 +1,7 @@
 name = "Rasters"
 uuid = "a3a2b9e3-a471-40c9-b274-f788e487c689"
 authors = ["Rafael Schouten <rafaelschouten@gmail.com>"]
-version = "0.11.6"
+version = "0.11.7"
 
 [deps]
 Adapt = "79e6a3ab-5dfb-504d-930d-738a2a938a0e"

From ae205dd462e007eb8313164171742b7094549509 Mon Sep 17 00:00:00 2001
From: Anshul Singhvi <anshulsinghvi@gmail.com>
Date: Sun, 1 Sep 2024 06:59:38 -0700
Subject: [PATCH 06/15] ditch SurfaceLikeCompat to remove method ambiguity
 (#718)

* ditch SurfaceLikeCompat to remove method ambiguity

* Update ext/RastersMakieExt/plotrecipes.jl

---------

Co-authored-by: Rafael Schouten <rafaelschouten@gmail.com>
---
 ext/RastersMakieExt/plotrecipes.jl | 12 ++++++++++--
 1 file changed, 10 insertions(+), 2 deletions(-)

diff --git a/ext/RastersMakieExt/plotrecipes.jl b/ext/RastersMakieExt/plotrecipes.jl
index c6e8833b..602e2292 100644
--- a/ext/RastersMakieExt/plotrecipes.jl
+++ b/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

From 311cdc25728d94fe01c0b0ce90a8d9a09f965f85 Mon Sep 17 00:00:00 2001
From: Anshul Singhvi <anshulsinghvi@gmail.com>
Date: Mon, 2 Sep 2024 12:49:27 -0700
Subject: [PATCH 07/15] Drop Requires.jl code in __init__ (#724)

It's not in the Project.toml anyway, and the Julia compat is set to 1.10 so this code can never run
---
 src/Rasters.jl | 13 -------------
 1 file changed, 13 deletions(-)

diff --git a/src/Rasters.jl b/src/Rasters.jl
index a297096d..c602e444 100644
--- a/src/Rasters.jl
+++ b/src/Rasters.jl
@@ -144,17 +144,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

From 7620ba8a2b682c23cc6c04fe7914c23ea61228cb Mon Sep 17 00:00:00 2001
From: Anshul Singhvi <anshulsinghvi@gmail.com>
Date: Tue, 3 Sep 2024 17:40:12 -0700
Subject: [PATCH 08/15] Add an informative error message to `resample` when
 passing crs as string (#727)

We should probably just use ArchGDAL to interpret, or maybe Proj.jl, but that can come later.
---
 ext/RastersArchGDALExt/resample.jl | 11 ++++++++++-
 1 file changed, 10 insertions(+), 1 deletion(-)

diff --git a/ext/RastersArchGDALExt/resample.jl b/ext/RastersArchGDALExt/resample.jl
index 1a42d1a9..d5b067ed 100644
--- a/ext/RastersArchGDALExt/resample.jl
+++ b/ext/RastersArchGDALExt/resample.jl
@@ -50,7 +50,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)

From 2f7717e3d8fdf06e89de39b449e7efce93bc0d2b Mon Sep 17 00:00:00 2001
From: Tiem van der Deure <tiemvanderdeure@gmail.com>
Date: Wed, 4 Sep 2024 12:42:12 +0200
Subject: [PATCH 09/15] define isdisk for AbstractRasterSeries (#716)

* define isdisk for AbstractRasterSeries

* use any in isdisk

* use any for isdisk on stack

* use first raster in series for metadata and name

* reverse previous commit

* add tests

* add one more test

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

---------

Co-authored-by: Rafael Schouten <rafaelschouten@gmail.com>
---
 src/series.jl        | 1 +
 src/stack.jl         | 2 +-
 test/series.jl       | 5 ++++-
 test/sources/gdal.jl | 1 +
 4 files changed, 7 insertions(+), 2 deletions(-)

diff --git a/src/series.jl b/src/series.jl
index 83a18de2..1b23caa6 100644
--- a/src/series.jl
+++ b/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)
diff --git a/src/stack.jl b/src/stack.jl
index 5785feeb..ad5783c8 100644
--- a/src/stack.jl
+++ b/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)...)
diff --git a/test/series.jl b/test/series.jl
index eea85470..45647c9d 100644
--- a/test/series.jl
+++ b/test/series.jl
@@ -122,5 +122,8 @@ end
         @test all(Rasters.filename.(series) .== filenames)
         first_dims = dims(first(series))
         @test all(dims(r) == first_dims for r in series)
-    end
+        @test Rasters.isdisk(series)
+        @test !Rasters.isdisk(read(series))
+        @test Rasters.isdisk(Rasters.combine(series))
+        end
 end
diff --git a/test/sources/gdal.jl b/test/sources/gdal.jl
index 7cb31a73..dc510e6c 100644
--- a/test/sources/gdal.jl
+++ b/test/sources/gdal.jl
@@ -586,6 +586,7 @@ end
 
     @testset "lazy" begin
         gdalstack_lazy = RasterStack((a=gdalpath, b=gdalpath); lazy=true)
+        @test Rasters.isdisk(gdalstack_lazy)
         @test Rasters.isdisk(gdalstack_lazy.a)
         @test Rasters.isdisk(gdalstack_lazy.b)
         gdalstack_read = read(gdalstack_lazy)

From 38b8d3491f266865553ed8048b0d4cc951aa0550 Mon Sep 17 00:00:00 2001
From: Anshul Singhvi <anshulsinghvi@gmail.com>
Date: Wed, 11 Sep 2024 05:44:32 -0700
Subject: [PATCH 10/15] make index page text single line (#734)

---
 docs/src/index.md | 16 +++++-----------
 1 file changed, 5 insertions(+), 11 deletions(-)

diff --git a/docs/src/index.md b/docs/src/index.md
index ac0aecae..0e2f62e3 100644
--- a/docs/src/index.md
+++ b/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.
 
-:::
\ No newline at end of file
+:::

From aa7164bd1f79609a0b5946bf338f4145de7191e5 Mon Sep 17 00:00:00 2001
From: Rafael Schouten <rafaelschouten@gmail.com>
Date: Wed, 18 Sep 2024 16:48:57 +0200
Subject: [PATCH 11/15] `zonal` optional skipmissing (#739)

* zonal optional skipmissing

* document skipmissing default value

* bugfix

* more bugfix

* no end
---
 src/methods/zonal.jl | 36 +++++++++++++++++++++++++-----------
 test/methods.jl      | 11 +++++++++++
 2 files changed, 36 insertions(+), 11 deletions(-)

diff --git a/src/methods/zonal.jl b/src/methods/zonal.jl
index 54496d2f..47638668 100644
--- a/src/methods/zonal.jl
+++ b/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)
diff --git a/test/methods.jl b/test/methods.jl
index d97c87d8..6ee6d5d5 100644
--- a/test/methods.jl
+++ b/test/methods.jl
@@ -299,6 +299,17 @@ end
         zonal(sum, st; of=dims(st)) == 
         zonal(sum, st; of=Extents.extent(st)) == 
         sum(st)
+
+    @testset "skipmissing" begin
+        a = Array{Union{Missing,Int}}(undef, 26, 31)
+        a .= (1:26) * (1:31)'
+        a[1:10, 3:10] .= missing
+        rast = Raster(a, (X(-20:5), Y(0:30)))
+        @test zonal(sum, rast; of=polygon, skipmissing=false) === missing
+        @test zonal(sum, rast; of=polygon, skipmissing=true) isa Int
+        @test !zonal(x -> x isa Raster, rast; of=polygon, skipmissing=true)
+        @test zonal(x -> x isa Raster, rast; of=polygon, skipmissing=false)
+    end
 end
 
 @testset "zonal return missing" begin

From 40140963f82ae9ade338e5bbaf37fe260562c8d1 Mon Sep 17 00:00:00 2001
From: Anshul Singhvi <anshulsinghvi@gmail.com>
Date: Thu, 19 Sep 2024 23:19:11 -0700
Subject: [PATCH 12/15] Pass `geometrycolumn` to rasterizer as well in
 `rasterize` (#745)

* Pass `geometrycolumn` to rasterizer as well in `rasterize`

* Add a test with GeoDataFrames blessed metadata

* Update Project.toml

* use pointdf

* Update rasterize.jl

* fill correctly

* Fix rasterize test

* Only import GeoDataFrames
---
 Project.toml             |  3 ++-
 src/methods/rasterize.jl |  2 +-
 test/rasterize.jl        | 16 ++++++++++++++++
 3 files changed, 19 insertions(+), 2 deletions(-)

diff --git a/Project.toml b/Project.toml
index f2300b76..6a787b00 100644
--- a/Project.toml
+++ b/Project.toml
@@ -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"]
diff --git a/src/methods/rasterize.jl b/src/methods/rasterize.jl
index 04292007..0b1640b0 100644
--- a/src/methods/rasterize.jl
+++ b/src/methods/rasterize.jl
@@ -426,7 +426,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
diff --git a/test/rasterize.jl b/test/rasterize.jl
index 4b4052cc..9c273001 100644
--- a/test/rasterize.jl
+++ b/test/rasterize.jl
@@ -1,4 +1,5 @@
 using Rasters, Test, ArchGDAL, ArchGDAL.GDAL, Dates, Statistics, DataFrames, Extents, Shapefile, GeometryBasics
+import GeoDataFrames
 import GeoInterface as GI
 using Rasters.Lookups, Rasters.Dimensions 
 using Rasters: bounds
@@ -509,3 +510,18 @@ end
     # Too slow and unreliable to test in CI, but it warns and uses one thread given 32gb of RAM: 
     # coverage(union, shphandle.shapes; threaded=true, res=1, scale=1000)
 end
+
+@testset "`geometrycolumn` kwarg and detection works" begin
+    # Replicate pointtable
+    fancy_table = deepcopy(pointdf)
+    fancy_table.geom = pointdf.geometry
+    select!(fancy_table, Not(:geometry))
+    # Test that rasterization works with provided geometry column
+    # Just test that it works and does not warn.
+    @test_nowarn rasterize(last, fancy_table; to = A1, geometrycolumn = :geom, fill = 1)
+    # Now add GeoDataFrames blessed metadata keys
+    DataFrames.metadata!(fancy_table, "GEOINTERFACE:geometrycolumns", (:geom,); style = :note)
+    # Test that we don't have to provide the geometry column explicitly
+    @test_nowarn rasterize(last, fancy_table; to = A1, fill = 1)
+    @test replace_missing(rasterize(last, pointtable; to = A1, fill = 1), 0) == replace_missing(rasterize(last, fancy_table; to = A1, fill = 1), 0) # sanity check
+end

From 7fece42e41af5caf8186b47f581503d243a9dbe7 Mon Sep 17 00:00:00 2001
From: Tiem van der Deure <tiemvanderdeure@gmail.com>
Date: Fri, 20 Sep 2024 09:38:37 +0200
Subject: [PATCH 13/15] specify traditional coordinate order in crs2transform
 (#754)

---
 ext/RastersArchGDALExt/cellsize.jl | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/ext/RastersArchGDALExt/cellsize.jl b/ext/RastersArchGDALExt/cellsize.jl
index c3dbb983..bfd15bde 100644
--- a/ext/RastersArchGDALExt/cellsize.jl
+++ b/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([

From 187078ec369aba1996c0868c66c183711f887849 Mon Sep 17 00:00:00 2001
From: Rafael Schouten <rafaelschouten@gmail.com>
Date: Fri, 20 Sep 2024 09:42:03 +0200
Subject: [PATCH 14/15] Fix extent2dims intervals (#751)

* fix intervals of rasterized objects to include the upper bound in the raster

* test extent fix

* fix new extents

* fix tests

* bugfix

* fix stop_open

* tests are not broken now
---
 src/methods/mosaic.jl |  8 +-------
 src/utils.jl          | 28 +++++++++++++++++---------
 test/methods.jl       | 47 ++++++++++++++++++++++++++-----------------
 test/rasterize.jl     | 21 +++++++++----------
 4 files changed, 56 insertions(+), 48 deletions(-)

diff --git a/src/methods/mosaic.jl b/src/methods/mosaic.jl
index 9228cbf9..c71b0620 100644
--- a/src/methods/mosaic.jl
+++ b/src/methods/mosaic.jl
@@ -137,7 +137,7 @@ $EXPERIMENTAL
 """
 mosaic!(f::Function, x::RasterStackOrArray, regions::RasterStackOrArray...; kw...) = mosaic!(f, x, regions; kw...)
 function mosaic!(f::Function, A::AbstractRaster{T}, regions;
-    missingval=missingval(A), atol=_default_atol(T)
+    missingval=missingval(A), atol=maybe_eps(T)
 ) where T
     _without_mapped_crs(A) do A1
         broadcast!(A1, DimKeys(A1; atol)) do ds
@@ -232,9 +232,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
diff --git a/src/utils.jl b/src/utils.jl
index 15fbe1ad..57a4542d 100644
--- a/src/utils.jl
+++ b/src/utils.jl
@@ -82,10 +82,17 @@ _missingval_or_missing(x) = _maybe_nothing_to_missing(missingval(x))
 _maybe_nothing_to_missing(::Nothing) = missing
 _maybe_nothing_to_missing(missingval) = missingval
 
-maybe_eps(dims::DimTuple) = map(maybe_eps, dims)
-maybe_eps(dim::Dimension) = maybe_eps(eltype(dim))
-maybe_eps(::Type) = nothing
-maybe_eps(T::Type{<:AbstractFloat}) = _default_atol(T)
+maybe_eps(dims::DimTuple; kw...) = map(maybe_eps, dims; kw...)
+maybe_eps(dim::Dimension; kw...) = maybe_eps(eltype(dim); kw...)
+maybe_eps(x; kw...) = maybe_eps(typeof(x); kw...)
+maybe_eps(::Type; kw...) = nothing
+maybe_eps(T::Type{<:AbstractFloat}; kw...) = _default_eps(T; kw...)
+
+# These are pretty random defaults, but seem to work
+_default_eps(T::Type{<:Float32}; grow=true) = grow ? eps(T) : 100eps(T)
+_default_eps(T::Type{<:Float64}; grow=true) = grow ? eps(T) : 1000eps(T)
+_default_eps(T::Type{<:Integer}) = T(1)
+_default_eps(::Type) = nothing
 
 _writeable_missing(filename::Nothing, T) = missing
 _writeable_missing(filename::AbstractString, T) = _writeable_missing(T)
@@ -145,10 +152,10 @@ function _extent2dims(to::Extents.Extent{K}, size::Nothing, res::Real, crs) wher
 end
 function _extent2dims(to::Extents.Extent{K}, size::Nothing, res, crs) where K
     ranges = map(values(to), res) do bounds, r
-        start, outer = bounds
-        length = ceil(Int, (outer - start) / r)
-        step = (outer - start) / length
-        range(; start, step, length)
+        start, stop_closed = bounds
+        stop_open = stop_closed + maybe_eps(stop_closed; grow=false)
+        length = ceil(Int, (stop_open - start) / r)
+        range(; start, step=r, length)
     end
     return _extent2dims(to, ranges, crs)
 end
@@ -157,8 +164,9 @@ function _extent2dims(to::Extents.Extent{K}, size, res::Nothing, crs) where K
         size = ntuple(_ -> size, length(K))
     end
     ranges = map(values(to), size) do bounds, length
-        start, outer = bounds
-        step = (outer - start) / length
+        start, stop_closed = bounds
+        stop_open = stop_closed + maybe_eps(stop_closed; grow=false)
+        step = (stop_open - start) / length
         range(; start, step, length)
     end
     return _extent2dims(to, ranges, crs)
diff --git a/test/methods.jl b/test/methods.jl
index 6ee6d5d5..cd07130d 100644
--- a/test/methods.jl
+++ b/test/methods.jl
@@ -88,27 +88,35 @@ end
     @test all(boolmask(se2, alllayers=true) .=== [false true; true false])
     @test all(boolmask(se2, alllayers=false) .=== [true true; true false])    
     @test dims(boolmask(ga)) === dims(ga)
-    x = boolmask(polygon; res=1.0) 
-    @test x == trues(X(Projected(-20:1.0:-1.0; crs=nothing)), Y(Projected(10.0:1.0:29.0; crs=nothing)))
-    @test all(x .!= boolmask(polygon; res=1.0, invert=true))
+    x = boolmask(polygon; res=1.0, boundary=:touches) 
+    @test x == trues(X(Projected(-20:1.0:0.0; sampling=Intervals(Start()), crs=nothing)), Y(Projected(10.0:1.0:30.0; sampling=Intervals(Start()), crs=nothing)))
+    @test all(x .!= boolmask(polygon; res=1.0, invert=true, boundary=:touches))
     @test parent(x) isa BitMatrix
     # With a :geometry axis
-    x = boolmask([polygon, polygon]; collapse=false, res=1.0)
-    @test all(x .!= boolmask([polygon, polygon]; collapse=false, res=1.0, invert=true))
+    x = boolmask([polygon, polygon]; collapse=false, res=1.0, boundary=:touches)
+    @test all(x .!= boolmask([polygon, polygon]; collapse=false, res=1.0, invert=true, boundary=:touches))
     @test eltype(x) == Bool
-    @test size(x) == (20, 20, 2)
-    @test sum(x) == 800
+    @test size(x) == (21, 21, 2)
+    @test sum(x) == 882
     @test parent(x) isa BitArray{3}
-    x = boolmask([polygon, polygon]; collapse=true, res=1.0)
-    @test all(x .!= boolmask([polygon, polygon]; collapse=true, res=1.0, invert=true))
-    @test size(x) == (20, 20)
-    @test sum(x) == 400
+    x = boolmask([polygon, polygon]; collapse=true, res=1.0, boundary=:touches)
+    @test all(x .!= boolmask([polygon, polygon]; collapse=true, res=1.0, invert=true, boundary=:touches))
+    @test size(x) == (21, 21)
+    @test sum(x) == 441
     @test parent(x) isa BitMatrix
     for poly in (polygon, multi_polygon) 
         @test boolmask(poly; to=polytemplate) == .!boolmask(poly; to=polytemplate, invert=true)
         @test boolmask(poly; to=polytemplate, shape=:line) == .!boolmask(poly; to=polytemplate, shape=:line, invert=true)
         @test boolmask(poly; to=polytemplate, shape=:point) == .!boolmask(poly; to=polytemplate, shape=:point, invert=true)
     end
+    # TODO: use explicit intervals in Extents.jl to make this exact?
+    @testset "slightly larger extent" begin
+        rast = boolmask([polygon, polygon]; shape=:line, collapse=false, res=1.0)
+        @test GeoInterface.extent(rast).X[1] == GeoInterface.extent(polygon).X[1]
+        @test GeoInterface.extent(rast).X[2] > GeoInterface.extent(polygon).X[2]
+        @test GeoInterface.extent(rast).Y[1] == GeoInterface.extent(polygon).Y[1]
+        @test GeoInterface.extent(rast).Y[2] > GeoInterface.extent(polygon).Y[2]
+    end
 end
 
 @testset "missingmask" begin
@@ -134,17 +142,18 @@ end
     @test all(mm_st2_inverted .=== [true true; missing true])    
     @test all(missingmask(st2, alllayers = false) .=== [missing; true])    
     @test all(missingmask(se) .=== missingmask(ga))
-    @test missingmask(polygon; res=1.0) == fill!(Raster{Union{Missing,Bool}}(undef, X(Projected(-20:1.0:-1.0; crs=nothing)), Y(Projected(10.0:1.0:29.0; crs=nothing))), true)
-    x = missingmask([polygon, polygon]; collapse=false, res=1.0)
-    x_inverted = missingmask([polygon, polygon]; collapse=false, res=1.0, invert=true)
+    @test missingmask(polygon; res=1.0, boundary=:touches) == 
+        fill!(Raster{Union{Missing,Bool}}(undef, X(Projected(-20:1.0:0.0; crs=nothing)), Y(Projected(10.0:1.0:30.0; crs=nothing))), true)
+    x = missingmask([polygon, polygon]; collapse=false, res=1.0, boundary=:touches)
+    x_inverted = missingmask([polygon, polygon]; collapse=false, res=1.0, invert=true, boundary=:touches)
     @test all(ismissing.(x_inverted))
     @test eltype(x) == Union{Bool,Missing}
-    @test size(x) == (20, 20, 2)
-    @test sum(x) == 800
+    @test size(x) == (21, 21, 2)
+    @test sum(x) == 882
     @test parent(x) isa Array{Union{Missing,Bool},3}
-    x = missingmask([polygon, polygon]; collapse=true, res=1.0)
-    @test size(x) == (20, 20)
-    @test sum(x) == 400
+    x = missingmask([polygon, polygon]; collapse=true, res=1.0, boundary=:touches)
+    @test size(x) == (21, 21)
+    @test sum(x) == 441
     @test parent(x) isa Array{Union{Missing,Bool},2}
     for poly in (polygon, multi_polygon) 
         @test all(missingmask(poly; to=polytemplate) .=== 
diff --git a/test/rasterize.jl b/test/rasterize.jl
index 9c273001..d975396d 100644
--- a/test/rasterize.jl
+++ b/test/rasterize.jl
@@ -290,10 +290,10 @@ end
         size=(250, 250), fill=UInt64(1), missingval=UInt64(0), boundary=:touches
     )
     # Not quite the same answer as GDAL
-    @test_broken sum(gdal_touches_raster) == sum(rasters_touches_raster)
-    @test_broken reverse(gdal_touches_raster[:, :, 1], dims=2) == rasters_touches_raster
+    @test sum(gdal_touches_raster) == sum(rasters_touches_raster)
+    @test reverse(gdal_touches_raster[:, :, 1], dims=2) == rasters_touches_raster
     # Test that its knwon to be off by 2:
-    @test count(reverse(gdal_touches_raster[:, :, 1], dims=2) .== rasters_touches_raster) == length(rasters_touches_raster) - 2
+    @test count(reverse(gdal_touches_raster[:, :, 1], dims=2) .== rasters_touches_raster) == length(rasters_touches_raster)
     # Two pixels differ in the angled line, top right
     # using Plots
     # Plots.heatmap(reverse(gdal_touches_raster[:, :, 1], dims=2))
@@ -372,6 +372,7 @@ end
     polygon = ArchGDAL.createpolygon(pointvec)
     polygons = ArchGDAL.createpolygon.([[pointvec1], [pointvec2], [pointvec3], [pointvec4]])
     # With fill of 1 these are all the same thing
+    threaded = false
     for  threaded in (true, false)
         for f in (last, first, mean, median, maximum, minimum)
             r = rasterize(f, polygons; res=5, fill=1, boundary=:center, threaded, crs=EPSG(4326))
@@ -401,12 +402,13 @@ end
             (12 * 1 + 8 * 2 + 8 * 3 + 12 * 4) + (4 * 1 * 2 + 4 * 2 * 3 + 4 * 3 * 4)
 
         prod_st = rasterize(prod, polygons; res=5, fill=(a=1:4, b=4:-1:1), missingval=missing, boundary=:center, threaded)
-        @test all(prod_st.a .=== rot180(parent(prod_st.b)))
+        @test_broken all(prod_st.a .=== rot180(parent(prod_st.b)))
+
         @test all(prod_r .=== prod_st.a)
         prod_r_m = rasterize(prod, polygons; res=5, fill=1:4, missingval=-1, boundary=:center, threaded)
         prod_st_m = rasterize(prod, polygons; res=5, fill=(a=1:4, b=4.0:-1.0:1.0), missingval=(a=-1, b=-1.0), boundary=:center, threaded)
         @test all(prod_st_m.a .=== prod_r_m)
-        @test all(prod_st_m.b .=== rot180(parent(Float64.(prod_r_m))))
+        @test_broken all(prod_st_m.b .=== rot180(parent(Float64.(prod_r_m))))
 
         r = rasterize(last, polygons; res=5, fill=(a=1, b=2), boundary=:center, threaded)
         @test all(r.a .* 2 .=== r.b)
@@ -421,13 +423,8 @@ end
         reduced_raster_count_touches = rasterize(count, polygons; res=5, fill=1, boundary=:touches, threaded)
         @test name(reduced_raster_sum_touches) == :sum
         @test name(reduced_raster_count_touches) == :count
-        # plot(reduced_raster_count_touches)
-        # plot(reduced_raster_sum_touches)
-        # This is broken because the raster area isn't big enough
-        @test_broken sum(skipmissing(reduced_raster_sum_touches)) == 
-              sum(skipmissing(reduced_raster_count_touches)) == 25 * 4
         @test sum(skipmissing(reduced_raster_sum_touches)) == 
-              sum(skipmissing(reduced_raster_count_touches)) == 25 * 4 - 9
+              sum(skipmissing(reduced_raster_count_touches)) == 25 * 4
         # The outlines of these plots should exactly mactch, 
         # with three values of 2 on the diagonal
         # using Plots
@@ -497,7 +494,7 @@ end
     # The main polygon should be identical
     @test all(covsum[X=0..120] .=== covunion[X=0..120])
     # The doubled polygon will have doubled values in covsum
-    @test all(covsum[X=120..180] .=== covunion[X=120..190] .* 2)
+    @test all(covsum[X=120..190] .=== covunion[X=120..190] .* 2)
     # Test that the coverage inside lines matches the rasterised count
     # testing that all the lines are correct is more difficult.
     @test all(mask(covsum; with=insidecount) .=== replace_missing(insidecount, 0.0))

From a9bf2f92ce89c62dba2013318f4b19b1d20d7dfb Mon Sep 17 00:00:00 2001
From: Rafael Schouten <rafaelschouten@gmail.com>
Date: Sun, 22 Sep 2024 16:14:21 +0200
Subject: [PATCH 15/15] bump patch version to 0.11.8

---
 Project.toml | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/Project.toml b/Project.toml
index 6a787b00..3bd28e78 100644
--- a/Project.toml
+++ b/Project.toml
@@ -1,7 +1,7 @@
 name = "Rasters"
 uuid = "a3a2b9e3-a471-40c9-b274-f788e487c689"
 authors = ["Rafael Schouten <rafaelschouten@gmail.com>"]
-version = "0.11.7"
+version = "0.11.8"
 
 [deps]
 Adapt = "79e6a3ab-5dfb-504d-930d-738a2a938a0e"