diff --git a/.github/ISSUE_TEMPLATE/bug-report.md b/.github/ISSUE_TEMPLATE/bug_report.md
similarity index 59%
rename from .github/ISSUE_TEMPLATE/bug-report.md
rename to .github/ISSUE_TEMPLATE/bug_report.md
index 9551724b0..dede49a57 100644
--- a/.github/ISSUE_TEMPLATE/bug-report.md
+++ b/.github/ISSUE_TEMPLATE/bug_report.md
@@ -1,15 +1,19 @@
 ---
-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.
-
+name: Bug report
+about: Describe your issue here, including a minimum working example all files.
 title: ''
-labels: ''
+labels: bug
 assignees: ''
 
 ---
 
-Your problem. You don't have to write anything here for a bug, the code and error are what is important.
+########################
+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`. 
+########################
 
 ```julia
 All the julia code needed to make the error goes here including download or generation of all data used
diff --git a/.github/workflows/ci.yml b/.github/workflows/ci.yml
index 6c8b91be8..644eca742 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/Project.toml b/Project.toml
index 3bd28e78f..7c6985c79 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.8"
+version = "0.12.1"
 
 [deps]
 Adapt = "79e6a3ab-5dfb-504d-930d-738a2a938a0e"
@@ -16,6 +16,8 @@ FillArrays = "1a297f60-69ca-5386-bcde-b61e274b549b"
 Flatten = "4c728ea3-d9ee-5c9a-9642-b6f7d7dc04fa"
 GeoFormatTypes = "68eda718-8dee-11e9-39e7-89f7f65f511f"
 GeoInterface = "cf35fbd7-0cd7-5166-be24-54bfbe79505f"
+GeometryOpsCore = "05efe853-fabf-41c8-927e-7063c8b9f013"
+LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 Missings = "e1d29d7a-bbdc-5cf2-9ac0-f12de2c33e28"
 Mmap = "a63ad114-7e13-5084-954f-fe012c677804"
 OffsetArrays = "6fe1bfb0-de20-5000-8ca7-80f57d26f881"
@@ -23,6 +25,7 @@ ProgressMeter = "92933f4c-e287-5a05-a399-4b506db050ca"
 RecipesBase = "3cdcf5f2-1ef4-517c-9805-6587b60abb01"
 Reexport = "189a3867-3050-52da-a836-e630ba90ab69"
 Setfield = "efcf1570-3423-57d1-acb7-fd33fddbac46"
+Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"
 
 [weakdeps]
 ArchGDAL = "c9ce4bd3-c3d5-55b8-8973-c0e20141b8c3"
@@ -30,7 +33,9 @@ CoordinateTransformations = "150eb455-5306-5404-9cee-2592286d6298"
 GRIBDatasets = "82be9cdb-ee19-4151-bdb3-b400788d9abc"
 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]
@@ -39,7 +44,9 @@ RastersCoordinateTransformationsExt = "CoordinateTransformations"
 RastersGRIBDatasetsExt = "GRIBDatasets"
 RastersMakieExt = "Makie"
 RastersNCDatasetsExt = "NCDatasets"
+RastersProjExt = "Proj"
 RastersRasterDataSourcesExt = "RasterDataSources"
+RastersStatsBaseExt = "StatsBase"
 RastersZarrDatasetsExt = "ZarrDatasets"
 
 [compat]
@@ -52,26 +59,32 @@ CommonDataModel = "0.2.3, 0.3"
 ConstructionBase = "1"
 CoordinateTransformations = "0.6.2"
 DataFrames = "1"
-DimensionalData = "0.27.3"
+DimensionalData = "0.29.4"
 DiskArrays = "0.3, 0.4"
 Extents = "0.1"
 FillArrays = "0.12, 0.13, 1"
 Flatten = "0.4"
 GRIBDatasets = "0.2, 0.3"
+GeoDataFrames = "0.3"
 GeoFormatTypes = "0.4"
-GeoInterface = "1"
-Makie = "0.19, 0.20, 0.21"
+GeoInterface = "1.0"
+GeometryBasics = "0.4"
+GeometryOpsCore = "0.1.1"
+Makie = "0.20, 0.21"
 Missings = "0.4, 1"
 NCDatasets = "0.13, 0.14"
 OffsetArrays = "1"
+Plots = "1"
 ProgressMeter = "1"
-RasterDataSources = "0.5.7, 0.6"
+Proj = "1.7.2"
+RasterDataSources = "0.7"
 RecipesBase = "0.7, 0.8, 1.0"
 Reexport = "0.2, 1.0"
 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"
@@ -82,17 +95,19 @@ 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"
+GeoDataFrames = "62cb38b5-d8d2-4862-a48e-6a340996859f"
 GeometryBasics = "5c1252a2-5f33-56bf-86c9-59e7332b4326"
 Makie = "ee78f7c6-11fb-53f2-987a-cfe4a2b5a57a"
 NCDatasets = "85f8d34a-cbdd-5861-8df4-14fed0d494ab"
 Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80"
+Proj = "c94c279d-25a6-4763-9509-64d165bea63e"
 RasterDataSources = "3cb90ccd-e1b6-4867-9617-4276c8b2ca36"
 SafeTestsets = "1bc83da4-3b8d-516f-aca4-4fe02f6d838f"
 Shapefile = "8e980c4a-a4fe-5da2-b3a7-4b4b0353a2f4"
-Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"
+StableRNGs = "860ef19b-820b-49d6-a774-d7a799459cd3"
+StatsBase = "2913bbd2-ae8a-5f71-8c99-4fb6c76f3a91"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 
 [targets]
-test = ["Aqua", "ArchGDAL", "CFTime", "CoordinateTransformations", "DataFrames", "GeoDataFrames", "GeometryBasics", "GRIBDatasets", "NCDatasets", "Plots", "RasterDataSources", "SafeTestsets", "Shapefile", "Statistics", "Test", "ZarrDatasets"]
+test = ["Aqua", "ArchGDAL", "CFTime", "CoordinateTransformations", "DataFrames", "GeoDataFrames", "GeometryBasics", "GRIBDatasets", "NCDatasets", "Plots", "Proj", "RasterDataSources", "SafeTestsets", "Shapefile", "StableRNGs", "StatsBase", "Test", "ZarrDatasets"]
diff --git a/README.md b/README.md
index f226aa91c..275826f39 100644
--- a/README.md
+++ b/README.md
@@ -16,36 +16,6 @@ manipulating rasterized spatial data.
 These currently include raster arrays like GeoTIFF and NetCDF, R grd files, 
 multi-layered stacks, and multi-file series of arrays and stacks. 
 
-
-## :warning: Packages extensions and Rasters 0.8 and onwards 
-
-On Julia 1.9 we can put additional packages in extensions, so the code only loads when
-you load a specific package. Rasters.jl was always intended to work like this,
-and its finally possible. This reduced package `using` time from many seconds to well under a second.
-
-But, it means you have to manually load packages you need for each backend or additional
-functionality.
-
-For example, to use the GDAL backend, and download files, you now need to do:
-
-```julia
-using Rasters, ArchGDAL, RasterDataSources
-```
-
-where previously it was just `using Rasters`.
-
-Sources and packages needed:
-- `:gdal`: `using ArchGDAL`
-- `:netcdf`: `using NCDatasets`
-- `:grd`: built-in.
-- `:smap`: `using HDF5`
-- `:grib`: `using GRIBDatasets`.
-
-Other functionality in extensions:
-- Raster data downloads, like `Worldclim{Climate}`: `using RasterDataSources`
-- Makie plots: `using Makie`
-- Coordinate transformations for gdal rasters: `using CoordinateTransformations`
-
 # Quick start
 Install the package by typing:
 
@@ -54,6 +24,7 @@ Install the package by typing:
 add Rasters
 ```
 
+Then to use it:
 ```julia
 using Rasters
 ```
@@ -87,6 +58,29 @@ values: [:, :, 1]
  30   0.334152   0.136551    0.183555    0.941133   0.450484    0.461862
 [and 12 more slices...]
 ```
+## Packages that work with Rasters
+
+Rasters reduces its dependencies to keep the `using` time low.
+But, it means you have to manually load packages you need for each 
+backend or additional functionality.
+
+For example, to use the GDAL backend, and download RasterDataSources files, you now need to do:
+
+```julia
+using Rasters, ArchGDAL, RasterDataSources
+```
+
+Sources and packages needed:
+- `:gdal`: `using ArchGDAL`
+- `:netcdf`: `using NCDatasets`
+- `:grd`: built-in.
+- `:grib`: `using GRIBDatasets`.
+- `:zarr`: `using ZarrDatasets`.
+
+Other functionality in extensions:
+- Raster data downloads, like `Worldclim{Climate}`: `using RasterDataSources`
+- Makie plots: `using GLMakie` (opengl interactive) or `using CairoMakie` (print) etc.
+- Coordinate transformations for gdal rasters: `using CoordinateTransformations`
 
 ## Getting the `lookup` array from dimensions
 
@@ -215,6 +209,13 @@ be used with identical syntax.
   projection, the conversion is handled automatically. This means lat/lon
   `EPSG(4326)` can be used seamlessly if you need that.
 
+# Performance
+
+Rasters should be the fastest tool available for most tasks. If you find 
+something is faster in another package, it's likely a bug - so make an issue!
+
+![image](https://github.com/user-attachments/assets/1c6c56ac-4c5a-4096-984d-15bf2783682c)
+
 
 # Bugs, errors and making issues for Rasters.jl
 
@@ -226,8 +227,11 @@ Because there are so many raster file types and variations of them, most of the
 
 To make an issue we can fix quickly (or at all) there are three key steps:
 
-1. Include the file in an accessible place on web *without authentication* or any other work on our part, so we can just get it and find your bug. You can put it on a file hosting platform (e.g. google drive, drop box, whatever you use) and share the url.
-2. Add a minimum working example to the issue template that first downloads the file, then runs the function that triggers the bug.
-3. Paste the complete stack trace of the error it produces, right to the bottom, into the issue template. Then we can be sure we reproduced the same problem.
+1. Use a RasterDataSources.jl file if you can, so there are no download hassles.
+   Otherwise store a file in an accessible place on web *without authentication* and preferably where you
+   can use `dowload` directly, so we just run the script can spend our time finding your bug.
+2. Add a minimum working example to the issue template that first downloads the file with `download`, then runs the function that triggers the bug.
+3. Paste the complete stack trace of the error it produces, right to the bottom, into the issue template.
+   Then we can be sure we have reproduced the same problem.
 
 Good issues are really appreciated, but they do take just a little extra effort with Rasters.jl because of this need for files.
diff --git a/docs/Project.toml b/docs/Project.toml
index 13e76347b..95ad821d4 100644
--- a/docs/Project.toml
+++ b/docs/Project.toml
@@ -12,9 +12,16 @@ GBIF2 = "dedd4f52-e074-43bf-924d-d6bce14ad628"
 GeoInterface = "cf35fbd7-0cd7-5166-be24-54bfbe79505f"
 HDF5 = "f67ccb44-e63f-5c2f-98bd-6dc0ccc4ba2f"
 Literate = "98b081ad-f1c9-55d3-8b20-4c87d4299306"
+MLJBase = "a7f614a8-145f-11e9-1d2a-a57a1082229d"
+MLJGLMInterface = "caf8df21-4939-456d-ac9c-5fefbfb04c0c"
 Makie = "ee78f7c6-11fb-53f2-987a-cfe4a2b5a57a"
+Maxnet = "81f79f80-22f2-4e41-ab86-00c11cf0f26f"
 NCDatasets = "85f8d34a-cbdd-5861-8df4-14fed0d494ab"
 Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80"
 RasterDataSources = "3cb90ccd-e1b6-4867-9617-4276c8b2ca36"
 Rasters = "a3a2b9e3-a471-40c9-b274-f788e487c689"
 Shapefile = "8e980c4a-a4fe-5da2-b3a7-4b4b0353a2f4"
+SpeciesDistributionModels = "3ef73bbf-0321-4d3b-9a2e-5fbebc8e35da"
+
+[sources]
+SpeciesDistributionModels = {url = "https://github.com/tiemvanderdeure/SpeciesDistributionModels.jl/"}
\ No newline at end of file
diff --git a/docs/src/.vitepress/config.mts b/docs/src/.vitepress/config.mts
index c94117066..bac1aef86 100644
--- a/docs/src/.vitepress/config.mts
+++ b/docs/src/.vitepress/config.mts
@@ -29,7 +29,7 @@ export default defineConfig({
   themeConfig: {
     outline: 'deep',
     // https://vitepress.dev/reference/default-theme-config
-    logo: 'REPLACE_ME_DOCUMENTER_VITEPRESS',
+    logo: '/logo.png',
     search: {
       provider: 'local',
       options: {
@@ -44,27 +44,30 @@ export default defineConfig({
           { text: 'Overview', link: '/methods' },
           { text: 'Array Operations', link: '/array_operations' },
         ]
-       },
-       { text: 'Data Sources',
-       items: [
-         { text: 'Overview', link: '/data_sources' },
-         { text: 'GBIF', link: '/gbif_wflow' }
-        ]
       },
+      { text: 'Data Sources', link: '/data_sources' },
       { text: 'Plots',
-      items: [
-        { text: 'Plots.jl', link: '/plotting' },
-        { text: 'Makie.jl', link: '/plot_makie' },
-      ]
-     },
-     { text: 'Ecosystem',
-      items: [
-        { text: 'DimensionalData.jl', link: 'https://rafaqz.github.io/DimensionalData.jl/dev/' },
-        { text: 'NCDatasets.jl', link: 'https://alexander-barth.github.io/NCDatasets.jl/stable/' },
-        { text: 'ArchGDAL.jl', link: 'https://yeesian.com/ArchGDAL.jl/stable/' },
-        { text: 'HDF5.jl', link: 'https://juliaio.github.io/HDF5.jl/stable/' },
-       ]
-     },
+        items: [
+          { text: 'Plots.jl', link: '/plotting' },
+          { text: 'Makie.jl', link: '/plot_makie' },
+        ]
+      },  
+      { text: 'Examples',
+        items: [
+          { text: 'Species Distribution Modelling', link: '/gbif_wflow' },
+        ]
+      },
+      { text: 'Ecosystem',
+        items: [
+          { text: 'DimensionalData.jl', link: 'https://rafaqz.github.io/DimensionalData.jl' },
+          { text: 'DiskArrays.jl', link: 'https://github.com/JuliaIO/DiskArrays.jl' },
+          { text: 'GeoInterface.jl', link: 'https://github.com/JuliaGeo/GeoInterface.jl' },
+          { text: 'NCDatasets.jl', link: 'https://alexander-barth.github.io/NCDatasets.jl/stable/' },
+          { text: 'ArchGDAL.jl', link: 'https://github.com/yeesian/ArchGDAL.jl' },
+          { text: 'GRIBDatasets.jl', link: 'https://github.com/JuliaGeo/GRIBDatasets.jl' },
+          { text: 'ZarrDatasets.jl', link: 'https://github.com/JuliaGeo/ZarrDatasets.jl' },
+        ]
+      },
       { text: 'API', link: '/api' }
     ],
 
@@ -75,20 +78,19 @@ export default defineConfig({
           { text: 'Overview', link: '/methods' },
           { text: 'Array Operations', link: '/array_operations' },
         ]
-       },
-       { text: 'Data Sources',
-       items: [
-         { text: 'Overview', link: '/data_sources' },
-         { text: 'GBIF', link: '/gbif_wflow' }
-        ]
       },
+      { text: 'Data Sources', link: '/data_sources' },
       { text: 'Plots',
-      items: [
-        { text: 'Plots.jl', link: '/plotting' },
-        { text: 'Makie.jl', link: '/plot_makie' },
-      ]
-     },
-
+        items: [
+          { text: 'Plots.jl', link: '/plotting' },
+          { text: 'Makie.jl', link: '/plot_makie' },
+        ]
+      },
+      { text: 'Examples',
+        items: [
+          { text: 'Species Distribution Modelling', link: '/gbif_wflow' },
+        ]
+      },
       { text: 'API', link: '/api' }
     ],
     editLink: {
@@ -102,4 +104,4 @@ export default defineConfig({
       copyright: `© Copyright ${new Date().getUTCFullYear()}. Released under the MIT License.`
     }
   }
-})
\ No newline at end of file
+})
diff --git a/docs/src/.vitepress/theme/style.css b/docs/src/.vitepress/theme/style.css
index 9816339b8..f89d0a690 100644
--- a/docs/src/.vitepress/theme/style.css
+++ b/docs/src/.vitepress/theme/style.css
@@ -279,4 +279,17 @@ kbd {
   .VPDoc {
     padding-left: 25px !important;  
   }
+}
+/* Component: Docstring Custom Block */
+
+.jldocstring.custom-block {
+  border: 1px solid var(--vp-c-gray-2);
+  color: var(--vp-c-text-1)
+}
+
+.jldocstring.custom-block summary {
+  font-weight: 700;
+  cursor: pointer;
+  user-select: none;
+  margin: 0 0 8px;
 }
\ No newline at end of file
diff --git a/docs/src/data_sources.md b/docs/src/data_sources.md
index 0b93ba996..2604d9ca0 100644
--- a/docs/src/data_sources.md
+++ b/docs/src/data_sources.md
@@ -3,13 +3,17 @@
 Rasters.jl uses a number of backends to load raster data. `Raster`, `RasterStack`
 and `RasterSeries` will detect which backend to use for you, automatically.
 
-## GRD
+## GDAL
 
-R GRD files can be loaded natively, using Julias `MMap` - which means they are very fast, but are not compressed. They are always 3 dimensional, and have `Y`, `X` and [`Band`](@ref) dimensions.
+All files GDAL can access, such as `.tiff` and `.asc` files, can be loaded,
+using [ArchGDAL.jl](https://github.com/yeesian/ArchGDAL.jl). These are
+generally best loaded as `Raster("filename.tif")`, but can be loaded as
+`RasterStack("filename.tif"; layersfrom=Band)`, taking layers from the `Band`
+dimension, which is also the default.
 
 ## NetCDF
 
-NetCDF `.nc` files are loaded using
+NetCDF `.nc` and some HDF5 `.h5` files cab be loaded using
 [NCDatasets.jl](https://github.com/Alexander-Barth/NCDatasets.jl). Layers from
 files can be loaded as `Raster("filename.nc"; name=:layername)`. Without `name`
 the first layer is used. `RasterStack("filename.nc")` will use all netcdf variables
@@ -19,26 +23,21 @@ NetCDF layers can have arbitrary dimensions. Known, common dimension names are
 converted to `X`, `Y` `Z`, and `Ti`, otherwise `Dim{:layername}` is used. Layers
 in the same file may also have different dimensions.
 
-NetCDF files still have issues loading directly from disk for some operations.
-Using `read(ncstack)` may help.
+## Zarr
 
-## GDAL
+Zarr files can be loaded with the [ZarrDatasets.jl](https://github.com/JuliaGeo/ZarrDatasets.jl) 
+backend. `Raster(filename; source=Zarrsource())` may be needed where the file type cant be detected 
+from the filename. `write` does not yet work for Zarr but will in future.
 
-All files GDAL can access, such as `.tiff` and `.asc` files, can be loaded,
-using [ArchGDAL.jl](https://github.com/yeesian/ArchGDAL.jl/issues). These are
-generally best loaded as `Raster("filename.tif")`, but can be loaded as
-`RasterStack("filename.tif"; layersfrom=Band)`, taking layers from the `Band`
-dimension, which is also the default.
+## GRIB
 
-## SMAP
+GRIB files can be loaded with the [ZarrDatasets.jl](https://github.com/JuliaGeo/GRIBDatasets.jl).
+`write` is not implemented for GRIB.
 
-The [Soil Moisture Active-Passive](https://smap.jpl.nasa.gov/) dataset provides
-global layers of soil moisture, temperature and other related data, in a custom
-HDF5 format. Layers are always 2 dimensional, with `Y` and `X` dimensions.
+## GRD
 
-These can be loaded as multi-layered `RasterStack("filename.h5")`. Individual
-layers can be loaded as `Raster("filename.h5"; name=:layername)`, without `name`
-the first layer is used.
+R GRD files can be loaded natively, using Julias `MMap` - which means they are very fast, but are not compressed. 
+They are always 3 dimensional, and have `Y`, `X` and [`Band`](@ref) dimensions.
 
 ````@example data_sources
 using Rasters
@@ -50,13 +49,12 @@ smapseries
 
 ## Writing file formats to disk
 
-Files can be written to disk in all formats other than SMAP HDF5 using
-`write("filename.ext", A)`. See the docs for [`write`](@ref). They can (with
-some caveats) be written to different formats than they were loaded in as,
-providing file-type conversion for spatial data.
+Files can be written to disk with ArchGDAL.jl and NCDatasets.jl backends using
+`write("filename.ext", raster)`. See the docs for [`write`](@ref). 
 
-Some metadata may be lost in formats that store little metadata, or where
-metadata conversion has not been completely implemented.
+They can (with some caveats) be written to different formats than they were loaded in as,
+providing file-type conversion for spatial data. Some metadata may be lost in formats that 
+store little metadata, or where metadata conversion has not been completely implemented.
 
 ## RasterDataSources.jl integration
 
@@ -75,4 +73,4 @@ Makie.plot(A)
 
 See the docs for [`Raster`](@ref), [`RasterStack`](@ref) and [`RasterSeries`](@ref),
 and the docs for `RasterDataSources.getraster` for syntax to specify various
-data sources.
\ No newline at end of file
+data sources.
diff --git a/docs/src/gbif_wflow.md b/docs/src/gbif_wflow.md
index de8474fac..5bd6fcedc 100644
--- a/docs/src/gbif_wflow.md
+++ b/docs/src/gbif_wflow.md
@@ -1,54 +1,119 @@
-Load occurrences for the Mountain Pygmy Possum using GBIF.jl
+# Species distribution modelling workflow
 
-## Load GBIF
+This example shows a full Species distribution modelling workflow, from loading data, to cleaning it, to fitting an ensemble and generating predictions.
 
-````@example gbif
-using Rasters, GBIF2
-using RasterDataSources
-const RS = Rasters
-````
+It uses GBIF and WorldClim data, which are common datasets in ecology.
+
+## Load Rasters, ArchGDAL, RasterDataSources and GBIF
+The GBIF2 library is used to download occurrence data, RasterDataSources to conveniently access Bioclim data. ArchGDAL is necessary to load in the Bioclim data.
 
 ````@example gbif
-records = GBIF2.occurrence_search("Burramys parvus"; limit=300)
+using Rasters, GBIF2
+using RasterDataSources, ArchGDAL
 ````
 
-## Extract coordinates
-
-Extract the longitude/latitude value to a `Vector` of points
-(a `Tuple` counts as a `(x, y)` point in GeoInterface.jl):
+Load occurrences for the Mountain Pygmy Possum using GBIF.jl
 
 ````@example gbif
-coords = [(r.decimalLongitude, r.decimalLatitude) for r in records if !ismissing(r.decimalLatitude)]
+records = GBIF2.occurrence_search("Burramys parvus"; limit=300)
 ````
 
-## Get layer / Band
-Get BioClim layers and subset to south-east Australia
+## Get Bioclimatic variables
+Get BioClim layers and subset to south-east Australia.
+The first time this is run, this will automatically download and save the files.
 
 ````@example gbif
 A = RasterStack(WorldClim{BioClim}, (1, 3, 7, 12))
-se_aus = A[X(138 .. 155), Y(-40 .. -25), RS.Band(1)]
+se_aus = A[X(138 .. 155), Y(-40 .. -25), Band(1)]
 ````
 Plot BioClim predictors and scatter occurrence points on all subplots
 
 ````@example gbif
-using Plots
-p = plot(se_aus);
-kw = (legend=:none, opacity=0.5, markershape=:cross, markercolor=:black)
-foreach(i -> scatter!(p, coords; subplot=i, kw...), 1:4)
+# The coordinates from the gbif table
+coords = collect(skipmissing(records.geometry))
+
+using CairoMakie
+p = Rasters.rplot(se_aus);
+for ax in p.content
+    if ax isa Axis
+        scatter!(ax, coords; alpha=0.5, marker='+', color=:black, markersize = 20)
+    end
+end
 p
 ````
 
-Then extract predictor variables and write to CSV.
+## Extract bioclim variables at occurrence points
+Then extract predictor variables and write to CSV. Use the skipmissing keyword to exclude both missing coordinates and coordinates with missing values in the RasterStack.
+
 ````@example gbif
 using CSV
-predictors = collect(extract(se_aus, coords))
-CSV.write("burramys_parvus_predictors.csv", predictors)
+presences = extract(se_aus, coords, skipmissing = true)
+CSV.write("burramys_parvus_predictors.csv", presences)
 ````
 
 Or convert them to a `DataFrame`:
 
 ````@example gbif
 using DataFrames
-df = DataFrame(predictors)
+df = DataFrame(presences)
 df[1:5,:]
 ````
+
+## Sample background points
+Next, sample random background points in the Raster. Rasters has a StatsBase extension to make this very straightforward. The syntax and output of `Rasters.sample` is very similar to that of `extract`.
+
+````@example gbif
+using StatsBase
+background = Rasters.sample(se_aus, 500, skipmissing = true)
+````
+
+## Fit a statistical ensemble
+In this example, we will [SpeciesDistributionModels.jl](https://github.com/tiemvanderdeure/SpeciesDistributionModels.jl) to fit a statistical ensemble to the occurrence and background data.
+
+First we need to load the models. SDM.jl integrates with MLJ - see the [model browser](https://juliaai.github.io/MLJ.jl/dev/model_browser/#Classification) for what models are available.
+
+````@example gbif
+import Maxnet: MaxnetBinaryClassifier
+import MLJGLMInterface: LinearBinaryClassifier
+# define the models in the ensemble
+models = (
+    maxnet = MaxnetBinaryClassifier(), 
+    maxnet2 = MaxnetBinaryClassifier(features = "lq"),
+    glm = LinearBinaryClassifier()
+)
+````
+
+Next, format the data using `sdmdata`. To test how rigurous our models are, we will use 3-fold cross-validation.
+
+````@example gbif
+using SpeciesDistributionModels
+const SDM = SpeciesDistributionModels
+data = sdmdata(presences, background; resampler = CV(; nfolds = 3))
+````
+
+Now, fit the ensemble, passing the data object and the `NamedTuple` of models!
+
+````@example gbif
+ensemble = sdm(data, models)
+````
+
+Use SDM.jl's evaluate function to see how this ensemble performs.
+
+````@example gbif
+SDM.evaluate(ensemble)
+````
+
+Not too bad!
+
+## Make predictions of climatic suitability
+Use the ensemble to 
+
+````@example gbif
+suitability = SDM.predict(ensemble, se_aus, reducer = mean)
+````
+
+And let's see what that looks like
+
+````@example gbif
+plot(suitability, colorrange = (0,1))
+````
diff --git a/docs/src/methods.md b/docs/src/methods.md
index fc7230331..9c1eb55ea 100644
--- a/docs/src/methods.md
+++ b/docs/src/methods.md
@@ -1,8 +1,16 @@
+## Point, polygon and table operations
+   
+| Methods                   | Description                                                                  |
+| :------------------------ | :--------------------------------------------------------------------------- |
+| [`rasterize`](@ref)       | rasterize points and geometries.                                             |
+| [`extract`](@ref)         | extract values from points or geometries.                                    |
+| [`zonal`](@ref)           | calculate zonal statistics for an object masked by geometries.               |    
+    
 ## Methods that change the resolution or extent of an object
 
 Click through to the function documentation for more in-depth descriptions and examples.
    
-| Methods |              Description                                        |
+| Methods                   | Description                                                                  |
 | :------------------------ | :--------------------------------------------------------------------------- |
 | [`aggregate`](@ref)       | aggregate data by the same or different amounts for each axis.               |
 | [`disaggregate`](@ref)    | similarly disaggregate data.                                                 |
@@ -16,32 +24,22 @@ Click through to the function documentation for more in-depth descriptions and e
     
 ## Methods that change an objects values
    
-!!! info
-      Note that most regular Julia methods, such as `replace`, work as for a standard
-      `Array`. These additional methods are commonly required in GIS applications.
-   
-| Methods |              Description                                        |
+| Methods                   | Description                                                                  |
 | :------------------------ | :--------------------------------------------------------------------------- |
 | [`classify`](@ref)        | classify values into categories.                                             |
 | [`mask`](@ref)            | mask an object by a polygon or `Raster` along `X/Y`, or other dimensions.    |
 | [`replace_missing`](@ref) | replace all missing values in an object and update `missingval`.             |
    
-   
-## Point, polygon and table operation
-   
-| Methods |              Description                                        |
-| :------------------------ | :--------------------------------------------------------------------------- |
-| [`rasterize`](@ref)       | rasterize points and geometries.                                             |
-| [`extract`](@ref)         | extract values from points or geometries.                                    |
-| [`zonal`](@ref)           | calculate zonal statistics for an object masked by geometries.               |    
+!!! info
+      Note that most regular Julia methods, such as `replace`, work as for a standard
+      `Array`. These additional methods are commonly required in GIS applications.
    
 ## Methods to load, write and modify data sources
    
-| Methods  |              Description                                   |
+| Methods                   | Description                                                             |
 | :------------------------ | :---------------------------------------------------------------------- |
 | [`modify`](@ref)          | replace the data in objects. Useful to e.g. move objects to/from a GPU. |
 | [`read`](@ref)            | read data to memory if it is on disk.                                   |
 | [`read!`](@ref)           | read data to predefined memory.                                         |
 | [`open`](@ref)            | open the underlying data for manually reading or writing.               |
-| [`write`](@ref)           | write objects to file.                                                  |
-    
+| [`write`](@ref)           | write objects to file.                                                  |
\ No newline at end of file
diff --git a/docs/src/plot_makie.md b/docs/src/plot_makie.md
index 94e82eda9..38b01eef0 100644
--- a/docs/src/plot_makie.md
+++ b/docs/src/plot_makie.md
@@ -17,14 +17,10 @@ fig, ax, _ = plot(A)
 contour(fig[1, 2], A)
 ax = Axis(fig[2, 1]; aspect = DataAspect())
 contourf!(ax, A)
-# surface(fig[2, 2], A) # even a 3D plot works! # broken mutating method
+surface(fig[2, 2], A) # even a 3D plot work!
 fig
 ````
 
-even a 3D plot works!
-````@example makie
-surface(A)
-````
 ## 3-D rasters in Makie
 
 !!! warning
@@ -59,10 +55,9 @@ Makie.set_theme!()
 `Rasters.rplot` should support Observable input out of the box, but the dimensions of that input
 must remain the same - i.e., the element names of a RasterStack must remain the same.
 
-````julia
+````@example makie
 Makie.set_theme!(Rasters.theme_rasters())
 # `stack` is the WorldClim climate data for January
-# observables not working here, same error as for contourf: ERROR: MethodError: no method matching typemax(::Type{ColorTypes.RGBA{Float32}})
 stack_obs = Observable(stack)
 fig = Rasters.rplot(stack_obs;
     Colorbar=(; height=Relative(0.75), width=5)
@@ -72,7 +67,9 @@ record(fig, "rplot.mp4", 1:12; framerate = 3) do i
 end 
 ````
 
-<!-- <video src="./rplot.mp4" controls="controls" autoplay="autoplay"></video> -->
+```@raw html
+<video src="./rplot.mp4" controls="controls" autoplay="autoplay"></video>
+```
 
 ````@example makie
 Makie.set_theme!() # reset theme
diff --git a/docs/src/plotting.md b/docs/src/plotting.md
index ffb57f9f3..99e9bcd2f 100644
--- a/docs/src/plotting.md
+++ b/docs/src/plotting.md
@@ -129,7 +129,7 @@ Then load raster data. We load some worldclim layers using `RasterDataSources` v
 
 ````@example plots
 using Rasters, RasterDataSources
-
+using Dates
 climate = RasterStack(WorldClim{Climate}, (:tmin, :tmax, :prec, :wind); month=July)
 ````
 
diff --git a/ext/RastersArchGDALExt/RastersArchGDALExt.jl b/ext/RastersArchGDALExt/RastersArchGDALExt.jl
index bf855bb9d..293013406 100644
--- a/ext/RastersArchGDALExt/RastersArchGDALExt.jl
+++ b/ext/RastersArchGDALExt/RastersArchGDALExt.jl
@@ -21,16 +21,17 @@ using Rasters: GDALsource, AbstractProjected, AbstractRaster, AbstractRasterStac
 
 import Rasters: reproject, resample, warp, cellsize, nokw, isnokw, isnokwornothing
 
+import LinearAlgebra: dot, cross
+
 const AG = ArchGDAL
 const RA = Rasters
+const AG = ArchGDAL
 const DD = DimensionalData
 const DA = DiskArrays
 const GI = GeoInterface
 const LA = Lookups
 
-include("cellsize.jl")
 include("gdal_source.jl")
-include("reproject.jl")
 include("resample.jl")
 include("warp.jl")
 
diff --git a/ext/RastersArchGDALExt/cellsize.jl b/ext/RastersArchGDALExt/cellsize.jl
deleted file mode 100644
index bfd15bdea..000000000
--- a/ext/RastersArchGDALExt/cellsize.jl
+++ /dev/null
@@ -1,87 +0,0 @@
-function _great_circle_bearing(lon1::AbstractFloat, lat1::AbstractFloat, lon2::AbstractFloat, lat2::AbstractFloat)
-    dLong = lon1 - lon2
-
-    s = cos(lat2)*sin(dLong)
-    c = cos(lat1)*sin(lat2) - sin(lat1)*cos(lat2)*cos(dLong)
-
-    return atan(s, c)
-end
-
-## Get the area of a LinearRing with coordinates in radians
-# Using Gidard's theorem
-function _area_from_rads(ring; R = 6371.0088)
-    n = GI.npoint(ring)
-    area = -(n-3)*pi
-
-    prevpoint = GI.getpoint(ring, n-1)
-    point = GI.getpoint(ring, 1)
-    
-    for i in 2:n
-        nextpoint = GI.getpoint(ring, i)
-
-        beta1 = _great_circle_bearing(GI.x(point), GI.y(point), GI.x(prevpoint), GI.y(prevpoint)) 
-        beta2 = _great_circle_bearing(GI.x(point), GI.y(point), GI.x(nextpoint), GI.y(nextpoint)) 
-        angle = acos(cos(-beta1)*cos(-beta2) + sin(-beta1)*sin(-beta2))
-        area += angle
-
-        prevpoint = point
-        point = nextpoint
-    end
-    
-    return area*R^2
-end
-
-_area_from_coords(transform, geom) = _area_from_coords(transform, GI.trait(geom), geom)
-_area_from_coords(geom) = _area_from_coords(GI.trait(geom), geom)
-function _area_from_coords(::GI.LinearRingTrait, ring) # no ArchGDAL, assumes degrees
-    points = map(GI.getpoint(ring)) do p 
-        (deg2rad(GI.x(p)), deg2rad(GI.y(p)))
-    end
-
-    return _area_from_rads(GI.LinearRing(points))
-end
-function _area_from_coords(transform::ArchGDAL.CoordTransform, ::GI.LinearRingTrait, ring)
-    points = map(GI.getpoint(ring)) do p 
-        t = ArchGDAL.transform!(ArchGDAL.createpoint(p...), transform)
-        (deg2rad(GI.x(t)), deg2rad(GI.y(t)))
-    end
-    return _area_from_rads(GI.LinearRing(points))
-end
-
-function cellsize(dims::Tuple{<:XDim, <:YDim})
-    # check the dimensions 
-    isnothing(crs(dims)) && _no_crs_error()
-    any(d -> d isa Points, sampling.(dims)) && throw(ArgumentError("Cannot calculate cell size for a `Raster` with `Points` sampling."))
-
-    xbnds, ybnds = DD.intervalbounds(dims)
-    if convert(CoordSys, crs(dims)) == CoordSys("Earth Projection 1, 104") # check if need to reproject
-        areas = [_area_from_coords(
-            GI.LinearRing([
-                (xb[1], yb[1]), 
-                (xb[2], yb[1]), 
-                (xb[2], yb[2]), 
-                (xb[1], yb[2]),
-                (xb[1], yb[1])
-            ]))
-            for xb in xbnds, yb in ybnds]
-    else 
-        areas = ArchGDAL.crs2transform(crs(dims), EPSG(4326), order = :trad) do transform
-            [_area_from_coords(
-                transform,         
-                GI.LinearRing([
-                    (xb[1], yb[1]), 
-                    (xb[2], yb[1]), 
-                    (xb[2], yb[2]), 
-                    (xb[1], yb[2]),
-                    (xb[1], yb[1])
-                ]))
-                for xb in xbnds, yb in ybnds]
-        end
-    end
-
-    return Raster(areas, dims)
-end
-
-function cellsize(x::Union{<:AbstractRaster, <:AbstractRasterStack, <:RA.DimTuple})
-    cellsize(dims(x, (XDim, YDim)))
-end
diff --git a/ext/RastersArchGDALExt/gdal_source.jl b/ext/RastersArchGDALExt/gdal_source.jl
index a906ce606..77341de3d 100644
--- a/ext/RastersArchGDALExt/gdal_source.jl
+++ b/ext/RastersArchGDALExt/gdal_source.jl
@@ -174,8 +174,8 @@ function RA._dims(raster::AG.RasterDataset, crs=nokw, mappedcrs=nokw)
         xorder = xstep > 0 ? ForwardOrdered() : ReverseOrdered()
         yorder = ystep > 0 ? ForwardOrdered() : ReverseOrdered()
         # Create lookup index. LinRange is easiest always the right size after fp error
-        xindex = LinRange(xmin, xmax, xsize)
-        yindex = LinRange(ymax, ymin, ysize)
+        xindex = range(; start=xmin, stop=xmax, length=xsize)
+        yindex = range(; start=ymax, stop=ymin, length=ysize)
 
         # Define `Projected` lookups fo X and Y dimensions
         xlookup = Projected(xindex;
diff --git a/ext/RastersArchGDALExt/reproject.jl b/ext/RastersArchGDALExt/reproject.jl
deleted file mode 100644
index e6046d999..000000000
--- a/ext/RastersArchGDALExt/reproject.jl
+++ /dev/null
@@ -1,21 +0,0 @@
-
-function Rasters._reproject(source::GeoFormat, target::GeoFormat, ::XDim, vals::AbstractVector) 
-    paired_vals = map(v -> (v, zero(v)), vals)
-    push!(paired_vals, (first(vals), one(first(vals))))
-    reps = AG.reproject(paired_vals, source, target; order=:trad)
-    at_one = pop!(reps)
-    # TODO is this test sufficient to reject the reprojection?
-    first(reps)[1] == at_one[1] || _reproject_crs_error(source, target)
-    return map(r -> r[1], reps)
-end
-function Rasters._reproject(source::GeoFormat, target::GeoFormat, dim::YDim, vals::AbstractVector)
-    paired_vals = map(v -> (zero(v), v), vals)
-    push!(paired_vals, (one(first(vals)), first(vals)))
-    reps = AG.reproject(paired_vals, source, target; order=:trad)
-    at_one = pop!(reps)
-    first(reps)[2] == at_one[2] || _reproject_crs_error(source, target)
-    return map(r -> r[2], reps)
-end
-
-_reproject_crs_error(source, target) = 
-    throw(ArgumentError("Cannot reproject from: \n $source \nto: \n $target, coordinate reference systems are not aligned with on all axes. You may need to use `resample` instead")) 
diff --git a/ext/RastersMakieExt/plotrecipes.jl b/ext/RastersMakieExt/plotrecipes.jl
index 602e22921..fac962fb9 100644
--- a/ext/RastersMakieExt/plotrecipes.jl
+++ b/ext/RastersMakieExt/plotrecipes.jl
@@ -1,8 +1,6 @@
 
 const HIDE_DEC = (; label=true, grid=false, minorgrid=false, minorticks=false)
 
-const SurfaceLikeCompat = isdefined(Makie, :SurfaceLike) ? Makie.SurfaceLike : Union{Makie.VertexGrid,Makie.CellGrid,Makie.ImageLike}
-
 function Rasters.style_rasters()
     Makie.Attributes(
         Axis=(
@@ -288,26 +286,22 @@ function Makie.plottype(raster::AbstractRaster{<:Union{Missing,Real},3})
     end
 end
 
+# ## `convert_arguments`
+# We need to handle missing values properly here, which DimensionalData.jl can't.  
+# That's why we have the `_prepare_dimarray` function, and also why this extension is necessary.
+
+# 1d
 function Makie.convert_arguments(t::Makie.PointBased, A::AbstractRaster{<:Any,1})
     return Makie.convert_arguments(t, _prepare_dimarray(A))
 end
 function Makie.convert_arguments(t::Makie.PointBased, A::AbstractRaster{<:Number,2})
     return Makie.convert_arguments(t, _prepare_dimarray(A))
 end
-@static if isdefined(Makie, :SurfaceLike)
 
-    function Makie.convert_arguments(t::SurfaceLike, A::AbstractRaster{var"#s115", 2, D}) where {var"#s115", D<:Tuple}
-        return Makie.convert_arguments(t, _prepare_dimarray(A))
-    end
-else # surfacelike is not a thing
-    Makie.convert_arguments(t::Makie.VertexGrid, A::AbstractRaster{<: Any, 2}) = Makie.convert_arguments(t, _prepare_dimarray(A))
-    Makie.convert_arguments(t::Makie.CellGrid, A::AbstractRaster{<: Any, 2}) = Makie.convert_arguments(t, _prepare_dimarray(A))
-    Makie.convert_arguments(t::Makie.ImageLike, A::AbstractRaster{<: Any, 2}) = Makie.convert_arguments(t, _prepare_dimarray(A))
-end
+Makie.convert_arguments(t::Makie.VertexGrid, A::AbstractRaster{<: Any, 2}) = Makie.convert_arguments(t, _prepare_dimarray(A))
+Makie.convert_arguments(t::Makie.CellGrid, A::AbstractRaster{<: Any, 2}) = Makie.convert_arguments(t, _prepare_dimarray(A))
+Makie.convert_arguments(t::Makie.ImageLike, A::AbstractRaster{<: Any, 2}) = Makie.convert_arguments(t, _prepare_dimarray(A))
 
-function Makie.convert_arguments(t::Makie.DiscreteSurface, A::AbstractRaster{<:Any,2})
-    return Makie.convert_arguments(t, _prepare_dimarray(A))
-end
 function Makie.convert_arguments(t::Makie.VolumeLike, A::AbstractRaster{<:Any,3}) 
     return Makie.convert_arguments(t, _prepare_dimarray(A))
 end
@@ -322,10 +316,13 @@ function Makie.convert_arguments(x::Makie.ConversionTrait, raster::AbstractRaste
         return Makie.convert_arguments(x, view(raster, rebuild(D, 1)))
     end
 end
+
 function Makie.convert_arguments(x::Makie.ConversionTrait, series::AbstractRasterSeries)
     return Makie.convert_arguments(x, first(series))
 end
 
+# ## Utility / helper functions
+
 function _series_dim(A)
     spatialdims = (X(), Y(), Z())
     last((dims(A, spatialdims)..., otherdims(A, spatialdims)...))
diff --git a/ext/RastersProjExt/RastersProjExt.jl b/ext/RastersProjExt/RastersProjExt.jl
new file mode 100644
index 000000000..bb9892946
--- /dev/null
+++ b/ext/RastersProjExt/RastersProjExt.jl
@@ -0,0 +1,33 @@
+module RastersProjExt
+
+using Rasters, Proj
+
+import DiskArrays,
+    Extents,
+    Missings
+
+using DimensionalData,
+    GeoFormatTypes,
+    GeoInterface
+
+
+using Rasters.Lookups
+using Rasters.Dimensions
+using Rasters: GDALsource, AbstractProjected, RasterStackOrArray, FileArray, NoKW,
+    RES_KEYWORD, SIZE_KEYWORD, CRS_KEYWORD, FILENAME_KEYWORD, SUFFIX_KEYWORD, EXPERIMENTAL,
+    GDAL_EMPTY_TRANSFORM, GDAL_TOPLEFT_X, GDAL_WE_RES, GDAL_ROT1, GDAL_TOPLEFT_Y, GDAL_ROT2, GDAL_NS_RES,
+    _no_crs_error
+
+import Rasters: reproject, _spherical_cellarea, nokw
+
+import LinearAlgebra: dot, cross
+
+const RA = Rasters
+const DD = DimensionalData
+const DA = DiskArrays
+const GI = GeoInterface
+const LA = Lookups
+
+include("cellarea.jl")
+include("reproject.jl")
+end
\ No newline at end of file
diff --git a/ext/RastersProjExt/cellarea.jl b/ext/RastersProjExt/cellarea.jl
new file mode 100644
index 000000000..7b2cda59d
--- /dev/null
+++ b/ext/RastersProjExt/cellarea.jl
@@ -0,0 +1,187 @@
+## Get the area of a LinearRing with coordinates in radians
+struct SphericalPoint{T <: Real}
+	data::NTuple{3, T}
+end
+SphericalPoint(x, y, z) = SphericalPoint((x, y, z))
+
+# define the 4 basic mathematical operators elementwise on the data tuple
+Base.:+(p::SphericalPoint, q::SphericalPoint) = SphericalPoint(p.data .+ q.data)
+Base.:-(p::SphericalPoint, q::SphericalPoint) = SphericalPoint(p.data .- q.data)
+Base.:*(p::SphericalPoint, q::SphericalPoint) = SphericalPoint(p.data .* q.data)
+Base.:/(p::SphericalPoint, q::SphericalPoint) = SphericalPoint(p.data ./ q.data)
+# Define sum on a SphericalPoint to sum across its data
+Base.sum(p::SphericalPoint) = sum(p.data)
+
+# define dot and cross products
+dot(p::SphericalPoint, q::SphericalPoint) = sum(p * q)
+function cross(a::SphericalPoint, b::SphericalPoint)
+	a1, a2, a3 = a.data
+    b1, b2, b3 = b.data
+	SphericalPoint((a2*b3-a3*b2, a3*b1-a1*b3, a1*b2-a2*b1))
+end
+
+function _spherical_quadrilateral_area(ring)
+    # don't assume the ring is a GI ring
+    (p1, p2, p3, p4) = _lonlat_to_sphericalpoint.(ring)
+    area = 0.0
+    area += _spherical_triangle_area(p1, p2, p3)
+    area += _spherical_triangle_area(p3, p4, p1)
+end
+
+# Using Eriksson's formula for the area of spherical triangles: https://www.jstor.org/stable/2691141
+function _spherical_triangle_area(a, b, c)
+    #t = abs(dot(a, cross(b, c)))
+    #t /= 1 + dot(b,c) + dot(c, a) + dot(a, b)
+    t = abs(dot(a, (cross(b - a, c - a))) / dot(b + a, c + a))
+    2*atan(t)
+end
+
+_lonlat_to_sphericalpoint(args) = _lonlat_to_sphericalpoint(args...)
+function _lonlat_to_sphericalpoint(lon, lat)
+    lonsin, loncos = sincosd(lon)
+    latsin, latcos = sincosd(lat)
+    x = latcos * loncos
+    y = latcos * lonsin
+    z = latsin
+    return SphericalPoint(x,y,z)
+end
+
+# _area_from_coords(transform, geom) = _area_from_coords(transform, geom)
+function _area_from_coords(transform::Proj.Transformation, ring_points)
+    t = (transform(GI.x(p), GI.y(p)) for p in ring_points)
+    return _spherical_quadrilateral_area(t)
+end
+
+# For lat-lon projections. Get the area of each latitudinal band, then multiply by the width
+function _area_from_lonlat(lon::XDim, lat::YDim; radius)
+    two_pi_R2 = 2 * pi * radius * radius
+    band_area = broadcast(DD.intervalbounds(lat)) do yb
+        two_pi_R2 * (sind(yb[2]) - sind(yb[1]))
+    end
+    
+    broadcast(DD.intervalbounds(lon), band_area') do xb, ba
+        abs((xb[2] - xb[1]) / 360 * ba)
+    end
+end
+
+function _spherical_cellarea(dims::Tuple{<:XDim, <:YDim}; radius = 6371008.8, use_area_of_use = true)
+    # check the dimensions 
+    isnothing(crs(dims)) && _no_crs_error()
+
+    areas = if _isdegrees(crs(dims)) # check if need to reproject
+        _area_from_lonlat(dims...; radius)
+    elseif !isnothing(mappedcrs(dims)) && _isdegrees(mappedcrs(dims))
+        _area_from_lonlat(reproject(dims; crs = mappedcrs(dims))...; radius)
+    else
+        xbnds, ybnds = DD.intervalbounds(dims)
+        R2 = radius * radius
+        area_of_use = if use_area_of_use
+            # Use a temporary transformation object to get the area of use,
+            # since the transformation object has to be recreated with the area of use
+            _get_area_of_use(Proj.Transformation(crs(dims), EPSG(4326); always_xy = true), Extents.Extent(X = extrema(dims[1]), Y = extrema(dims[2])))
+        else
+            C_NULL
+        end
+
+        transform = Proj.Transformation(crs(dims), EPSG(4326); always_xy = true, area = area_of_use)
+
+        result = [
+                _area_from_coords(
+                    transform,         
+                    (
+                        (xb[1], yb[1]), 
+                        (xb[2], yb[1]), 
+                        (xb[2], yb[2]), 
+                        (xb[1], yb[2]),
+                    )
+                ) * R2
+                for xb in xbnds, yb in ybnds
+            ]
+        
+        if use_area_of_use
+            Proj.proj_area_destroy(area_of_use)
+        end
+
+        return result
+    end
+end
+
+function _get_area_of_use(transform::Proj.Transformation, extent::Extents.Extent; densify_pts = 21)
+    # Transform the extent using `proj_trans_bounds`
+    (xmin, xmax), (ymin, ymax) = Proj.bounds(transform, extent.X, extent.Y; densify_pts)
+
+    # Create an area of use object 
+    # This MUST be destroyed by the caller
+    area = Proj.proj_area_create()
+
+    # Set the bounding box of the area of use
+    Proj.proj_area_set_bbox(
+        area, 
+        xmin, # west_lon_degree
+        ymin, # south_lat_degree
+        xmax, # east_lon_degree
+        ymax  # north_lat_degree
+    )
+
+    return area
+end
+
+# TODO: put these in Proj (specifically the dispatches on GFT types)
+_isgeographic(crs) = _isgeographic(Proj.CRS(crs))
+_isgeographic(crs::Proj.CRS) = Proj.is_geographic(crs)
+
+_isdegrees(crs) = _isdegrees(Proj.CRS(crs))
+function _isdegrees(crs::Proj.CRS)
+    _isgeographic(crs) || return false
+    # This is a tiny bit inefficient, but it takes 500ns on my machine, 
+    # so I think we can disregard the inefficiency...
+    return axis_is_degrees(crs, 0) && axis_is_degrees(crs, 1)
+end
+
+function axis_is_degrees(crs::Proj.CRS, axis_index::Int; context::Ptr{Proj.PJ_CONTEXT} = C_NULL)
+    @assert axis_index in (0, 1)
+    # In Proj, the `CoordinateSystem` object is contained within the `CRS` object.
+    # So we need to extract it explicitly, since Proj doesn't provide utilities for this.
+    cs = Proj.proj_crs_get_coordinate_system(crs.pj, context)
+
+    # Instantiate some refs that we'll use to get information out of the PJ struct
+    auth_name = Ref{Cstring}()
+    code = Ref{Cstring}()
+    unitname = Ref{Cstring}()
+    
+    # Load unit info for the given axis into the pointers
+    Proj.proj_cs_get_axis_info(
+        cs, 
+        axis_index, 
+        C_NULL, # out_name
+        C_NULL, # out_abbrev
+        C_NULL, # out_direction
+        C_NULL, # out_unit_conv_factor
+        C_NULL, # out_unit_name
+        auth_name,
+        code,
+        context
+    )
+
+    # We don't `unsafe_string` the C strings, because we're just going to pass them to Proj's unit lookup function
+    Proj.proj_uom_get_info_from_database(
+        auth_name[],
+        code[],
+        unitname,   # out_name
+        C_NULL,     # out_conv_factor
+        C_NULL,     # out_category
+        context
+    )
+    # Destroy the coordinate system object
+    Proj.proj_destroy(cs)
+
+    unit_str = unsafe_string(unitname[])
+
+    # TODO: check if this is the correct way to check if the unit is degrees
+    # We can also check if the unit category is "angular", but I chose to replicate
+    # the original test from ArchGDAL here.
+    # If the unit is not "degree", we could still have an angular unit (radians or some linearly scaled thing),
+    # in which case we should technically return true.
+    # We'd also have to return the conversion factor in this case, and maybe the category (radians or degrees)...
+    return isequal(unit_str, "degree")
+end
diff --git a/ext/RastersProjExt/reproject.jl b/ext/RastersProjExt/reproject.jl
new file mode 100644
index 000000000..52f9eb471
--- /dev/null
+++ b/ext/RastersProjExt/reproject.jl
@@ -0,0 +1,53 @@
+
+Rasters._reproject(source::GeoFormat, target::GeoFormat, dim::XDim, vals::AbstractVector) = _reproject!(source, target, dim, Float64.(vals))
+Rasters._reproject(source::GeoFormat, target::GeoFormat, dim::YDim, vals::AbstractVector) = _reproject!(source, target, dim, Float64.(vals))
+
+function _reproject!(source::GeoFormat, target::GeoFormat, ::XDim, vals::Vector{Float64}) 
+    # First, construct a transformation from `source` to `target`.
+    # TODO: add area of use, so that the transformation is more accurate.
+    trans = Proj.Transformation(source, target; always_xy = true, direction = Proj.PJ_FWD)
+    # Check that the reprojection is valid, i.e., trans((x, 0))[1] == trans((x, 1))[1]
+    first_x = first(vals)
+    zero_y = (first_x, zero(first_x))
+    one_y = (first_x, one(first_x))
+    # Check that the reprojection is valid, i.e., trans((x, 0))[1] == trans((x, 1))[1]
+    trans(zero_y)[1] == trans(one_y)[1] || _reproject_crs_error(source, target)
+    # Now that we've proved that the reprojection is valid, we can proceed.
+    # Here, for efficiency, and to avoid allocations, we'll use `proj_trans_generic`,
+    # which mutates values in-place.
+    Proj.proj_trans_generic(
+        trans.pj,                                    # pointer to transformation
+        Proj.PJ_FWD,                                 # direction of transformation
+        vals, sizeof(typeof(first_x)), length(vals), # input x values
+        C_NULL, 0, 0,                                # default (assumed 0) y values
+        C_NULL, 0, 0,                                # default (assumed 0) z values
+        C_NULL, 0, 0                                 # default (assumed 0) t values
+    )
+    return vals
+end
+function _reproject!(source::GeoFormat, target::GeoFormat, dim::YDim, vals::AbstractVector)
+    # First, construct a transformation from `source` to `target`.
+    # TODO: add area of use, so that the transformation is more accurate.
+    trans = Proj.Transformation(source, target; always_xy = true, direction = Proj.PJ_FWD)
+    # Check that the reprojection is valid, i.e., trans((0, y))[2] == trans((1, y))[2]
+    first_y = first(vals)
+    zero_x = (zero(first_y), first_y)
+    one_x = (one(first_y), first_y)
+    # Check that the reprojection is valid, i.e., trans((0, y))[2] == trans((1, y))[2]
+    trans(zero_x)[2] == trans(one_x)[2] || _reproject_crs_error(source, target)
+    # Now that we've proved that the reprojection is valid, we can proceed.
+    # Here, for efficiency, and to avoid allocations, we'll use `proj_trans_generic`,
+    # which mutates values in-place.
+    Proj.proj_trans_generic(
+        trans.pj,                                    # pointer to transformation
+        Proj.PJ_FWD,                                 # direction of transformation
+        C_NULL, 0, 0,                                # default (assumed 0) x values
+        vals, sizeof(typeof(first_y)), length(vals), # input y values
+        C_NULL, 0, 0,                                # default (assumed 0) z values
+        C_NULL, 0, 0                                 # default (assumed 0) t values
+    )
+    return vals
+end
+
+_reproject_crs_error(source, target) = 
+    throw(ArgumentError("Cannot reproject from: \n $source \nto: \n $target, coordinate reference systems are not aligned on all axes. You may need to use `resample` instead")) 
diff --git a/ext/RastersRasterDataSourcesExt/constructors.jl b/ext/RastersRasterDataSourcesExt/constructors.jl
index 24bafb87c..b6623402c 100644
--- a/ext/RastersRasterDataSourcesExt/constructors.jl
+++ b/ext/RastersRasterDataSourcesExt/constructors.jl
@@ -26,7 +26,21 @@ function RA.Raster(T::Type{<:RDS.RasterDataSource}, layer; crs=_source_crs(T), k
     filename = getraster(T, layer; rds_kw...)
     Raster(filename; name=RDS.layerkeys(T, layer), crs, ra_kw...)
 end
-
+function RA.Raster(T::Type{<:WorldClim{<:Future{BioClim, CMIP6}}}; crs=_source_crs(T), kw...)
+    rds_kw, ra_kw = _filterkw(T, kw)
+    filename = getraster(T; rds_kw...)
+    Raster(filename; crs, ra_kw...)
+end
+function RA.Raster(T::Type{<:WorldClim{<:Future{BioClim, CMIP6}}}, layer; crs=_source_crs(T), lazy = false, kw...)
+    rds_kw, ra_kw = _filterkw(T, kw)
+    filename = getraster(T, layer; rds_kw...)
+    ras_all = Raster(filename; name=RDS.bioclim_key(layer), crs, ra_kw...)
+    ras = view(ras_all, Band(RDS.bioclim_int(layer)))
+    if ~lazy
+        read(ras)
+    end
+    return ras
+end
 """
     RasterStack(T::Type{<:RasterDataSource}, [layers::Union{Symbol,AbstractArray,Tuple}]; kw...) => RasterStack
 
@@ -57,6 +71,18 @@ function RA.RasterStack(T::Type{<:RDS.RasterDataSource}, layers::Tuple; crs=_sou
     filenames = map(l -> RDS.getraster(T, l; rds_kw...), layers)
     RasterStack(filenames; name=RDS.layerkeys(T, layers), crs, ra_kw...)
 end
+# WorldClim future has all layers in one .tif file - format this as a rasterdatastack in a dedicated dispatch
+function RA.RasterStack(T::Type{<:WorldClim{<:Future{BioClim, CMIP6}}}, layers::Tuple; crs=_source_crs(T), lazy = false, kw...)
+    rds_kw, ra_kw = _filterkw(T, kw)
+    keys = map(RDS.bioclim_key, layers)
+    filename = RDS.getraster(T, layers; rds_kw...)
+    raster = Raster(filename; lazy = true, ra_kw...)
+    stack = RasterStack(eachslice(raster; dims = Band)...; name=RDS.layerkeys(T), crs)[keys]
+    if ~lazy
+        stack = read(stack)
+    end
+    return stack
+end
 
 """
     RasterSeries(T::Type{<:RasterDataSource}, [layers::Union{Symbol,AbstractArray,Tuple}]; kw...) => AbstractRasterSeries
diff --git a/ext/RastersStatsBaseExt/RastersStatsBaseExt.jl b/ext/RastersStatsBaseExt/RastersStatsBaseExt.jl
new file mode 100644
index 000000000..fd6f1e918
--- /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 000000000..f54f9d919
--- /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/Rasters.jl b/src/Rasters.jl
index a4d4a7921..e660e7c86 100644
--- a/src/Rasters.jl
+++ b/src/Rasters.jl
@@ -21,29 +21,33 @@ import Adapt,
        FillArrays,
        Flatten,
        GeoInterface,
+       GeometryOpsCore,
        OffsetArrays,
        ProgressMeter,
        Missings,
        Mmap,
        RecipesBase,
        Reexport,
-       Setfield
+       Setfield,
+       Statistics
 
 Reexport.@reexport using DimensionalData, GeoFormatTypes
 
 using DimensionalData.Tables,
       DimensionalData.Lookups,
-      DimensionalData.Dimensions
+      DimensionalData.Dimensions,
       DimensionalData.Lookups.IntervalSets
 
 using DimensionalData: Name, NoName
 using .Dimensions: StandardIndices, DimTuple
 using .Lookups: LookupTuple
 
+using Statistics: mean
 using RecipesBase: @recipe, @series
 using Base: tail, @propagate_inbounds
 
 import GeoInterface: crs
+import Extents: Extent, extent
 
 using Setfield: @set, @set!
 using ColorTypes: RGB
@@ -52,6 +56,9 @@ using CommonDataModel: AbstractDataset, AbstractVariable
 
 using DiskArrays: @implement_diskarray
 
+using GeometryOpsCore: Planar, Spherical
+export Planar, Spherical
+
 export AbstractRaster, Raster
 export AbstractRasterStack, RasterStack
 export AbstractRasterSeries, RasterSeries
@@ -61,10 +68,10 @@ export missingval, boolmask, missingmask, replace_missing, replace_missing!,
        aggregate, aggregate!, disaggregate, disaggregate!, mask, mask!,
        resample, warp, zonal, crop, extend, trim, slice, combine, points,
        classify, classify!, mosaic, mosaic!, extract, rasterize, rasterize!,
-       coverage, coverage!, setcrs, setmappedcrs, smapseries, cellsize
+       coverage, coverage!, setcrs, setmappedcrs, smapseries, cellsize, cellarea
 export crs, mappedcrs, mappedindex, mappedbounds, projectedindex, projectedbounds
 export reproject, convertlookup
-
+export Extent, extent
 
 const DD = DimensionalData
 const DA = DiskArrays
diff --git a/src/array.jl b/src/array.jl
index 1fdd590b2..32d356376 100644
--- a/src/array.jl
+++ b/src/array.jl
@@ -104,20 +104,17 @@ function DD.modify(f, A::AbstractRaster)
     return rebuild(A, newdata)
 end
 
-function DD.DimTable(As::Tuple{<:AbstractRaster,Vararg{<:AbstractRaster}}...)
-    DD.DimTable(DimStack(map(read, As...)))
-end
+DD.DimTable(As::Tuple{<:AbstractRaster,Vararg{AbstractRaster}}) =
+    DD.DimTable(DimStack(map(read, As)))
 
 # DiskArrays methods
 
 DiskArrays.eachchunk(A::AbstractRaster) = DiskArrays.eachchunk(parent(A))
 DiskArrays.haschunks(A::AbstractRaster) = DiskArrays.haschunks(parent(A))
-function DA.readblock!(A::AbstractRaster, dst, r::AbstractUnitRange...)
+DA.readblock!(A::AbstractRaster, dst, r::AbstractUnitRange...) =
     DA.readblock!(parent(A), dst, r...)
-end
-function DA.writeblock!(A::AbstractRaster, src, r::AbstractUnitRange...)
+DA.writeblock!(A::AbstractRaster, src, r::AbstractUnitRange...) =
     DA.writeblock!(parent(A), src, r...)
-end
 
 # Base methods
 
@@ -260,6 +257,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...
@@ -354,4 +352,4 @@ end
 filekey(ds, name) = name
 filekey(filename::String) = Symbol(splitext(basename(filename))[1])
 
-DD.dimconstructor(::Tuple{<:Dimension{<:AbstractProjected},Vararg{<:Dimension}}) = Raster
+DD.dimconstructor(::Tuple{<:Dimension{<:AbstractProjected},Vararg{Dimension}}) = Raster
\ No newline at end of file
diff --git a/src/extensions.jl b/src/extensions.jl
index 36f20a024..b335fdadc 100644
--- a/src/extensions.jl
+++ b/src/extensions.jl
@@ -159,44 +159,121 @@ $EXPERIMENTAL
 """
 warp(args...; kw...) = throw_extension_error(warp, "ArchGDAL", :RastersArchGDALExt, args)
 
-"""
-    cellsize(x)
+# stubs that need Proj
 
-Gives the approximate size of each cell in square km.
-This function works for any projection, using an algorithm for polygons on a sphere. It approximates the true size to about 0.1%, depending on latitude.
+"""
+    cellarea([method], x)
 
-Run `using ArchGDAL` to make this method available.
+Gives the approximate area of each gridcell of `x`.
+By assuming the earth is a sphere, it approximates the true size to about 0.1%, depending on latitude.
 
-# Arguments
+Run `using ArchGDAL` to make this method fully available.
 
-- `x`: A `Raster` or a `Tuple` of `X` and `Y` dimensions.
+`method` can be `Spherical(; radius)` (the default) or `Planar()`.
+- `Spherical` will compute cell area on the sphere, by transforming all points back to long-lat.  You can specify the radius by the `radius` keyword argument here.  By default, this is `6371008.8`, the mean radius of the Earth.
+- `Planar` will compute cell area in the plane of the CRS you have chosen.  Be warned that this will likely be incorrect for non-equal-area projections.
 
 ## Example
 
 ```julia
 using Rasters, ArchGDAL, Rasters.Lookups
-dimz = X(Projected(90.0:10.0:120; sampling=Intervals(Start()), order=ForwardOrdered(), span=Regular(10.0), crs=EPSG(4326))),
-       Y(Projected(0.0:10.0:50; sampling=Intervals(Start()), order=ForwardOrdered(), span=Regular(10.0), crs=EPSG(4326)))
-
-cs = cellsize(dimz)
+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)
+cs = cellarea(myraster)
 
 # output
-4×6 Raster{Float64,2} with dimensions:
-  X Projected{Float64} 90.0:10.0:120.0 ForwardOrdered Regular Intervals{Start} crs: EPSG,
-  Y Projected{Float64} 0.0:10.0:50.0 ForwardOrdered Regular Intervals{Start} crs: EPSG
-extent: Extent(X = (90.0, 130.0), Y = (0.0, 60.0))
-missingval: missing
-crs: EPSG:4326
-parent:
-        0.0       10.0       20.0        30.0        40.0            50.0
-  90.0  1.2332e6   1.1952e6   1.12048e6   1.01158e6   8.72085e5  706488.0
- 100.0  1.2332e6   1.1952e6   1.12048e6   1.01158e6   8.72085e5  706488.0
- 110.0  1.2332e6   1.1952e6   1.12048e6   1.01158e6   8.72085e5  706488.0
- 120.0  1.2332e6   1.1952e6   1.12048e6   1.01158e6   8.72085e5  706488.0
+╭───────────────────────╮
+│ 4×6 Raster{Float64,2} │
+├───────────────────────┴─────────────────────────────────────────────────── dims ┐
+  ↓ X Projected{Float64} 90.0:10.0:120.0 ForwardOrdered Regular Intervals{Start},
+  → Y Projected{Float64} 0.0:10.0:50.0 ForwardOrdered Regular Intervals{Start}
+├───────────────────────────────────────────────────────────────────────── raster ┤
+  extent: Extent(X = (90.0, 130.0), Y = (0.0, 60.0))
+
+  crs: EPSG:4326
+└─────────────────────────────────────────────────────────────────────────────────┘
+   ↓ →  0.0        10.0        20.0        30.0            40.0      50.0
+  90.0  1.23017e6   1.19279e6   1.11917e6   1.01154e6  873182.0  708290.0
+ 100.0  1.23017e6   1.19279e6   1.11917e6   1.01154e6  873182.0  708290.0
+ 110.0  1.23017e6   1.19279e6   1.11917e6   1.01154e6  873182.0  708290.0
+ 120.0  1.23017e6   1.19279e6   1.11917e6   1.01154e6  873182.0  708290.0
 ```
 $EXPERIMENTAL
 """
-cellsize(args...; kw...) = throw_extension_error(cellsize, "ArchGDAL", :RastersArchGDALExt, args)
+cellarea(x; kw...) = cellarea(Spherical(), x; kw...)
+cellarea(method::GeometryOpsCore.Manifold, x; kw...) = cellarea(method, dims(x, (XDim, YDim)); kw...)
+
+function cellarea(method::GeometryOpsCore.Planar, dims::Tuple{<:XDim, <:YDim}; kw...)
+    isintervals(dims) || throw(ArgumentError("Cannot calculate cell size for a `Raster` with `Points` sampling."))
+    areas = _planar_cellarea(dims; kw...)
+    return Raster(areas; dims)
+end
+
+function cellarea(method::GeometryOpsCore.Spherical, dims::Tuple{<:XDim, <:YDim}; kw...)
+    isintervals(dims) || throw(ArgumentError("Cannot calculate cell size for a `Raster` with `Points` sampling."))
+    areas = _spherical_cellarea(dims; radius = method.radius, kw...)
+    return Raster(areas; dims)
+end
+
+_spherical_cellarea(args...; kw...) = throw_extension_error(_spherical_cellarea, "Proj", :RastersProjExt, args)
+
+function _planar_cellarea(dims::Tuple{<:XDim, <:YDim})
+    xbnds, ybnds = DD.intervalbounds(dims)
+    broadcast(xb -> xb[2] - xb[1], xbnds) .* broadcast(yb -> yb[2] - yb[1], ybnds)'
+end
+
+function cellsize(args...; kw...) 
+    @warn """
+cellsize is deprecated and will be removed in a future version, use cellarea instead. 
+Note that cellarea returns the area in square m, while cellsize still uses square km.
+"""
+    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
diff --git a/src/lookup.jl b/src/lookup.jl
index 0c0c05fa1..07d1a68e0 100644
--- a/src/lookup.jl
+++ b/src/lookup.jl
@@ -89,18 +89,30 @@ GeoInterface.crs(lookup::Projected) = lookup.crs
 mappedcrs(lookup::Projected) = lookup.mappedcrs
 dim(lookup::Projected) = lookup.dim
 
-function LA.selectindices(l::Projected, sel::Contains)
+@inline function LA.selectindices(l::Projected, sel::LA.Selector; kw...)
     selval = reproject(mappedcrs(l), crs(l), dim(l), val(sel))
-    LA.contains(l, rebuild(sel; val=selval))
+    LA._selectindices(l, rebuild(sel; val=selval); kw...)
 end
-function LA.selectindices(l::Projected, sel::At)
+@inline function LA.selectindices(l::Projected, sel::LA.Selector{<:AbstractVector}; kw...)
     selval = reproject(mappedcrs(l), crs(l), dim(l), val(sel))
-    LA.at(l, rebuild(sel; val=selval))
+    _selectvec(l, rebuild(sel; val=selval)sel; kw...)
 end
-function LA.selectindices(l::Projected, sel::Between)
+@inline function LA.selectindices(l::Projected, sel::LA.IntSelector{<:Tuple}; kw...)
+    selval = reproject(mappedcrs(l), crs(l), dim(l), val(sel))
+    _selecttuple(l, rebuild(sel; val=selval); kw...)
+end
+@inline LA.selectindices(l::Projected{<:Tuple}, sel::LA.IntSelector{<:Tuple}; kw...) = LA._selectindices(l, sel; kw...)
+@inline LA.selectindices(l::Projected{<:Tuple}, sel::LA.IntSelector{<:Tuple{<:Tuple,<:Tuple}}; kw...) = 
+    _selecttuple(l, sel; kw...)
+
+function LA.selectindices(l::Projected, sel::Between{<:Tuple})
     selval = map(v -> reproject(mappedcrs(l), crs(l), dim(l), v), val(sel))
     LA.between(l, rebuild(sel; val=selval))
 end
+function LA.selectindices(l::Projected, sel::T) where T<:DD.IntervalSets.Interval
+    left, right = map(v -> reproject(mappedcrs(l), crs(l), dim(l), v), (sel.left, sel.right))
+    LA.between(l, T(left, right))
+end
 
 """
     Mapped <: AbstractProjected
diff --git a/src/methods/aggregate.jl b/src/methods/aggregate.jl
index ffa5baabd..91686fa83 100644
--- a/src/methods/aggregate.jl
+++ b/src/methods/aggregate.jl
@@ -1,6 +1,6 @@
 
-const DimOrDimTuple = Union{Dimension,Tuple{Vararg{<:Dimension}}}
-const IntOrIntTuple = Union{Int,Tuple{Vararg{<:Int}}}
+const DimOrDimTuple = Union{Dimension,Tuple{Vararg{Dimension}}}
+const IntOrIntTuple = Union{Int,Tuple{Vararg{Int}}}
 
 struct Ag end
 struct DisAg end
@@ -38,6 +38,7 @@ $PROGRESS_KEYWORD
 
 ```jldoctest
 using Rasters, RasterDataSources, Statistics, Plots
+import ArchGDAL
 using Rasters: Center
 st = read(RasterStack(WorldClim{Climate}; month=1))
 ag = aggregate(Center(), st, (Y(20), X(20)); skipmissingval=true, progress=false)
@@ -150,7 +151,7 @@ function aggregate!(f, dst::AbstractRaster, src, scale; skipmissingval=false)
         end
         if skipmissingval
             # All missing values return a missing value
-            if all(map(x -> x === missingval(src), block))
+            if all(x -> x === missingval(src), block)
                 _missingval_or_missing(dst)
             else
                 # Skip missing values
@@ -158,7 +159,7 @@ function aggregate!(f, dst::AbstractRaster, src, scale; skipmissingval=false)
             end
         else
             # Any missing values return a missing value
-            if any(map(x -> x === missingval(src), block))
+            if any(x -> x === missingval(src), block)
                 _missingval_or_missing(dst)
             else
                 f(block)
diff --git a/src/methods/burning/edges.jl b/src/methods/burning/edges.jl
index 58e34319c..50ee09c29 100644
--- a/src/methods/burning/edges.jl
+++ b/src/methods/burning/edges.jl
@@ -42,6 +42,7 @@ _x_at_y(e::Edge, y) = (y - e.start[2]) * e.gradient + e.start[1]
 struct Edges <: AbstractVector{Edge}
     edges::Vector{Edge}
     max_ylen::Int
+    min_y::Int
     edge_count::Int
 end
 Edges(geom, dims; kw...) = Edges(GI.geomtrait(geom), geom, dims; kw...)
@@ -55,12 +56,14 @@ function Edges(
     # TODO fix bug that requires this to be redefined
     edges = Vector{Edge}(undef, 0)
     local edge_count = max_ylen = 0
+    local min_y = typemax(Int)
     if tr isa GI.AbstractCurveTrait
-        edge_count, max_ylen = _to_edges!(edges, geom, dims, edge_count)
+        edge_count, max_ylen, min_y = _to_edges!(edges, geom, dims, edge_count)
     else
         for ring in GI.getring(geom)
-             edge_count, ring_max_ylen = _to_edges!(edges, ring, dims, edge_count)
+             edge_count, ring_max_ylen, ring_min_y = _to_edges!(edges, ring, dims, edge_count)
              max_ylen = max(max_ylen, ring_max_ylen)
+             min_y = min(min_y, ring_min_y)
         end
     end
 
@@ -72,7 +75,7 @@ function Edges(
         sort!(edges1; scratch)
     end
 
-    return Edges(edges, max_ylen, edge_count)
+    return Edges(edges, max_ylen, min_y, edge_count)
 end
 
 Base.parent(edges::Edges) = edges.edges
@@ -99,8 +102,9 @@ end
     local firstpos = prevpos = nextpos = Position((0.0, 0.0), 0)
     isfirst = true
     local max_ylen = 0
+    local min_y = typemax(Int)
 
-    GI.npoint(geom) > 0 || return edge_count, max_ylen
+    GI.npoint(geom) > 0 || return edge_count, max_ylen, min_y
     xlookup, ylookup = lookup(dims, (X(), Y())) 
     (length(xlookup) > 0 && length(ylookup) > 0) || return edge_count, max_ylen
 
@@ -136,6 +140,7 @@ end
         edge = Edge(prevpos, nextpos)
         _add_edge!(edges, edge, edge_count)
         max_ylen = max(max_ylen, edge.iystop - edge.iystart)
+        min_y = min(min_y, edge.iystart)
         prevpos = nextpos
         prevpoint = p
     end
@@ -145,11 +150,12 @@ end
         edge = Edge(prevpos, firstpos)
         # Update the longest y distance of any edge
         max_ylen = max(max_ylen, edge.iystop - edge.iystart)
+        min_y = min(min_y, edge.iystart)
         # assign/push the edge to edges
         _add_edge!(edges, edge, edge_count)
     end
 
-    return edge_count, max_ylen
+    return edge_count, max_ylen, min_y
 end
 
 function _add_edge!(edges, edge, edge_count)
diff --git a/src/methods/burning/line.jl b/src/methods/burning/line.jl
index f7932589d..9155479bb 100644
--- a/src/methods/burning/line.jl
+++ b/src/methods/burning/line.jl
@@ -178,7 +178,7 @@ function _burn_line!(A::AbstractRaster, line, fill)
     end
     return n_on_line
 end
-function _burn_line!(A::AbstractRaster, line, fill, order::Tuple{Vararg{<:Dimension}})
+function _burn_line!(A::AbstractRaster, line, fill, order::Tuple{Vararg{Dimension}})
     msg = """"
         Converting a `:line` geometry to raster is currently only implemented for 2d lines.
         Make a Rasters.jl github issue if you need this for more dimensions.
diff --git a/src/methods/burning/polygon.jl b/src/methods/burning/polygon.jl
index 2def91428..f6d4512e4 100644
--- a/src/methods/burning/polygon.jl
+++ b/src/methods/burning/polygon.jl
@@ -67,11 +67,11 @@ function _set_crossings!(crossings::Vector{Float64}, edges::Edges, iy::Int, prev
     # max_ylen tells us how big the largest y edge is.
     # We can use this to jump back from the last y position
     # rather than iterating from the start of the edges
-    ypos = max(1, prev_ypos - edges.max_ylen - 1)
-    ncrossings = 0
-    # We know the maximum size on y, so we can start from ypos 
+    ypos = max(edges.min_y, prev_ypos - edges.max_ylen - 1)
     start_ypos = searchsortedfirst(edges, ypos)
+    # We know the maximum size on y, so we can start from ypos 
     prev_ypos = start_ypos
+    ncrossings = 0
     for i in start_ypos:lastindex(edges)
         e = @inbounds edges[i]
         # Edges are sorted on y, so we can skip
diff --git a/src/methods/burning/utils.jl b/src/methods/burning/utils.jl
index d51af3c8f..0e87d3e73 100644
--- a/src/methods/burning/utils.jl
+++ b/src/methods/burning/utils.jl
@@ -4,16 +4,16 @@ _geomindices(::GI.FeatureCollectionTrait, geoms) = 1:GI.nfeature(geoms)
 _geomindices(::GI.FeatureTrait, geoms) = _geomindices(GI.geometry(geoms))
 _geomindices(::GI.AbstractGeometryTrait, geoms) = 1:GI.ngeom(geoms)
 
-_getgeom(geoms, i::Integer) = _getgeom(GI.trait(geoms), geoms, i)
-_getgeom(::GI.FeatureCollectionTrait, geoms, i::Integer) = GI.geometry(GI.getfeature(geoms, i))
-_getgeom(::GI.FeatureTrait, geoms, i::Integer) = GI.getgeom(GI.geometry(geoms), i)
-_getgeom(::GI.AbstractGeometryTrait, geoms, i::Integer) = GI.getgeom(geom, i)
-_getgeom(::GI.PointTrait, geom, i::Integer) = error("PointTrait should not be reached")
-_getgeom(::Nothing, geoms, i::Integer) = geoms[i] # Otherwise we can probably just index?
+_getgeom(geoms, i::Integer) = __getgeom(GI.trait(geoms), geoms, i)
+__getgeom(::GI.FeatureCollectionTrait, geoms, i::Integer) = GI.geometry(GI.getfeature(geoms, i))
+__getgeom(::GI.FeatureTrait, geoms, i::Integer) = GI.getgeom(GI.geometry(geoms), i)
+__getgeom(::GI.AbstractGeometryTrait, geoms, i::Integer) = GI.getgeom(geom, i)
+__getgeom(::GI.PointTrait, geom, i::Integer) = error("PointTrait should not be reached")
+__getgeom(::Nothing, geoms, i::Integer) = geoms[i] # Otherwise we can probably just index?
 
-_getgeom(geoms) = _getgeom(GI.trait(geoms), geoms)
-_getgeom(::GI.FeatureCollectionTrait, geoms) = (GI.geometry(f) for f in GI.getfeature(geoms))
-_getgeom(::GI.FeatureTrait, geoms) = GI.getgeom(GI.geometry(geoms))
-_getgeom(::GI.AbstractGeometryTrait, geoms) = GI.getgeom(geoms)
-_getgeom(::GI.PointTrait, geom) = error("PointTrait should not be reached")
-_getgeom(::Nothing, geoms) = geoms
+_getgeom(geoms) = __getgeom(GI.trait(geoms), geoms)
+__getgeom(::GI.FeatureCollectionTrait, geoms) = (GI.geometry(f) for f in GI.getfeature(geoms))
+__getgeom(::GI.FeatureTrait, geoms) = GI.getgeom(GI.geometry(geoms))
+__getgeom(::GI.AbstractGeometryTrait, geoms) = GI.getgeom(geoms)
+__getgeom(::GI.PointTrait, geom) = error("PointTrait should not be reached")
+__getgeom(::Nothing, geoms) = geoms
diff --git a/src/methods/classify.jl b/src/methods/classify.jl
index 300770bd1..dbd501e4b 100644
--- a/src/methods/classify.jl
+++ b/src/methods/classify.jl
@@ -150,7 +150,11 @@ function classify!(A::AbstractRaster, pairs;
     end
     return rebuild(out; missingval=missingval)
 end
-function classify!(xs::RasterSeriesOrStack, pairs...; kw...)
+function classify!(xs::RasterStack, pairs...; kw...)
+    maplayers(x -> classify!(x, pairs...; kw...),  xs)
+    return xs
+end
+function classify!(xs::RasterSeries, pairs...; kw...)
     map(x -> classify!(x, pairs...; kw...),  xs)
     return xs
 end
diff --git a/src/methods/extract.jl b/src/methods/extract.jl
index 898bae166..ef41b04e7 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
@@ -88,11 +82,7 @@ function _extract(A::RasterStackOrArray, ::Nothing, data;
     names, skipmissing, geometrycolumn, kw...
 )
     geoms = _get_geometries(data, geometrycolumn)
-    T = if istrue(skipmissing)
-        _rowtype(A, nonmissingtype(eltype(geoms)); names, skipmissing, kw...)
-    else
-        _rowtype(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[]
     
@@ -101,7 +91,7 @@ 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))
+    rows = Vector{T}(undef, length(geoms))
         if istrue(skipmissing)
             j = 1
             for i in eachindex(geoms)
@@ -287,72 +277,7 @@ 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
-
-_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}
-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}
+    end |> T
 end
 
 @inline _skip_missing_rows(rows, ::Missing, names) = 
diff --git a/src/methods/mask.jl b/src/methods/mask.jl
index 232ed3b4e..1401a990c 100644
--- a/src/methods/mask.jl
+++ b/src/methods/mask.jl
@@ -175,7 +175,7 @@ function _mask!(x::RasterStackOrArray, geom; kw...)
 end
 # Array mask
 function _mask!(st::RasterStack, with::AbstractRaster; kw...)
-    map(A -> mask!(A; with, kw...), st)
+    maplayers(A -> mask!(A; with, kw...), st)
     return st
 end
 
@@ -431,7 +431,9 @@ end
 
 _false_to_missing(b::Bool) = (b ? true : missing)::Union{Missing,Bool}
 
-function _mask_multilayer(layers::Union{<:AbstractRasterStack,<:AbstractRasterSeries}, to; 
+_mask_multilayer(st::AbstractRasterStack, to; kw...) =
+    _mask_multilayer(layers(st), to; kw...)
+function _mask_multilayer(layers::Union{<:NamedTuple,<:AbstractRasterSeries}, to; 
     _dest_presentval,
     _dest_missingval, 
     missingval=nothing, 
diff --git a/src/methods/mosaic.jl b/src/methods/mosaic.jl
index 6e4b14ff8..8eaae3d61 100644
--- a/src/methods/mosaic.jl
+++ b/src/methods/mosaic.jl
@@ -163,11 +163,11 @@ nothing
 
 $EXPERIMENTAL
 """
-mosaic!(f::Function, x::RasterStackOrArray, regions::RasterStackOrArray...; kw...) =
-    mosaic!(f, x, regions; kw...)
-function mosaic!(f::Function, A::AbstractRaster{T}, regions;
-    missingval=Rasters.missingval(A),
-    atol=maybe_eps(T),
+mosaic!(f::Function, dest::RasterStackOrArray, regions::RasterStackOrArray...; kw...) =
+    _mosaic!(f, dest, regions; kw...)
+
+function _mosaic!(f::Function, A::AbstractRaster{T}, regions::Union{Tuple,AbstractArray};
+    missingval=missingval(A), atol=maybe_eps(T)
 ) where T
     isnokwornothing(missingval) && throw(ArgumentError("destination array must have a `missingval`"))
     _without_mapped_crs(A) do A1
@@ -200,13 +200,14 @@ function mosaic!(f::Function, A::AbstractRaster{T}, regions;
     end
     return A
 end
-function mosaic!(f::Function, st::AbstractRasterStack, regions::Tuple; kw...)
-    map(st, regions...) do A, r...
-        mosaic!(f, A, r; kw...)
+function _mosaic!(f::Function, st::AbstractRasterStack, regions::Union{Tuple,AbstractArray}; kw...)
+    map(values(st), map(values, regions)...) do A, r...
+        mosaic!(f, A, r...; kw...)
     end
+    return st
 end
 
-_mosaic(alldims::Tuple{<:DimTuple,Vararg{<:DimTuple}}) = map(_mosaic, alldims...)
+_mosaic(alldims::Tuple{<:DimTuple,Vararg{DimTuple}}) = map(_mosaic, alldims...)
 function _mosaic(dims::Dimension...)
     map(dims) do d
         DD.comparedims(first(dims), d; val=false, length=false, valtype=true)
diff --git a/src/methods/rasterize.jl b/src/methods/rasterize.jl
index 1e835c6e6..6c12a2eae 100644
--- a/src/methods/rasterize.jl
+++ b/src/methods/rasterize.jl
@@ -423,7 +423,7 @@ function rasterize(reducer::typeof(count), data; fill=nothing, init=nothing, kw.
 end
 # `mean` is sum ./ count. This is actually optimal with threading, 
 # as its means order is irrelevant so its threadsafe.
-function rasterize(reducer::typeof(DD.Statistics.mean), data; fill, kw...)
+function rasterize(reducer::typeof(Statistics.mean), data; fill, kw...)
     sums = rasterize(sum, data; kw..., fill)
     counts = rasterize(count, data; kw..., fill=nothing)
     rebuild(sums ./ counts; name=:mean)
@@ -830,6 +830,7 @@ end
 # We get 64 Bool values to a regular `Int` meaning this doesn't scale too
 # badly for large tables of geometries. 64k geometries and a 1000 * 1000
 # raster needs 1GB of memory just for the `BitArray`.
+# TODO combine these theyre nearly the same
 function _reduce_bitarray!(f, st::AbstractRasterStack, geoms, fill::NamedTuple, r::Rasterizer, allocs)
     (; lock, shape, boundary, verbose, progress, threaded) = r
     # Define mask dimensions, the same size as the spatial dims of x
@@ -839,7 +840,7 @@ function _reduce_bitarray!(f, st::AbstractRasterStack, geoms, fill::NamedTuple,
     # Use a generator over the array axis in case the iterator has no length
     geom_axis = axes(masks, Dim{:geometry}())
     fill = map(itr -> [v for (_, v) in zip(geom_axis, itr)], fill)
-    T = NamedTuple{keys(st),Tuple{map(eltype, st)...}}
+    T = eltype(st)
     range = axes(st, Y())
     _run(range, threaded, progress, "Reducing...") do y
         _reduce_bitarray_loop(f, st, T, fill, masks, y)
diff --git a/src/methods/reproject.jl b/src/methods/reproject.jl
index a4f417edb..7ce5aef46 100644
--- a/src/methods/reproject.jl
+++ b/src/methods/reproject.jl
@@ -74,7 +74,7 @@ function _reproject(source::Nothing, target::GeoFormat, dim::Union{XDim,YDim}, v
     reshape(reproject(source, target, dim, vec(vals)), size(vals))
 end
 function _reproject(source::GeoFormat, target::GeoFormat, dim::Union{XDim,YDim}, vals::AbstractVector) 
-    error("Rasters.jl requires backends to be loaded externally as of version 0.8. Run `using ArchGDAL` to use `reproject`")
+    throw_extension_error(reproject, "Proj", :RastersProjExt, (source, target, dim, vals))
 end
 
 # Guess the step for arrays
diff --git a/src/methods/slice_combine.jl b/src/methods/slice_combine.jl
index 724aeda46..98bb9f553 100644
--- a/src/methods/slice_combine.jl
+++ b/src/methods/slice_combine.jl
@@ -58,9 +58,9 @@ function _combine(ser::AbstractRasterSeries; lazy = isdisk(ser))
     ras1 = first(ser)
     alldims = (dims(ras1)..., dims(ser)...)
     ser_res = DD._insert_length_one_dims(ser, alldims)
-    data = DA.ConcatDiskArray(ser_res)
+    data = DA.ConcatDiskArray(map(parent, parent(ser_res)))
     data = lazy ? data : collect(data)
-    rebuild(ras1; data, dims = alldims, refdims = otherdims(dims(ras1, alldims)))
+    rebuild(ras1; data, dims=alldims, refdims=otherdims(dims(ras1, alldims)))
 end
 function _combine(ser::AbstractRasterSeries{<:AbstractRasterStack{K}}; kw...) where K
     r1 = first(ser)
@@ -69,4 +69,4 @@ function _combine(ser::AbstractRasterSeries{<:AbstractRasterStack{K}}; kw...) wh
     end |> NamedTuple{K}
     newdims = combinedims(new_layers...)
     DimensionalData.rebuild_from_arrays(r1, new_layers; refdims=otherdims(dims(r1), newdims))
-end
\ No newline at end of file
+end
diff --git a/src/methods/trim.jl b/src/methods/trim.jl
index 21830f8ef..588eb7dcf 100644
--- a/src/methods/trim.jl
+++ b/src/methods/trim.jl
@@ -131,5 +131,5 @@ end
 
 # Broadcast over the array and tracker to mark axis indices as being missing or not
 _update!(tr::AxisTrackers, A::AbstractRaster) = tr .= A .!== missingval(A)
-_update!(tr::AxisTrackers, st::AbstractRasterStack) = map(A -> tr .= A .!== missingval(A), st)
+_update!(tr::AxisTrackers, st::AbstractRasterStack) = maplayers(A -> tr .= A .!== missingval(A), st)
 
diff --git a/src/methods/zonal.jl b/src/methods/zonal.jl
index 509a47c33..c08b62cc5 100644
--- a/src/methods/zonal.jl
+++ b/src/methods/zonal.jl
@@ -84,7 +84,7 @@ _zonal(f, x::Raster, of::Extents.Extent; skipmissing=true) =
 function _zonal(f, x::RasterStack, ext::Extents.Extent; skipmissing=true)
     cropped = crop(x; to=ext, touches=true)
     length(cropped) > 0 || return missing
-    return map(cropped) do A
+    return maplayers(cropped) do A
         _maybe_skipmissing_call(f, A, skipmissing)
     end
 end
@@ -109,7 +109,7 @@ function _zonal(f, st::AbstractRasterStack, ::GI.AbstractGeometryTrait, geom;
     cropped = crop(st; to=geom, touches=true)
     length(cropped) > 0 || return map(_ -> missing, st)
     masked = mask(cropped; with=geom, kw...)
-    return map(masked) do A
+    return maplayers(masked) do A
         length(A) > 0 || return missing
         _maybe_skipmissing_call(f, A, skipmissing)
     end
diff --git a/src/read.jl b/src/read.jl
index eb6c7f622..af01f16ef 100644
--- a/src/read.jl
+++ b/src/read.jl
@@ -50,6 +50,10 @@ function Base.read!(src::AbstractRasterSeries, dst::AbstractRasterSeries)
     map(read!, src, dst)
     return dst
 end
+Base.read!(::AbstractRaster{<:Union{AbstractString,NamedTuple},1}, ::AbstractRasterSeries) =
+    error("Cant read Raster to RasterSeries")
+Base.read!(::AbstractRaster, ::AbstractRasterSeries) =
+    error("Cant read Raster to RasterSeries")
 
 # Filename methods
 function Base.read!(filename::AbstractString, dst::AbstractRaster)
diff --git a/src/series.jl b/src/series.jl
index 1b23caa6f..337eb8711 100644
--- a/src/series.jl
+++ b/src/series.jl
@@ -52,6 +52,8 @@ DD.modify(f, A::AbstractRasterSeries) = map(child -> modify(f, child), values(A)
     RasterSeries(paths::AbstractArray{<:AbstractString}, dims; child, duplicate_first, kw...)
     RasterSeries(path:::AbstractString, dims; ext, separator, child, duplicate_first, kw...)
 
+    RasterSeries(objects::AbstractBasicDimArray; kw...)
+
 Concrete implementation of [`AbstractRasterSeries`](@ref).
 
 A `RasterSeries` is an array of `Raster`s or `RasterStack`s, along some dimension(s).
@@ -111,11 +113,24 @@ When loading a series from a single `String` path:
 Others:
 - `refdims`: existing reference dimension/s, normally not required.
 """
-struct RasterSeries{T,N,D,R,A<:AbstractArray{T,N}} <: AbstractRasterSeries{T,N,D,A}
+struct RasterSeries{T<:Union{AbstractRaster,AbstractRasterStack},N,D<:Tuple,R<:Tuple,A<:AbstractArray{T,N}} <: AbstractRasterSeries{T,N,D,A}
     data::A
     dims::D
     refdims::R
+    # Needed for ambiguity with DD
+    function RasterSeries{T,N,D,R,A}(
+        data::AbstractArray, dims::Tuple, refdims::Tuple
+    ) where {T,N,D,R,A}
+        new{T,N,D,R,A}(data, dims, refdims)
+    end
 end
+function RasterSeries(
+    data::A, dims::D, refdims::R
+) where {A<:AbstractArray{T,N},D<:Tuple,R<:Tuple} where {T,N}
+    RasterSeries{T,N,D,R,A}(data, dims, refdims)
+end
+RasterSeries(data::DD.AbstractBasicDimArray; kw...) = 
+    RasterSeries(data, dims(data); kw...)
 function RasterSeries(data::AbstractArray{<:Union{AbstractRasterStack,AbstractRaster}}, dims;
     refdims=()
 )
@@ -205,6 +220,17 @@ end
 
 @inline function DD.rebuild(
     A::RasterSeries, data, dims::Tuple, refdims=(), name=nothing, metadata=nothing,
+)
+    # if `data` is not an AbstractArray of Raster or RasterStacks, return a Raster
+    Raster(data, dims; refdims, name, metadata)
+end
+@inline function DD.rebuild(
+    A::RasterSeries, 
+    data::AbstractArray{<:Union{AbstractRaster,AbstractRasterStack}}, 
+    dims::Tuple,
+    refdims=(),
+    name=nothing,
+    metadata=nothing,
 )
     RasterSeries(data, dims, refdims)
 end
@@ -212,7 +238,7 @@ end
     A::RasterSeries;
     data=parent(A), dims=dims(A), refdims=refdims(A), name=nothing, metadata=nothing,
 )
-    RasterSeries(data, dims, refdims)
+    rebuild(A, data, dims, refdims)
 end
 
 function Base.map(f, series::RasterSeries)
diff --git a/src/sources/commondatamodel.jl b/src/sources/commondatamodel.jl
index ded7d64dd..30c17cc09 100644
--- a/src/sources/commondatamodel.jl
+++ b/src/sources/commondatamodel.jl
@@ -166,7 +166,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
diff --git a/src/sources/grd.jl b/src/sources/grd.jl
index 2b53d7b46..ee1ba4a2d 100644
--- a/src/sources/grd.jl
+++ b/src/sources/grd.jl
@@ -79,8 +79,8 @@ function _dims(A::RasterDiskArray{GRDsource}, crs=nokw, mappedcrs=nokw)
     # Not fully implemented yet
     xy_metadata = _metadatadict(GRDsource())
 
-    xindex = LinRange(xbounds[1], xbounds[2] - xspan, xsize)
-    yindex = LinRange(ybounds[2] + yspan, ybounds[1], ysize)
+    xindex = range(; start=xbounds[1], stop=xbounds[2] - xspan, length=xsize)
+    yindex = range(; start=ybounds[2] + yspan, stop=ybounds[1], length=ysize)
 
     xlookup = Projected(xindex;
         order=GRD_X_ORDER,
diff --git a/src/sources/sources.jl b/src/sources/sources.jl
index 4c981478d..0c90cef65 100644
--- a/src/sources/sources.jl
+++ b/src/sources/sources.jl
@@ -58,7 +58,10 @@ end
 
 # error message to show when backend is not loaded
 function Base.showerror(io::IO, e::BackendException)
-    print(io, "`Rasters.jl` requires backends to be loaded externally as of version 0.8. Run `import $(e.backend)` to fix this error.")
+    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.")
 end
 
 # Get the source backend for a file extension, falling back to GDALsource
diff --git a/src/stack.jl b/src/stack.jl
index c16b445af..c80009ec5 100644
--- a/src/stack.jl
+++ b/src/stack.jl
@@ -31,8 +31,10 @@ abstract type AbstractRasterStack{K,T,N,L} <: AbstractDimStack{K,T,N,L} end
 missingval(stack::AbstractRasterStack) = getfield(stack, :missingval)
 missingval(s::AbstractRasterStack, name::Symbol) = _singlemissingval(missingval(s), name)
 
-filename(stack::AbstractRasterStack{<:Any,<:Any,<:Any,<:NamedTuple}) = map(s -> filename(s), stack)
-filename(stack::AbstractRasterStack{<:Any,<:Any,<:Any,<:Union{FileStack,OpenStack}}) = filename(parent(stack))
+filename(stack::AbstractRasterStack{<:Any,<:Any,<:Any,<:NamedTuple}) = 
+    map(s -> filename(s), layers(stack))
+filename(stack::AbstractRasterStack{<:Any,<:Any,<:Any,<:Union{FileStack,OpenStack}}) = 
+    filename(parent(stack))
 
 isdisk(st::AbstractRasterStack) = any(isdisk, layers(st))
 
@@ -92,12 +94,15 @@ function DD.rebuild(s::AbstractRasterStack;
 end
 
 function DD.rebuild_from_arrays(
-    s::AbstractRasterStack{<:Union{FileStack{<:Any,Keys},OpenStack{<:Any,Keys}}}, das::Tuple{Vararg{<:AbstractDimArray}}; kw...
+    s::AbstractRasterStack{<:Union{FileStack{<:Any,Keys},OpenStack{<:Any,Keys}}}, 
+    das::Tuple{Vararg{AbstractDimArray}}; 
+    kw...
 ) where Keys
     DD.rebuild_from_arrays(s, NamedTuple{Keys}(das); kw...)
 end
 function DD.rebuild_from_arrays(
-    s::AbstractRasterStack, das::NamedTuple{<:Any,<:Tuple{Vararg{AbstractDimArray}}};
+    s::AbstractRasterStack, 
+    das::NamedTuple{<:Any,<:Tuple{Vararg{AbstractDimArray}}};
     refdims=refdims(s),
     metadata=DD.metadata(s),
     dims=nothing,
@@ -121,12 +126,13 @@ end
 
 DD.name(s::AbstractRasterStack) = keys(s)
 Base.names(s::AbstractRasterStack) = keys(s)
-Base.copy(stack::AbstractRasterStack) = map(copy, stack)
+Base.copy(stack::AbstractRasterStack) = maplayers(copy, stack)
 
 #### Stack getindex ####
 # Different to DimensionalData as we construct a Raster
-Base.getindex(s::AbstractRasterStack, name::AbstractString) = s[Symbol(name)]
-function Base.getindex(s::AbstractRasterStack, name::Symbol)
+Base.@constprop :aggressive @propagate_inbounds Base.getindex(s::AbstractRasterStack, name::AbstractString) = 
+    s[Symbol(name)]
+Base.@constprop :aggressive @propagate_inbounds function Base.getindex(s::AbstractRasterStack, name::Symbol)
     data_ = parent(s)[name]
     dims_ = dims(s, DD.layerdims(s, name))
     metadata = DD.layermetadata(s, name)
@@ -215,7 +221,7 @@ function RasterStack(
     data::Union{FileStack,OpenStack,NamedTuple};
     dims::Tuple,
     refdims::Tuple=(),
-    layerdims,
+    layerdims::NamedTuple,
     metadata=nokw,
     layermetadata=nokw,
     missingval=nokw,
@@ -243,7 +249,7 @@ function RasterStack(
     return _postprocess_stack(st, dims; kw...)
 end
 # Convert Tuple/Array of array to NamedTuples using name/key
-function RasterStack(data::Tuple{Vararg{<:AbstractArray}}, dims::Tuple;
+function RasterStack(data::Tuple{Vararg{AbstractArray}}, dims::Tuple;
     name::Union{Tuple,AbstractArray,NamedTuple,Nothing}=nokw,
     kw...
 )
@@ -251,7 +257,7 @@ function RasterStack(data::Tuple{Vararg{<:AbstractArray}}, dims::Tuple;
     return RasterStack(NamedTuple{cleankeys(name)}(data), dims; kw...)
 end
 # Multi Raster stack from NamedTuple of AbstractArray
-function RasterStack(data::NamedTuple{<:Any,<:Tuple{Vararg{<:AbstractArray}}}, dims::Tuple;
+function RasterStack(data::NamedTuple{<:Any,<:Tuple{Vararg{AbstractArray}}}, dims::Tuple; 
     layerdims=nokw,
     kw...
 )
@@ -270,7 +276,7 @@ end
 # Multi Raster stack from AbstractDimArray splat
 RasterStack(layers::AbstractDimArray...; kw...) = RasterStack(layers; kw...)
 # Multi Raster stack from tuple with `name` keyword
-function RasterStack(layers::Tuple{Vararg{<:AbstractDimArray}};
+function RasterStack(layers::Tuple{Vararg{AbstractDimArray}};
     name=map(name, layers),
     kw...
 )
@@ -278,7 +284,7 @@ function RasterStack(layers::Tuple{Vararg{<:AbstractDimArray}};
 end
 # Multi RasterStack from NamedTuple
 # This method is called after most other RasterStack methods.
-function RasterStack(layers::NamedTuple{K,<:Tuple{Vararg{<:AbstractDimArray}}};
+function RasterStack(layers::NamedTuple{K,<:Tuple{Vararg{AbstractDimArray}}};
     resize::Union{Function,NoKW}=nokw,
     _layers=resize isa NoKW ? layers : resize(layers),
     dims::Tuple=DD.combinedims(_layers...),
@@ -428,10 +434,10 @@ function RasterStack(filename::AbstractString;
                 source, name, lazy, group, missingval, scaled, coerce, kw...
             )
             # Maybe split the stack into separate arrays to remove extra dims.
-            if !isnokw(name)
-                map(identity, l_st)
-            else
+            if isnokw(name)
                 l_st
+            else
+                maplayers(identity, l_st)
             end
         else
             # With bands actings as layers
@@ -446,12 +452,14 @@ function RasterStack(filename::AbstractString;
     return _maybe_drop_single_band(st, dropband, lazy)
 end
 
+# TODO test this properly
 function DD.modify(f, s::AbstractRasterStack{<:FileStack{<:Any,K}}) where K
-    open(s) do o
+    data = open(s) do o
         map(K) do k
-            Array(parent(ost)[k])
+            f(parent(ost)[k])
         end
     end
+    rebuild(s; data)
 end
 
 # Open a single file stack
diff --git a/src/table_ops.jl b/src/table_ops.jl
index 44213d128..84f9a1ae5 100644
--- a/src/table_ops.jl
+++ b/src/table_ops.jl
@@ -8,7 +8,7 @@ function _auto_dim_columns(table, dims::Tuple)
 end
 
 _not_a_dimcol(table, dimcols::DimTuple) = _not_a_dimcol(table, map(name, dimcols))
-_not_a_dimcol(table, dimcols::Tuple{Vararg{<:Pair}}) = _not_a_dimcol(table, map(last, dimcols))
+_not_a_dimcol(table, dimcols::Tuple{Vararg{Pair}}) = _not_a_dimcol(table, map(last, dimcols))
 _not_a_dimcol(table, dimcols::Tuple{}) = ()
 function _not_a_dimcol(table, dimcols::Tuple{Vararg{Symbol}})
     dimcols = (dimcols..., :Band) # TODO deal with how annoying `Band` is 
diff --git a/src/utils.jl b/src/utils.jl
index 29c613b0b..d1a50767f 100644
--- a/src/utils.jl
+++ b/src/utils.jl
@@ -24,6 +24,16 @@
 # const URLREGEX = r"^[a-zA-Z][a-zA-Z\d+\-.]*:"
 
 # _isurl(str::AbstractString) = !occursin(WINDOWSREGEX, str) && occursin(URLREGEX, str)
+function _maybe_use_type_missingval(A::AbstractRaster{T}, source::Source, missingval=nokw) where T
+    if ismissing(Rasters.missingval(A))
+        newmissingval = missingval isa NoKW ? _type_missingval(Missings.nonmissingtype(T)) : missingval
+        A1 = replace_missing(A, newmissingval)
+        @warn "`missing` cant be written with $(SOURCE2SYMBOL[source]), missingval for `$(eltype(A1))` of `$newmissingval` used instead"
+        return A1
+    else
+        return A
+    end
+end
 
 # cleankeys(name) = (_cleankey(name),)
 # function cleankeys(keys::Union{NamedTuple,Tuple,AbstractArray})
@@ -482,3 +492,73 @@ _filenotfound_error(filename) = throw(ArgumentError("file \"$filename\" not foun
 _size_and_res_error() = throw(ArgumentError("Both `size` and `res` keywords are passed, but only one can be used"))
 
 _no_crs_error() = throw(ArgumentError("The provided object does not have a CRS. Use `setcrs` to set one."))
+
+
+# Rows for extraction etc
+
+# _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
+
+# 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)
+    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/array.jl b/test/array.jl
index e7c4ab231..fbd52580d 100644
--- a/test/array.jl
+++ b/test/array.jl
@@ -23,6 +23,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
diff --git a/test/cellarea.jl b/test/cellarea.jl
new file mode 100644
index 000000000..8f4c66ac3
--- /dev/null
+++ b/test/cellarea.jl
@@ -0,0 +1,73 @@
+using Rasters, DimensionalData, Rasters.Lookups, Proj
+using Test
+using DimensionalData: @dim, YDim
+include(joinpath(dirname(pathof(Rasters)), "../test/test_utils.jl"))
+
+@testset "cellarea" begin
+    x = X(Projected(90.0:0.1:99.9; sampling=Intervals(Start()), order = ForwardOrdered(), span = Regular(0.1), crs=EPSG(4326), dim = X()))
+    y = Y(Projected(0.0:0.1:89.9; sampling=Intervals(Start()), order = ForwardOrdered(), span = Regular(0.1), crs=EPSG(4326), dim = Y()))
+    y_rev = Y(Projected(89.9:-0.1:0; sampling=Intervals(Start()), order = ReverseOrdered(), span = Regular(-0.1), crs=EPSG(4326), dim = Y()))
+    dimz = (x,y)
+    ras = ones(dimz)
+
+    dimz_25832 = X(Projected(0.0:100:10000.0; sampling=Intervals(Start()), order = ForwardOrdered(), span = Regular(100), crs=EPSG(25832))),
+       Y(Projected(0.0:100:10000.0; sampling=Intervals(Start()), order = ForwardOrdered(), span = Regular(100), crs=EPSG(25832)))
+    # to test that this is still stable on as 5x5m grid
+    smallspan_25832 = X(Projected(0.0:5:100.0; sampling=Intervals(Start()), order = ForwardOrdered(), span = Regular(5), crs=EPSG(25832))),
+       Y(Projected(0.0:5:100.0; sampling=Intervals(Start()), order = ForwardOrdered(), span = Regular(5), crs=EPSG(25832)))
+
+    cs = cellarea(dimz)
+    cs2 = cellarea(dimz_25832)
+    cs3 = cellarea((x, y_rev))
+    cs4 = cellarea(smallspan_25832)
+
+    # check cellsize still works
+    csold = @test_warn "deprecated" cellsize(dimz)
+    @test all(csold .≈ (cs ./ 1e6))
+
+    # Check the output is a raster 
+    @test cs isa Raster
+    @test cs2 isa Raster
+    # Test that the total area matches the expected area (1/72th of the Earth surface)
+    @test sum(cs) ≈ sum(cs3)
+    @test sum(cs) ≈ 510.1e12/72 rtol = 0.01
+    # Test all areas are about 10,000 m2 for the 100x100m grid
+    @test maximum(cs2) ≈ 1e4 rtol = 0.01
+    @test minimum(cs2) ≈ 1e4 rtol = 0.01
+    # Test all areas are about 25 m2 on the 5x5m grid
+    @test maximum(cs4) ≈ 25 rtol = 0.01
+    @test minimum(cs4) ≈ 25 rtol = 0.01
+    # test passing in a raster or dims gives the same result
+    cs_ras = cellarea(ras)
+    @test cs == cs_ras
+
+    ## test with area_in_crs true
+    cs_planar = cellarea(Planar(), dimz)
+    cs_planar2 = cellarea(Planar(), dimz_25832)
+    cs_planar3 = cellarea(Planar(), (x, y_rev))
+    @test all(≈(0.01), cs_planar)
+    @test all(≈(10_000), cs_planar2)
+    @test all(≈(0.01), cs_planar3)
+
+    # test a YDim other than Y is handled correctly
+    @dim Lat YDim "latitude"
+    lat = Lat(Projected(0.0:0.1:89.9; sampling=Intervals(Start()), order = ForwardOrdered(), span = Regular(0.1), crs=EPSG(4326)))
+    cs_lat = cellarea((x, lat))
+    @test parent(cs_lat) == parent(cs)
+    @test dims(cs_lat) == (x, lat)
+
+    # test using mapped dims to calculate cellsize
+    x_3857 = X(Projected(0:1e5:1e7; sampling=Intervals(Start()), order = ForwardOrdered(), span = Regular(1e5), crs=EPSG(3857), dim = X()))
+    y_3857 = Y(Projected(0:1e5:1e7; sampling=Intervals(Start()), order = ForwardOrdered(), span = Regular(1e5), crs=EPSG(3857), dim = Y()))
+    dimz_3857 = (x_3857, y_3857)
+    mappeddimz = setmappedcrs(dimz_3857, EPSG(4326))
+    @test all(isapprox.(cellarea((dimz_3857)), cellarea((mappeddimz)); rtol = 0.001))
+    @test parent(cellarea((dimz_3857))) != parent(cellarea((mappeddimz))) # to make sure it actually used the mapped dims
+    # test point sampling throws an error
+    pointsy = set(y, Points())
+    @test_throws ArgumentError cellarea((x, pointsy))
+
+    # test missing crs throws an error
+    nocrsdimz = setcrs(dimz, nothing)
+    @test_throws ArgumentError cellarea(nocrsdimz)
+end
diff --git a/test/cellsize.jl b/test/cellsize.jl
deleted file mode 100644
index aa74202bb..000000000
--- a/test/cellsize.jl
+++ /dev/null
@@ -1,45 +0,0 @@
-using Rasters, DimensionalData, Rasters.Lookups, ArchGDAL
-using Test
-using DimensionalData: @dim, YDim
-include(joinpath(dirname(pathof(Rasters)), "../test/test_utils.jl"))
-
-@testset "cellsize" begin
-    x = X(Projected(90.0:0.1:99.9; sampling=Intervals(Start()), order = ForwardOrdered(), span = Regular(0.1), crs=EPSG(4326)))
-    y = Y(Projected(0.0:0.1:89.9; sampling=Intervals(Start()), order = ForwardOrdered(), span = Regular(0.1), crs=EPSG(4326)))
-    dimz = (x,y)
-
-    dimz_25832 = X(Projected(0.0:100:10000.0; sampling=Intervals(Start()), order = ForwardOrdered(), span = Regular(100), crs=EPSG(25832))),
-       Y(Projected(0.0:100:10000.0; sampling=Intervals(Start()), order = ForwardOrdered(), span = Regular(100), crs=EPSG(25832)))
-    ras = ones(dimz)
-
-    cs = cellsize(dimz)
-    cs2 = cellsize(dimz_25832)
-
-    # Check the output is a raster 
-    @test cs isa Raster
-    @test cs2 isa Raster
-    # Test that the total area matches the expected area (1/72th of the Earth surface)
-    @test sum(cs) ≈ 510.1e6/72 rtol = 0.01
-    # Test all areas are about 0.01 km2
-    @test maximum(cs2) ≈ 0.01 rtol = 0.01
-    @test minimum(cs2) ≈ 0.01 rtol = 0.01
-    # test passing in a raster or dims gives the same result
-    cs_ras = cellsize(ras)
-    @test cs == cs_ras
-
-    # test a YDim other than Y is handled correctly
-    @dim Lat YDim "latitude"
-    lat = Lat(Projected(0.0:0.1:89.9; sampling=Intervals(Start()), order = ForwardOrdered(), span = Regular(0.1), crs=EPSG(4326)))
-    cs_lat = cellsize((x, lat))
-    @test parent(cs_lat) == parent(cs)
-    @test dims(cs_lat) == (x, lat)
-
-    # test point sampling throws an error
-    pointsy = set(y, Points())
-    @test_throws ArgumentError cellsize((x, pointsy))
-
-    # test missing crs throws an error
-    nocrsdimz = setcrs(dimz, nothing)
-    @test_throws ArgumentError cellsize(nocrsdimz)
-
-end
\ No newline at end of file
diff --git a/test/methods.jl b/test/methods.jl
index 71ee990fa..3abb0a5d1 100644
--- a/test/methods.jl
+++ b/test/methods.jl
@@ -60,7 +60,7 @@ gaMi = replace_missing(ga)
     read(dNaN)
     @test all(isequal.(dNaN, [NaN32 7.0f0; 2.0f0 NaN32]))
     stNaN = replace_missing(st, NaN32; filename="teststack.tif")
-    @test all(map(stNaN[Band(1)], (a=[NaN32 7.0f0; 2.0f0 NaN32], b=[1.0 0.4; 2.0 NaN])) do x, y
+    @test all(maplayers(stNaN[Band(1)], (a=[NaN32 7.0f0; 2.0f0 NaN32], b=[1.0 0.4; 2.0 NaN])) do x, y
         all(x .=== y)
     end)
     rm("teststack_a.tif")
@@ -221,6 +221,12 @@ end
             @test sum(skipmissing(mask(polytemplate; with=polygon, boundary=:touches, invert=true))) == prod(size(polytemplate)) - 21 * 21
         end
     end
+
+    @testset "geometry encompassing raster" begin
+        geom = GeoInterface.Polygon([GeoInterface.LinearRing([(0.0, 0.0), (0.0, 10.0), (10.0, 10.0), (10.0, 0.0), (0.0, 0.0)])])
+        raster = Raster(ones(X(1:0.1:2), Y(1:0.1:2)), missingval=false)
+        @test sum(mask(raster; with=geom)) == sum(raster)
+    end
 end
 
 @testset "mask" begin
@@ -304,7 +310,7 @@ end
         zonal(sum, st; of=(geometry=polygon, x=:a, y=:b)) ==
         zonal(sum, st; of=[(geometry=polygon, x=:a, y=:b)])[1] ==
         zonal(sum, st; of=[(geometry=polygon, x=:a, y=:b)])[1] ==
-        map(sum ∘ skipmissing, mask(st; with=polygon))  
+        maplayers(sum ∘ skipmissing, mask(st; with=polygon))  
     @test zonal(sum, st; of=st) == 
         zonal(sum, st; of=dims(st)) == 
         zonal(sum, st; of=Extents.extent(st)) == 
@@ -401,7 +407,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)
@@ -414,7 +420,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)
@@ -449,7 +455,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)
@@ -558,11 +564,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)
@@ -576,56 +587,56 @@ 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))
         filename = tempname() * ".tif"
         extended_d = extend(cropped; to=ga, filename)
-        @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
@@ -674,3 +685,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
diff --git a/test/rasterize.jl b/test/rasterize.jl
index c74a30b6f..919583b87 100644
--- a/test/rasterize.jl
+++ b/test/rasterize.jl
@@ -198,7 +198,7 @@ end
             @test eltype(rst) == @NamedTuple{fill1::Union{Missing,Int64}, fill2::Union{Missing,Float32}}
             @test keys(rst) == (:fill1, :fill2)
             @test dims(rst) == dims(A)
-            @test map(sum ∘ skipmissing, rst) === (fill1=15, fill2=30.0f0)
+            @test maplayers(sum ∘ skipmissing, rst) === (fill1=15, fill2=30.0f0)
         end
         @testset "Single value fill makes an array (ignoring table vals)" begin
             ra = rasterize(sum, data; to=A, fill=0x03, missingval=0x00)
@@ -216,7 +216,7 @@ end
             @test parent(rst.val2) isa Array{Union{Missing,Float32},2}
             @test keys(rst) == (:val1, :val2)
             @test dims(rst) == dims(A)
-            @test map(sum ∘ skipmissing, rst) === (val1=4, val2=8.0f0)
+            @test maplayers(sum ∘ skipmissing, rst) === (val1=4, val2=8.0f0)
         end
         @testset "Symbol fill makes an array" begin
             ra = rasterize(feature; to=A, fill=:val1)
@@ -235,7 +235,7 @@ end
                 @test parent(rst.val1) isa Array{Union{Missing,Int},2}
                 @test parent(rst.val2) isa Array{Union{Missing,Float32},2}
                 @test keys(rst) == (:val1, :val2)
-                @test map(sum ∘ skipmissing, rst) === (val1=15, val2=30.0f0)
+                @test maplayers(sum ∘ skipmissing, rst) === (val1=15, val2=30.0f0)
                 @test_throws ArgumentError rasterize(data; to=A, fill=(:val1, :not_a_column), threaded)
             end
             @testset "Symbol fill makes an array" begin
diff --git a/test/reproject.jl b/test/reproject.jl
index e471361d8..3bb337b6c 100644
--- a/test/reproject.jl
+++ b/test/reproject.jl
@@ -1,13 +1,13 @@
-using Rasters, Test, ArchGDAL
+using Rasters, Test, Proj
 using Rasters.Lookups, Rasters.Dimensions
 using Rasters: reproject, convertlookup
 
 @testset "reproject" begin
-    cea = ProjString("+proj=cea +lon_0=0 +lat_ts=30 +x_0=0 +y_0=0 +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0")
-    wktcea = convert(WellKnownText, cea)
-    projcea = convert(ProjString, cea)
-    wkt4326 = convert(WellKnownText, EPSG(4326))
-    proj4326 = convert(ProjString, EPSG(4326))
+    cea = ProjString("+proj=cea +lon_0=0 +lat_ts=30 +x_0=0 +y_0=0 +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0 +type=crs")
+    wktcea = convert(WellKnownText, Proj.CRS(cea))
+    projcea = convert(ProjString, Proj.CRS(cea))
+    wkt4326 = convert(WellKnownText, Proj.CRS(EPSG(4326)))
+    proj4326 = convert(ProjString, Proj.CRS(EPSG(4326)))
 
 
     @test reproject(proj4326, cea, Y(), -30.0) ≈ -3.658789324855012e6
@@ -36,8 +36,9 @@ using Rasters: reproject, convertlookup
     reprojected_raster = reproject(raster; crs=projcea)
     # Dims have changed
     @test dims(reprojected_raster) == (x1, y1)
-    # Raster has not changed
-    @test reprojected_raster == raster
+    # Parent data has not changed
+    # NOTE: this test only works because both rasters are in memory.
+    @test parent(reprojected_raster) == parent(raster)
 
     @test_throws ArgumentError reproject(cea, EPSG(32618), Y(), [-3.658789324855012e6])
     @test_throws ArgumentError reproject(cea, EPSG(32618), X(), [-3.658789324855012e6])
@@ -48,8 +49,8 @@ using Rasters: reproject, convertlookup
 end
 
 @testset "convertlookup" begin
-    projcea = ProjString("+proj=cea")
-    proj4326 = convert(ProjString, EPSG(4326))
+    projcea = ProjString("+proj=cea +type=crs")
+    proj4326 = convert(ProjString, Proj.CRS(EPSG(4326)))
 
     lonstart, lonend = 0.5, 179.5
     cealonstart, cealonend = reproject(proj4326, projcea, X(), [lonstart, lonend])
diff --git a/test/resample.jl b/test/resample.jl
index d1234d349..c24179906 100644
--- a/test/resample.jl
+++ b/test/resample.jl
@@ -1,4 +1,5 @@
-using Rasters, ArchGDAL, GeoInterface, Extents
+using Rasters, ArchGDAL, GeoInterface
+using Extents
 using Test
 using Rasters.Lookups
 import DimensionalData: @dim, YDim
@@ -115,31 +116,6 @@ include(joinpath(dirname(pathof(Rasters)), "../test/test_utils.jl"))
         @test eltype(r1) == UInt8
     end
 
-    @testset "resample to the same size or resolution leaves raster unchanged" begin
-        res = 2
-        ys = LinRange((90 - res / 2):-res:(-90 + res / 2))
-        xs = LinRange((-180 + res / 2):res:(180 - res / 2))
-
-        point_dims = Y(ys), X(xs)
-        interval_dims = Y(ys; sampling=Intervals(Center())), X(xs; sampling=Intervals(Center()))
-        rev_y_point_dims = Y(reverse(ys)), X(xs)
-        rev_y_interval_dims = reverse(interval_dims[1]), interval_dims[2]
-        rev_x_point_dims = Y(ys), X(reverse(xs))
-        rev_x_interval_dims = interval_dims[1], reverse(interval_dims[2])
-        test_dims = (point_dims, interval_dims, rev_x_point_dims, rev_x_interval_dims, rev_y_point_dims, rev_y_interval_dims)
-        ds_fwd = point_dims; f = identity
-        ds_fwd = point_dims; f = reverse
-
-        for ds_fwd in test_dims, f in (identity, reverse)
-            ds = f(ds_fwd)
-            raster = Raster(rand(ds), crs=EPSG(4326), missingval=NaN)
-            resampled_res = resample(raster; res)
-            resampled_size = resample(raster; size=size(raster), method=:near)
-            @test resampled_res == resampled_size == raster
-            @test dims(resampled_res) == dims(resampled_size) == dims(raster)
-        end
-    end
-
     @testset "dimensions matcha after resampling with only `to`" begin
         # some weird dimensions
         to = (
@@ -160,7 +136,7 @@ include(joinpath(dirname(pathof(Rasters)), "../test/test_utils.jl"))
     maskingval = Rasters.nokw
     for maskingval in (nothing, missing, Rasters.nokw)
         # Resample cea.tif using resample
-        cea = Raster(raster_path; missingval=0x00=>maskingval, name=:cea)
+        cea = Raster(raster_path; missingval=0x00 => maskingval, name=:cea)
         raster_output = resample(cea; res=output_res, crs=output_crs, method, missingval=0x00 => maskingval)
         disk_output = resample(cea; res=output_res, crs=output_crs, method, missingval=0x00 => maskingval, filename="resample.tif")
 
@@ -183,4 +159,35 @@ include(joinpath(dirname(pathof(Rasters)), "../test/test_utils.jl"))
 
         rm("resample.tif")
     end
+
+    @testset "resample to the same size or resolution leaves raster unchanged" begin
+        res = 2
+        ys = range(; start=(90 - res / 2), step=-res, stop=(-90 + res / 2))
+        xs = range(; start=(-180 + res / 2), step=res, stop=(180 - res / 2))
+
+        point_dims = Y(ys), X(xs)
+        interval_dims = Y(ys; sampling=Intervals(Center())), X(xs; sampling=Intervals(Center()))
+        rev_y_point_dims = Y(reverse(ys)), X(xs)
+        rev_y_interval_dims = reverse(interval_dims[1]), interval_dims[2]
+        rev_x_point_dims = Y(ys), X(reverse(xs))
+        rev_x_interval_dims = interval_dims[1], reverse(interval_dims[2])
+        test_dims = (point_dims, interval_dims, rev_x_point_dims, rev_x_interval_dims, rev_y_point_dims, rev_y_interval_dims)
+        ds_fwd = point_dims; f = identity
+        ds_fwd = point_dims; f = reverse
+
+        ds = ds_fwd
+        raster = Raster(rand(ds), crs=EPSG(4326), missingval=NaN)
+        resampled_res = resample(raster; res)
+        resampled_size = resample(raster; size=size(raster), method=:near)
+        @test resampled_res == resampled_size == raster
+        @test dims(resampled_res) == dims(resampled_size) == dims(raster)
+
+        ds = reverse(ds_fwd)
+        raster = Raster(rand(ds), crs=EPSG(4326), missingval=NaN)
+        resampled_res = resample(raster; res)
+        resampled_size = resample(raster; size=size(raster), method=:near)
+        @test resampled_res == resampled_size == raster
+        @test dims(resampled_res) == dims(resampled_size) == dims(raster)
+    end
+
 end
diff --git a/test/runtests.jl b/test/runtests.jl
index 34b8a6354..55f13eb5e 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -1,48 +1,43 @@
 using Rasters, Test, Aqua, SafeTestsets
 
-@testset "Rasters" begin
-    @testset "Aqua" begin
-        # Aqua.test_ambiguities([Rasters, Base, Core])
-        Aqua.test_unbound_args(Rasters)
-        Aqua.test_stale_deps(Rasters)
-        Aqua.test_undefined_exports(Rasters)
-        Aqua.test_project_extras(Rasters)
-        # Aqua.test_deps_compat(Rasters) # This breaks GDAL downstream tests
-        # Aqua.test_project_toml_formatting(Rasters) # This seems to change between versions for extensions
-    end
-
-    @time @safetestset "extensions" begin include("extensions.jl") end
-    @time @safetestset "methods" begin include("methods.jl") end
-    @time @safetestset "array" begin include("array.jl") end
-    @time @safetestset "stack" begin include("stack.jl") end
-    @time @safetestset "series" begin include("series.jl") end
-    @time @safetestset "create" begin include("create.jl") end
-    @time @safetestset "utils" begin include("utils.jl") end
-    @time @safetestset "set" begin include("set.jl") end
-    @time @safetestset "aggregate" begin include("aggregate.jl") end
-    @time @safetestset "rasterize" begin include("rasterize.jl") 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 "cellsize" begin include("cellsize.jl") end
-
-    # CommondataModel sources
-    @time @safetestset "ncdatasets" begin include("sources/ncdatasets.jl") end
-    if !Sys.iswindows()
-        # GRIBDatasets doesn't work on Windows for now
-        @time @safetestset "gribdatasets" begin include("sources/gribdatasets.jl") end
-    end
+@testset "Aqua" begin
+    Aqua.test_ambiguities(Rasters)
+    Aqua.test_unbound_args(Rasters)
+    Aqua.test_stale_deps(Rasters)
+    Aqua.test_undefined_exports(Rasters)
+    Aqua.test_project_extras(Rasters)
+    # Aqua.test_deps_compat(Rasters) # This breaks GDAL downstream tests
+end
 
-    # Only test SMAP locally for now, also RasterDataSources because CI downloads keep breaking
-    if !haskey(ENV, "CI")
-        @time @safetestset "rasterdatasources" begin include("sources/rasterdatasources.jl") end
-    end
+@time @safetestset "extensions" begin include("extensions.jl") end
+@time @safetestset "methods" begin include("methods.jl") end
+@time @safetestset "array" begin include("array.jl") end
+@time @safetestset "stack" begin include("stack.jl") end
+@time @safetestset "series" begin include("series.jl") end
+@time @safetestset "create" begin include("create.jl") end
+@time @safetestset "utils" begin include("utils.jl") end
+@time @safetestset "set" begin include("set.jl") end
+@time @safetestset "aggregate" begin include("aggregate.jl") end
+@time @safetestset "rasterize" begin include("rasterize.jl") 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 "cellarea" begin include("cellarea.jl") end
 
-    if !Sys.iswindows()
-        # GDAL Environment vars need to be set manually for windows, so skip for now
-        @time @safetestset "gdal" begin include("sources/gdal.jl") end
-        @time @safetestset "grd" begin include("sources/grd.jl") end
-    end
-    @time @safetestset "plot recipes" begin include("plotrecipes.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
+    # GDAL Environment vars need to be set manually for windows, so skip for now
+    @time @safetestset "gdal" begin include("sources/gdal.jl") end
+    @time @safetestset "grd" begin include("sources/grd.jl") end
+end
+# Only test SMAP locally for now, also RasterDataSources because CI downloads keep breaking
+if !haskey(ENV, "CI")
+    @time @safetestset "rasterdatasources" begin include("sources/rasterdatasources.jl") end
 end
+@time @safetestset "plot recipes" begin include("plotrecipes.jl") end
+@time @safetestset "resample" begin include("resample.jl") end
diff --git a/test/series.jl b/test/series.jl
index 45647c9dc..e5ebb5ed8 100644
--- a/test/series.jl
+++ b/test/series.jl
@@ -19,6 +19,7 @@ stack1 = RasterStack(r1, r2; name=(:r1, :r2))
 stack2 = RasterStack(r1a, r2a; name=(:r1, :r2))
 dates = [DateTime(2017), DateTime(2018)]
 ser = RasterSeries([stack1, stack2], Ti(dates))
+@test RasterSeries(DimArray([stack1, stack2],Ti(dates))) == ser
 @test issorted(dates)
 
 @testset "getindex returns the currect types" begin
@@ -76,6 +77,7 @@ end
     r1 = Raster(ones(4, 5, 10), (X(), Y(), Ti(10:10:100))) .* reshape(1:10, (1, 1, 10))
     r2 = r1 .* 2
     ser = slice(r1, Ti)
+
     @test size(ser) == (10,)
     combined = Rasters.combine(ser, Ti())
     @test combined == r1
diff --git a/test/sources/commondatamodel.jl b/test/sources/commondatamodel.jl
index 451d18b31..6e2460092 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/gdal.jl b/test/sources/gdal.jl
index f3986e659..50271eabd 100644
--- a/test/sources/gdal.jl
+++ b/test/sources/gdal.jl
@@ -450,7 +450,7 @@ gdalpath = maybedownload(url)
             @test (@allocations write(filename2, gdalarray; force=true)) < 1e4
             saved = Raster(filename2; crs=crs(gdalarray), mappedcrs=crs(gdalarray))
             @test size(saved) == size(gdalarray)
-            @test saved ≈ gdalarray
+            @test parent(saved) ≈ parent(gdalarray)
             clat, clon = DimensionalData.shiftlocus.(Ref(Center()), dims(gdalarray, (Y, X)))
             @test index(clat) ≈ index(saved, Y)
             @test index(clon) ≈ index(saved, X)
@@ -580,7 +580,7 @@ gdalpath = maybedownload(url)
         @test order(dims(rast)) == (ForwardOrdered(), ForwardOrdered())
         @test span(rast) == (Regular(1.0), Regular(1.0))
         @test sampling(rast) == (Intervals(Start()), Intervals(Start()))
-        @test index(rast) == (LinRange(0.0, 239.0, 240), LinRange(0.0, 179.0, 180))
+        @test index(rast) == (range(; start=0.0, stop=239.0, length=240), range(start=0.0, stop=179.0, length=180))
     end
 
 end
@@ -653,20 +653,20 @@ end
 
     @testset "methods" begin
         @testset "mean" begin
-            means = map(A -> mean(parent(A); dims=2), gdalstack)
-            @test map((a, b) -> all(a .== b), mean(gdalstack; dims=Y), means) |> all
+            means = maplayers(A -> mean(parent(A); dims=2), gdalstack)
+            @test maplayers((a, b) -> all(a .== b), mean(gdalstack; dims=Y), means) |> all
         end
         @testset "trim, crop, extend" begin
             mv = zero(eltype(gdalstack[:a]))
             st = read(replace_missing(gdalstack, mv))
-            st = map(A -> (view(A, X(1:100)) .= mv; A), st)
+            st = maplayers(A -> (view(A, X(1:100)) .= mv; A), st)
             trimmed = trim(st)
             @test size(trimmed) == (414, 514)
             cropped = crop(st; to=trimmed)
             @test size(cropped) == (414, 514)
-            @test map((c, t) -> all(collect(c .=== t)), cropped, trimmed) |> all
+            @test all(maplayers((c, t) -> all(collect(c .=== t)), cropped, trimmed))
             extended = extend(read(cropped); to=st)
-            @test all(map((s, e) -> all(s .=== e), st, extended))
+            @test all(maplayers((s, e) -> all(s .=== e), st, extended))
         end
         @testset "mask and mask!" begin
             st = read(gdalstack)
@@ -729,7 +729,7 @@ end
             # Test forcing
             write(filename, gdalstack; force=true);
             saved = RasterStack(filename);
-            @test all(read(saved[:a]) .== geoA)
+            @test all(parent(read(saved[:a])) .== parent(geoA))
             rm(filename)
         end
 
@@ -888,7 +888,7 @@ end
         read!([(a=gdalpath, b=gdalpath), (a=gdalpath, b=gdalpath)], ser2)
         read!(ser1, ser3)
         @test map(ser1, ser2, ser3) do st1, st2, st3
-            map(st1, st2, st3) do A1, A2, A3
+            maplayers(st1, st2, st3) do A1, A2, A3
                 (A2 .=== A2 .=== A3) |> all
             end |> all
         end |> all
diff --git a/test/sources/grd.jl b/test/sources/grd.jl
index dedf915f6..214b2e6a4 100644
--- a/test/sources/grd.jl
+++ b/test/sources/grd.jl
@@ -228,8 +228,8 @@ grdpath = stem * ".gri"
             # 1 band is added again on save
             @test size(saved) == size(grdarray[Band(1)])
             @test parent(saved) == parent(grdarray[Band(1)])
-            write(filename2, grdarray[Band(1)]; force=true, verbose=false)
-            @test (@allocations write(filename2, grdarray[Band(1)]; force=true, verbose=false)) < 1e3
+            write(filename2, grdarray; force=true)
+            @test (@allocations write(filename2, grdarray; force=true, verbose=false)) < 3e3
         end
 
         @testset "3d with subset" begin
@@ -255,7 +255,8 @@ grdpath = stem * ".gri"
             @test all(parent(saved) .=== parent(geoA))
             @test saved isa typeof(geoA)
             @test parent(saved) == parent(geoA)
-            @test (@allocations write(filename, GRDsource(), geoA; force = true)) < 1e3
+            write(filename, GRDsource(), geoA; force = true)
+            @test (@allocations write(filename, GRDsource(), geoA; force = true)) < 3e3
         end
 
         @testset "to netcdf" begin
@@ -264,20 +265,21 @@ grdpath = stem * ".gri"
             write(filename2, grdarray[Band(1)]; force = true)
             saved = Raster(filename2; crs=crs(grdarray))
             @test size(saved) == size(grdarray[Band(1)])
-            @test all(replace_missing(saved, missingval(grdarray)) .≈ grdarray[Band(1)])
+            @test all(parent(replace_missing(saved, missingval(grdarray))) .≈ parent(grdarray[Band(1)]))
             @test index(saved, X) ≈ index(grdarray, X) .+ 0.5
             @test index(saved, Y) ≈ index(grdarray, Y) .+ 0.5
             @test bounds(saved, Y) == bounds(grdarray, Y)
             @test bounds(saved, X) == bounds(grdarray, X)
-            @test (@allocations write(filename2, grdarray[Band(1)]; force = true)) < 1e3
+            write(filename2, grdarray; force = true)
+            @test (@allocations write(filename2, grdarray; force = true)) < 3e3
         end
 
         @testset "to gdal" begin
             # No Band
             gdalfilename = tempname() * ".tif"
             write(gdalfilename, GDALsource(), grdarray[Band(1)]; force = true)
-            @test (@allocations write(gdalfilename, GDALsource(), grdarray[Band(1)]; force = true)) < 1e4
-            gdalarray = Raster(gdalfilename; maskingval=nothing)
+            @test (@allocations write(gdalfilename, GDALsource(), grdarray[Band(1)]; force = true)) < 3e4
+            gdalarray = Raster(gdalfilename)
             # @test convert(ProjString, crs(gdalarray)) == convert(ProjString, EPSG(4326))
             @test val(dims(gdalarray, X)) ≈ val(dims(grdarray, X))
             @test val(dims(gdalarray, Y)) ≈ val(dims(grdarray, Y))
diff --git a/test/sources/ncdatasets.jl b/test/sources/ncdatasets.jl
index b505ab396..74c3144a9 100644
--- a/test/sources/ncdatasets.jl
+++ b/test/sources/ncdatasets.jl
@@ -364,7 +364,7 @@ end
             # Tiff locus = Start, Netcdf locus = Center
             @test index(gdalarray, Y) .+ 0.5 ≈ index(nccleaned, Y)
             @test index(gdalarray, X) .+ 1.0  ≈ index(nccleaned, X)
-            @test gdalarray ≈ nccleaned
+            @test parent(gdalarray) ≈ parent(nccleaned)
         end
 
         @testset "to grd" begin
@@ -376,7 +376,7 @@ end
             @test bounds(grdarray) == bounds(nccleaned)
             @test index(grdarray, Y) ≈ reverse(index(nccleaned, Y)) .- 0.5
             @test index(grdarray, X) ≈ index(nccleaned, X) .- 1.0
-            @test reverse(grdarray; dims=Y) ≈ nccleaned
+            @test parent(reverse(grdarray; dims=Y)) ≈ parent(nccleaned)
             rm("testgrd.gri")
             rm("testgrd.grd")
         end
@@ -494,7 +494,7 @@ end
         cp(ncmulti, ncmulti_custom, force=true)
         @time ncstack_custom = RasterStack(ncmulti_custom, source=Rasters.NCDsource())
         @test ncstack_custom isa RasterStack
-        @test map(read(ncstack_custom), read(ncstack)) do a, b
+        @test maplayers(read(ncstack_custom), read(ncstack)) do a, b
             all(a .=== b)
         end |> all
     end
diff --git a/test/sources/rasterdatasources.jl b/test/sources/rasterdatasources.jl
index 6cd00cf99..5c3049f32 100644
--- a/test/sources/rasterdatasources.jl
+++ b/test/sources/rasterdatasources.jl
@@ -1,4 +1,4 @@
-using Rasters, RasterDataSources, Test, Dates, ArchGDAL, NCDatasets
+using Rasters, RasterDataSources, Test, Dates, ArchGDAL, NCDatasets, Proj
 
 # Too big to test on CI
 # if !haskey(ENV, "CI")
@@ -33,6 +33,20 @@ end
     @test all(st.bio1 .=== A)
     @test st isa RasterStack
     @test A isa Raster
+    # Future Bioclim works
+    st = RasterStack(WorldClim{Future{BioClim, CMIP6, GFDL_ESM4, SSP370}}, (1, 2); 
+        date = Date(2050), res = "10m",
+        lazy=true, 
+        missingval=Inf, 
+        crs=nothing, 
+        mappedcrs=EPSG(4326),
+    )
+    @test missingval(st) === Inf32
+    @test missingval(st.bio1) === Inf32
+    ra = Raster(WorldClim{Future{BioClim, CMIP6, GFDL_ESM4, SSP370}}, 2; 
+        date = Date(2050), res = "10m"
+    )
+    @test Rasters.name(ra) == :bio2    
 end
 
 @testset "load CHELSA BioClim" begin
@@ -53,7 +67,7 @@ end
     # Allow forcing keywords
     st = RasterStack(CHELSA{BioClim}, (1, 2); 
          lazy=true, 
-         missingval=-Int16(9999), 
+         missingval= Int16(9999), 
          metadata=Rasters.NoMetadata(), 
          crs=nothing, 
          mappedcrs=EPSG(4326),
diff --git a/test/sources/zarr.jl b/test/sources/zarr.jl
index 2cac0ead1..9a3da8cb2 100644
--- a/test/sources/zarr.jl
+++ b/test/sources/zarr.jl
@@ -1,6 +1,6 @@
 using Rasters
-using Zarr
 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"
diff --git a/test/stack.jl b/test/stack.jl
index c8acc6e06..56f8dbc68 100644
--- a/test/stack.jl
+++ b/test/stack.jl
@@ -58,7 +58,7 @@ end
     @test collect(values(st)) == [raster1, raster2]
 end
 
-@testset "st fields " begin
+@testset "stack fields" begin
     @test DimensionalData.layerdims(st, :r1) == DimensionalData.format(dims1, data1)
     @test metadata(st) == NoMetadata()
     @test metadata(st, :r1) == NoMetadata()
@@ -87,7 +87,7 @@ end
         @test dims(s[:r2]) == (X(Sampled(10.0:10.0:100.0, ForwardOrdered(), Regular(10.0), Points(), NoMetadata())),
                                Y(Sampled(-10.0:10.0:10.0, ForwardOrdered(), Regular(10.0), Points(), NoMetadata())))
         @test refdims(s[:r2]) == 
-            (Ti(Sampled([DateTime(2019)], ForwardOrdered(), Irregular((DateTime(2019), DateTime(2019))), Points(), NoMetadata())),)
+            (Ti(Sampled([DateTime(2019)], ForwardOrdered(), Irregular(nothing, nothing), Points(), NoMetadata())),)
         @test isnothing(missingval(s, :r2)) && isnothing(missingval(s[:r2]))
     end
 
@@ -101,7 +101,7 @@ end
         @test dims(sv.r2) == (X(Sampled(10.0:10:100.0, ForwardOrdered(), Regular(10.0), Points(), NoMetadata())),
                                Y(Sampled(0.0:10:20.0, ForwardOrdered(), Regular(10.0), Points(), NoMetadata())))
         @test refdims(sv[:r2])[1] == 
-            Ti(Sampled(view([DateTime(2019)], 1:1), ForwardOrdered(), Irregular((DateTime(2019), DateTime(2019))), Points(), NoMetadata()))
+            Ti(Sampled(view([DateTime(2019)], 1:1), ForwardOrdered(), Irregular(nothing, nothing), Points(), NoMetadata()))
         # Stack of view-based Rasters
         v = view(st, X(2:4), Y(5:6))
         # @inferred view(st, X(2:4), Y(5:6))
diff --git a/test/test_utils.jl b/test/test_utils.jl
index 1557dfc76..688f993d2 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