Cdm overhaul

Project.toml

@@ -1,7 +1,7 @@
 name = "Rasters"
 uuid = "a3a2b9e3-a471-40c9-b274-f788e487c689"
 authors = ["Rafael Schouten <[email protected]>"]
-version = "0.14.3"
+version = "0.14.6"
 
 [deps]
 Adapt = "79e6a3ab-5dfb-504d-930d-738a2a938a0e"
@@ -22,6 +22,7 @@ LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 Missings = "e1d29d7a-bbdc-5cf2-9ac0-f12de2c33e28"
 Mmap = "a63ad114-7e13-5084-954f-fe012c677804"
 OffsetArrays = "6fe1bfb0-de20-5000-8ca7-80f57d26f881"
+OrderedCollections = "bac558e1-5e72-5ebc-8fee-abe8a469f55d"
 ProgressMeter = "92933f4c-e287-5a05-a399-4b506db050ca"
 RecipesBase = "3cdcf5f2-1ef4-517c-9805-6587b60abb01"
 Reexport = "189a3867-3050-52da-a836-e630ba90ab69"
@@ -61,7 +62,7 @@ CommonDataModel = "0.2.3, 0.3"
 ConstructionBase = "1"
 CoordinateTransformations = "0.6.2"
 DataFrames = "1"
-DimensionalData = "0.29.4"
+DimensionalData = "0.30.0"
 DiskArrays = "0.4"
 Extents = "0.1"
 FillArrays = "0.12, 0.13, 1"
@@ -71,12 +72,13 @@ GeoDataFrames = "0.3"
 GeoFormatTypes = "0.4"
 GeoInterface = "1.0"
 GeometryBasics = "0.4"
-GeometryOps = "0.1.19"
-GeometryOpsCore = "0.1.1"
-Makie = "0.20, 0.21, 0.22"
+GeometryOps = "0.1"
+GeometryOpsCore = "0.1"
+Makie = "0.20, 0.21, 0.22, 0.23, 0.24"
 Missings = "0.4, 1"
 NCDatasets = "0.13, 0.14"
 OffsetArrays = "1"
+OrderedCollections = "1.8"
 Plots = "1"
 ProgressMeter = "1"
 Proj = "1.7.2"
@@ -119,5 +121,4 @@ test = ["Aqua", "ArchGDAL", "CFTime", "CoordinateTransformations", "DataFrames",
 
 
 [sources]
-DimensionalData = {url = "https://github.com/rafaqz/DimensionalData.jl", rev = "as/individualindexing"}
-GeometryOps = {url = "https://github.com/JuliaGeo/GeometryOps.jl", rev = "as/extentforwarding_for_predicates"}
+DimensionalData = {url = "https://github.com/rafaqz/DimensionalData.jl", rev = "as/multidimtrait"}

README.md

@@ -1,22 +1,26 @@
 # Rasters
-
+<img src="docs/src/assets/logo.png" align="right" width="30%"></img>
 [![License: MIT](https://img.shields.io/badge/License-MIT-green.svg)](https://github.com/rafaqz/Rasters.jl/blob/main/LICENSE)
 [![](https://img.shields.io/badge/docs-stable-blue.svg)](https://rafaqz.github.io/Rasters.jl/stable)
 [![](https://img.shields.io/badge/docs-dev-blue.svg)](https://rafaqz.github.io/Rasters.jl/dev)
 [![CI](https://github.com/rafaqz/Rasters.jl/actions/workflows/ci.yml/badge.svg)](https://github.com/rafaqz/Rasters.jl/actions/workflows/ci.yml)
 [![Codecov](https://codecov.io/gh/rafaqz/Rasters.jl/branch/main/graph/badge.svg)](https://codecov.io/gh/rafaqz/Rasters.jl)
-[![Aqua.jl Quality Assurance](https://img.shields.io/badge/Aquajl-%F0%9F%8C%A2-aqua.svg)](https://github.com/JuliaTesting/Aqua.jl)
-[![Downloads](https://shields.io/endpoint?url=https://pkgs.genieframework.com/api/v1/badge/Rasters&label=Downloads)](https://pkgs.genieframework.com?packages=Rasters)
-
-<img src="/docs/src/assets/logo.png" align="right" width="30%"></img>
+[![Aqua QA](https://raw.githubusercontent.com/JuliaTesting/Aqua.jl/master/badge.svg)](https://github.com/JuliaTesting/Aqua.jl)
+[![Downloads](https://img.shields.io/badge/dynamic/json?url=http%3A%2F%2Fjuliapkgstats.com%2Fapi%2Fv1%2Fmonthly_downloads%2FRasters&query=total_requests&suffix=%2Fmonth&label=Downloads)](https://juliapkgstats.com/pkg/Rasters)
 
-[Rasters.jl](https://rafaqz.github.io/Rasters.jl/dev) defines common types and methods for reading, writing and
-manipulating rasterized spatial data. 
+[Rasters.jl](https://rafaqz.github.io/Rasters.jl/dev) is a powerful Julia package for working with spatial raster data. It provides a unified interface for reading, writing, and manipulating raster data. The package extends [DimensionalData.jl](https://rafaqz.github.io/DimensionalData.jl/dev/) to enable intuitive spatial indexing and manipulation of raster data.
 
-These currently include raster arrays like GeoTIFF and NetCDF, R grd files, 
-multi-layered stacks, and multi-file series of arrays and stacks. 
+Key features:
+- Support for multiple raster formats (e.g. GeoTIFF, NetCDF, GRD)
+- Support for multi-layered stacks and multi-file series of arrays
+- Lazy loading of large datasets
+- Intuitive spatial indexing with named dimensions (X, Y, Time)
+- Efficient handling of multi-layered stacks and time series
+- Built-in support for coordinate reference systems (CRS)
+- High-performance operations optimized for spatial data
 
 # Quick start
+
 Install the package by typing:
 
 ```julia
@@ -31,7 +35,7 @@ using Rasters
 
 Using `Rasters` to read GeoTiff or NetCDF files will output something similar to the
 following toy examples. This is possible because Rasters.jl extends
-[DimensionalData.jl](https://github.com/rafaqz/DimensionalData.jl) so that
+[DimensionalData.jl](https://rafaqz.github.io/DimensionalData.jl/dev/) so that
 spatial data can be indexed using named dimensions like `X`, `Y` and `Ti` (time)
 and e.g. spatial coordinates.
 
@@ -41,7 +45,7 @@ lon, lat = X(25:1:30), Y(25:1:30)
 ti = Ti(DateTime(2001):Month(1):DateTime(2002))
 ras = Raster(rand(lon, lat, ti)) # this generates random numbers with the dimensions given
 ```
-```
+```julia
 6×6×13 Raster{Float64,3} with dimensions: 
   X Sampled{Int64} 25:1:30 ForwardOrdered Regular Points,
   Y Sampled{Int64} 25:1:30 ForwardOrdered Regular Points,
@@ -64,7 +68,7 @@ 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:
+For example, to use the GDAL backend, and download RasterDataSources files, you need to do:
 
 ```julia
 using Rasters, ArchGDAL, RasterDataSources
@@ -80,7 +84,7 @@ Sources and packages needed:
 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`
+- Coordinate transformations for GDAL rasters: `using CoordinateTransformations`
 
 ## Getting the `lookup` array from dimensions
 
@@ -100,7 +104,7 @@ Selecting a time slice by `index` is done via
 ```julia
 ras[Ti(1)]
 ```
-```
+```julia
 6×6 Raster{Float64,2} with dimensions: 
   X Sampled{Int64} 25:1:30 ForwardOrdered Regular Points,
   Y Sampled{Int64} 25:1:30 ForwardOrdered Regular Points
@@ -120,7 +124,7 @@ values:      25         26          27          28         29          30
 ```julia
 ras[Ti=1]
 ```
-```
+```julia
 6×6 Raster{Float64,2} with dimensions: 
   X Sampled{Int64} 25:1:30 ForwardOrdered Regular Points,
   Y Sampled{Int64} 25:1:30 ForwardOrdered Regular Points
@@ -142,7 +146,7 @@ or and interval of indices using the syntax `=a:b` or `(a:b)`
 ```julia
 ras[Ti(1:10)]
 ```
-```
+```julia
 6×6×10 Raster{Float64,3} with dimensions: 
   X Sampled{Int64} 25:1:30 ForwardOrdered Regular Points,
   Y Sampled{Int64} 25:1:30 ForwardOrdered Regular Points,
@@ -165,7 +169,7 @@ values: [:, :, 1]
 ```julia
 ras[Ti=At(DateTime(2001))]
 ```
-```
+```julia
 6×6 Raster{Float64,2} with dimensions: 
   X Sampled{Int64} 25:1:30 ForwardOrdered Regular Points,
   Y Sampled{Int64} 25:1:30 ForwardOrdered Regular Points

docs/make.jl

@@ -3,6 +3,7 @@ using Documenter, Rasters, Plots, Logging, Statistics, Dates,
 import Makie, CairoMakie
 using DocumenterVitepress
 using Rasters.LookupArrays, Rasters.Dimensions
+import Shapefile, DataFrames, NaturalEarth # to avoid precompilation in doctests
 
 # Don't output huge svgs for Makie plots
 CairoMakie.activate!(type = "png")
@@ -43,17 +44,15 @@ makedocs(
         devbranch = "main",
         devurl = "dev";
     ),
-    draft = false,
     source = "src",
     build = "build",
-    warnonly=true,
+    warnonly=false,
 )
 
 # Enable logging to console again
 Logging.disable_logging(Logging.BelowMinLevel)
 
-deploydocs(; repo="github.com/rafaqz/Rasters.jl",
-    target = "build", # this is where Vitepress stores its output
+DocumenterVitepress.deploydocs(; repo="github.com/rafaqz/Rasters.jl",
     branch = "gh-pages",
     devbranch = "main",
     push_preview = true

docs/src/.vitepress/config.mts

@@ -99,6 +99,9 @@ export default defineConfig({
   ],
   ignoreDeadLinks: false,
   vite: {
+    define: {
+      __DEPLOY_ABSPATH__: JSON.stringify('REPLACE_ME_DOCUMENTER_VITEPRESS_DEPLOY_ABSPATH'),
+    },
     resolve: {
       alias: {
         '@': path.resolve(__dirname, '../components')

docs/src/tutorials/spatial_mean.md

@@ -147,7 +147,7 @@ As a next step, we would like to know how precipitation will change in Chile unt
 
 To start, we define a simple function that takes an SSP (socioeconomic scenario) and a GCM (climate model) as input, and return the appropriate climate data. 
 
-````@example zonal
+````@example cellarea
 using Dates
 getfutureprec(ssp, gcm) = Raster(WorldClim{Future{Climate, CMIP6, gcm, ssp}}, :prec, date = Date(2090))
 ````
@@ -163,11 +163,10 @@ GCMs = Dim{:gcm}([GFDL_ESM4, IPSL_CM6A_LR]) # These are different general circul
 precip_future = (@d getfutureprec.(SSPs, GCMs)) |> RasterSeries |> Rasters.combine
 ````
 
-Since the format of WorldClim's datasets for future climate is slightly different from the dataset for the historical period, this actually returned a 5-dimensional raster, with a `Band` dimension that represents months. Here we'll just select the 6th month, matching the selection above (but note that the analysis would also work for all Bands simultaneously). We will also replace the `NaN` missing value by the more standard `missing` using [`replace_missing`](@ref). 
+Since the format of WorldClim's datasets for future climate is slightly different from the dataset for the historical period, this actually returned a 5-dimensional raster, with a `Band` dimension that represents months. Here we'll just select the 6th month, matching the selection above (but note that the analysis would also work for all Bands simultaneously).
 
 ````@example cellarea
 precip_future = precip_future[Band = 6]
-precip_future = replace_missing(precip_future)
 ````
 
 On our 4-dimensional raster, functions like `crop` and `mask`, as well as broadcasting, will still work.

ext/RastersGRIBDatasetsExt/gribdatasets_source.jl

@@ -21,8 +21,7 @@ Base.close(os::RA.OpenStack{GRIBsource}) = nothing
 RA.missingval(var::GDS.Variable, ::RA.Metadata{<:RA.CDMsource}) = _missingval(var)
 RA.missingval(var::GDS.Variable, args...) = _missingval(var)
 
-function _missingval(var::GDS.Variable{T}) where T
+function _missingval(var::GDS.Variable)
     mv = GDS.missing_value(var)
-    T1 = promote_type(typeof(mv), T)
-    return T1(mv)
+    RA._fix_missingval(var, mv)
 end

src/Rasters.jl

@@ -1,12 +1,5 @@
 module Rasters
 
-# Use the README as the module docs
-@doc let
-    # path = joinpath(dirname(@__DIR__), "README.md")
-    # include_dependency(path)
-    # read(path, String)
-end Rasters
-
 using Dates
 
 # Load first to fix StaticArrays invalidations
@@ -24,6 +17,7 @@ import Adapt,
        GeometryOps,
        GeometryOpsCore,
        OffsetArrays,
+       OrderedCollections,
        ProgressMeter,
        Missings,
        Mmap,
@@ -44,6 +38,7 @@ using DimensionalData: Name, NoName
 using .Dimensions: StandardIndices, DimTuple
 using .Lookups: LookupTuple
 
+using OrderedCollections: OrderedDict
 using Statistics: mean
 using RecipesBase: @recipe, @series
 using Base: tail, @propagate_inbounds
@@ -123,7 +118,6 @@ include("utils.jl")
 include("skipmissing.jl")
 
 include("geometry_lookup/geometry_lookup.jl")
-include("geometry_lookup/lookups.jl")
 include("geometry_lookup/methods.jl")
 include("geometry_lookup/io.jl")
 

src/array.jl

@@ -329,7 +329,7 @@ function Raster(filename::AbstractString; source=nokw, kw...)
     end::Raster
 end
 # Load a Raster from an opened Dataset
-# We need the inner method for AbstractArray ambiguit
+# We need the inner method for AbstractArray ambiguity
 Raster(ds; kw...) = _raster(ds; kw...)
 function _raster(ds;
     dims=nokw,
@@ -401,4 +401,4 @@ filekey(filename::String) = Symbol(splitext(basename(filename))[1])
 
 # Add a `dimconstructor` method so `AbstractProjected` lookups create a Raster
 # TODO this should be unwrapped to `DD.lookupconstructor` to avoid future ambiguities
-DD.dimconstructor(::Tuple{<:Dimension{<:AbstractProjected},Vararg{Dimension}}) = Raster
+DD.dimconstructor(::Tuple{<:Dimension{<:AbstractProjected},Vararg{Dimension}}) = Raster

src/create.jl

@@ -248,8 +248,10 @@ function create(filename::Nothing, types::NamedTuple, dims::Tuple;
     layerdims=nokw,
     layermetadata=nokw,
     f=identity,
+    lazy=false,
     kw...
 )
+    lazy && throw(ArgumentError("`lazy` cannot be `true` without passing a `filename` keyword"))
     layerdims = isnokw(layerdims) ? map(_ -> basedims(dims), types) : layerdims
     layermetadata = layermetadata isa NamedTuple ? layermetadata : map(_ -> layermetadata, types)
     layerfill = fill isa NamedTuple ? fill : map(_ -> fill, types)
@@ -355,4 +357,4 @@ end
 
 _warn_keyword_not_used(label, obj) = @warn "`$label` of `$obj` found. But `chunks` are not used for in-memory rasters"
 _missingval_pair(missingval::Pair) = missingval
-_missingval_pair(missingval) = missingval => missingval
+_missingval_pair(missingval) = missingval => missingval

src/dimensions.jl

@@ -19,6 +19,28 @@ mean(A; dims=Band)
 """
 @dim Band
 
+"""
+    Lat <: Dimension
+
+    Lat(val=:)
+
+Used for holding degrees north lookups.
+
+Will error on lookup construction if metadata of `units="degrees_north"` does not exist.
+"""
+@dim Lat
+
+"""
+    Lon <: Dimension
+
+    Lon(val=:)
+
+Used for holding degrees east lookups.
+
+Will error on lookup construction if metadata of `units="degrees_east"` does not exist.
+"""
+@dim Lon
+
 """
     Geometry <: Dimension
 
@@ -37,4 +59,4 @@ val = A[Geometry(Touches(other_geom))] # this is automatically accelerated by sp
 mean(A; dims=Geometry)
 ```
 """
-@dim Geometry
+@dim Geometry

src/filearray.jl

@@ -49,6 +49,18 @@ function FileArray{S}(
     T = _mod_eltype(var, mod)
     return FileArray{S,T,N}(filename, size(var); eachchunk, haschunks, mod, kw...)
 end
+function FileArray{S}(
+    var::AbstractArray{Char,N}, filename; mod, kw...
+) where {S<:CDMsource,N}
+    eachchunk = DA.eachchunk(var)
+    haschunks = DA.haschunks(var)
+    if Missings.nonmissingtype(eltype(mod)) isa AbstractString
+        return FileArray{S,eltype(mod),N-1}(filename, size(var); eachchunk, haschunks, mod, kw...)
+    else
+        T = _mod_eltype(var, mod)
+        return FileArray{S,T,N}(filename, size(var); eachchunk, haschunks, mod, kw...)
+    end
+end
 
 # FileArray has S, T and N parameters not recoverable from fields
 ConstructionBase.constructorof(::Type{<:FileArray{S,T,N}}) where {S,T,N} = FileArray{S,T,N}
@@ -68,7 +80,12 @@ end
 function DA.readblock!(A::FileArray, dst, r::AbstractUnitRange...)
     open(A) do O
         if isdisk(O)
-            DA.readblock!(O, dst, r...)
+            # Handle CF 2d Char arrays that are really 1d strings
+            if eltype(O) <: Char && eltype(A) <: String
+                DA.readblock!(DiskCharToString(O), dst, r...)
+            else
+                DA.readblock!(O, dst, r...)
+            end
         else
             dest[r...] .= view(parent(O), r...)
         end
@@ -77,7 +94,11 @@ end
 function DA.writeblock!(A::FileArray, src, r::AbstractUnitRange...)
     open(A; write=A.write) do O
         if isdisk(A)
-            DA.writeblock!(O, src, r...)
+            if eltype(O) <: Char && eltype(A) <: String
+                DA.writeblock!(DiskCharToString(O), src, r...)
+            else
+                DA.writeblock!(O, src, r...)
+            end
         else
             parent(O)[r...] .= src
         end