Merge branch 'main' into test_unitful_cellarea

.github/ISSUE_TEMPLATE/bug-report.md → .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

.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:

Project.toml

@@ -34,6 +34,7 @@ Makie = "ee78f7c6-11fb-53f2-987a-cfe4a2b5a57a"
 NCDatasets = "85f8d34a-cbdd-5861-8df4-14fed0d494ab"
 Proj = "c94c279d-25a6-4763-9509-64d165bea63e"
 RasterDataSources = "3cb90ccd-e1b6-4867-9617-4276c8b2ca36"
+StatsBase = "2913bbd2-ae8a-5f71-8c99-4fb6c76f3a91"
 ZarrDatasets = "519a4cdf-1362-424a-9ea1-b1d782dbb24b"
 
 [extensions]
@@ -44,6 +45,7 @@ RastersMakieExt = "Makie"
 RastersNCDatasetsExt = "NCDatasets"
 RastersProjExt = "Proj"
 RastersRasterDataSourcesExt = "RasterDataSources"
+RastersStatsBaseExt = "StatsBase"
 RastersZarrDatasetsExt = "ZarrDatasets"
 
 [compat]
@@ -81,6 +83,7 @@ SafeTestsets = "0.1"
 Setfield = "0.6, 0.7, 0.8, 1"
 Shapefile = "0.10, 0.11"
 Statistics = "1"
+StatsBase = "0.34"
 Test = "1"
 ZarrDatasets = "0.1"
 julia = "1.10"
@@ -101,9 +104,12 @@ Proj = "c94c279d-25a6-4763-9509-64d165bea63e"
 RasterDataSources = "3cb90ccd-e1b6-4867-9617-4276c8b2ca36"
 SafeTestsets = "1bc83da4-3b8d-516f-aca4-4fe02f6d838f"
 Shapefile = "8e980c4a-a4fe-5da2-b3a7-4b4b0353a2f4"
+StableRNGs = "860ef19b-820b-49d6-a774-d7a799459cd3"
+StatsBase = "2913bbd2-ae8a-5f71-8c99-4fb6c76f3a91"
 Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"
 Unitful = "1986cc42-f94f-5a68-af5c-568840ba703d"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 
 [targets]
-test = ["Aqua", "ArchGDAL", "CFTime", "CoordinateTransformations", "DataFrames", "GeoDataFrames", "GeometryBasics", "GRIBDatasets", "NCDatasets", "Plots", "Proj", "RasterDataSources", "SafeTestsets", "Shapefile", "Statistics", "Test", "Unitful", "ZarrDatasets"]
+test = ["Aqua", "ArchGDAL", "CFTime", "CoordinateTransformations", "DataFrames", "GeoDataFrames", "GeometryBasics", "GRIBDatasets", "NCDatasets", "Plots", "Proj", "RasterDataSources", "SafeTestsets", "Shapefile", "StableRNGs", "Statistics", "StatsBase", "Test", "Unitful", "ZarrDatasets"]
+

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

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

src/array.jl

@@ -269,6 +269,7 @@ function Raster(A::AbstractArray{T,1}, dims::Tuple{<:Dimension,<:Dimension,Varar
 )::Raster{T,length(dims)} where T
     Raster(reshape(A, map(length, dims)), dims; kw...)
 end
+Raster(A::AbstractArray{<:Any,1}, dim::Dimension; kw...) = Raster(A, (dim,); kw...)
 function Raster(table, dims::Tuple;
     name=nokw,
     kw...

src/extensions.jl

@@ -228,6 +228,50 @@ Note that cellarea returns the area in square m, while cellsize still uses squar
     return cellarea(args...; kw..., radius = 6371.0088)
 end
 
+"""
+Rasters.sample([rng], x, [n::Integer]; kw...)
+
+Sample `n` random and optionally weighted points from from a `Raster` or `RasterStack`.
+Returns a `Vector` of `NamedTuple`, closely resembling the return type of [`extract`](@ref).
+
+Run `using StatsBase` to make this method available.
+Note that this function is not exported to avoid confusion with StatsBase.sample
+
+# Keywords
+
+- `geometry`: include `:geometry` in returned `NamedTuple`. Specify the type and dimensions of the returned geometry by
+    providing a `Tuple` or `NamedTuple` of dimensions. Defaults to `(X,Y)`
+- `index`: include `:index` of the `CartesianIndex` in returned `NamedTuple`, `false` by default.
+- `name`: a `Symbol` or `Tuple` of `Symbol` corresponding to layer/s of a `RasterStack` to extract. All layers by default.
+- `skipmissing`: skip missing points automatically.
+- `weights`: A DimArray that matches one or more of the dimensions of `x` with weights for sampling.
+- `weightstype`: a `StatsBase.AbstractWeights` specifying the type of weights. Defaults to `StatsBase.Weights`.
+- `replace`: sample with replacement, `true` by default. See `StatsBase.sample`
+- `ordered`: sample in order, `false` by default. See `StatsBase.sample`
+
+# Example
+This code draws 5 random points from a raster, weighted by cell area.
+```julia
+using Rasters, Rasters.Lookups, Proj, StatsBase
+xdim = X(Projected(90.0:10.0:120; sampling=Intervals(Start()), crs=EPSG(4326)))
+ydim = Y(Projected(0.0:10.0:50; sampling=Intervals(Start()), crs=EPSG(4326)))
+myraster = rand(xdim, ydim)
+Rasters.sample(myraster, 5; weights = cellarea(myraster))
+
+# output
+
+5-element Vector{@NamedTuple{geometry::Tuple{Float64, Float64}, ::Union{Missing, Float64}}}:
+ @NamedTuple{geometry::Tuple{Float64, Float64}, ::Union{Missing, Float64}}(((90.0, 10.0), 0.7360504790189618))
+ @NamedTuple{geometry::Tuple{Float64, Float64}, ::Union{Missing, Float64}}(((90.0, 30.0), 0.5447657183842469))
+ @NamedTuple{geometry::Tuple{Float64, Float64}, ::Union{Missing, Float64}}(((90.0, 30.0), 0.5447657183842469))
+ @NamedTuple{geometry::Tuple{Float64, Float64}, ::Union{Missing, Float64}}(((90.0, 10.0), 0.7360504790189618))
+ @NamedTuple{geometry::Tuple{Float64, Float64}, ::Union{Missing, Float64}}(((110.0, 10.0), 0.5291143028176258))
+```
+"""
+sample(args...; kw...) = throw_extension_error(sample, "StatsBase", :RastersStatsBaseExt, args)
+
+
+
 # Other shared stubs
 function layerkeys end
 function smapseries end

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) = 

src/sources/commondatamodel.jl

@@ -241,7 +241,14 @@ function _layermetadata(ds::AbstractDataset; layers)
     map(layers.attrs) do attr
         md = _metadatadict(_sourcetrait(ds), attr)
         if haskey(attr, "grid_mapping")
-            md["grid_mapping"] = Dict(CDM.attribs(ds[attr["grid_mapping"]]))
+            if haskey(ds, attr["grid_mapping"])
+                md["grid_mapping"] = Dict(CDM.attribs(ds[attr["grid_mapping"]]))
+            else
+                global_attrs = CDM.attribs(ds)
+                if haskey(global_attrs, attr["grid_mapping"])
+                    md["grid_mapping"] = global_attrs["grid_mapping"]
+                end
+            end     
         end
         md
     end
@@ -571,4 +578,4 @@ function _def_dim_var!(ds::AbstractDataset, dim::Dimension)
     end
     CDM.defVar(ds, dimname, Vector(index(dim)), (dimname,); attrib=attrib)
     return nothing
-end
+end

src/sources/sources.jl

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