diff --git a/.github/dependabot.yml b/.github/dependabot.yml new file mode 100644 index 000000000..6867e71eb --- /dev/null +++ b/.github/dependabot.yml @@ -0,0 +1,11 @@ +# To get started with Dependabot version updates, you'll need to specify which +# package ecosystems to update and where the package manifests are located. +# Please see the documentation for all configuration options: +# https://docs.github.com/code-security/dependabot/dependabot-version-updates/configuration-options-for-the-dependabot.yml-file + +version: 2 +updates: + - package-ecosystem: "github-actions" + directory: "/" + schedule: + interval: "weekly" diff --git a/Project.toml b/Project.toml index 8ef210f98..5b580da91 100644 --- a/Project.toml +++ b/Project.toml @@ -60,7 +60,7 @@ ConstructionBase = "1" CoordinateTransformations = "0.6.2" DataFrames = "1" DimensionalData = "0.29.4" -DiskArrays = "0.3, 0.4" +DiskArrays = "0.4" Extents = "0.1" FillArrays = "0.12, 0.13, 1" Flatten = "0.4" diff --git a/ext/RastersStatsBaseExt/sample.jl b/ext/RastersStatsBaseExt/sample.jl index f54f9d919..80180fa6f 100644 --- a/ext/RastersStatsBaseExt/sample.jl +++ b/ext/RastersStatsBaseExt/sample.jl @@ -52,6 +52,7 @@ function _getindex(::Type{T}, x::AbstractRasterStack{<:Any, NT}, dims, idx) wher RA._maybe_add_fields( T, NT(x[RA.commondims(idx, x)]), + nothing, DimPoints(dims)[RA.commondims(idx, dims)], val(idx) ) diff --git a/src/create.jl b/src/create.jl index 25cdba8ae..a321a9ad2 100644 --- a/src/create.jl +++ b/src/create.jl @@ -291,15 +291,15 @@ function create(filename::AbstractString, source::Source, ::Type{T}, dims::Tuple end # Create layers of zero arrays rast = Raster(A, dims; name, missingval=mv_inner) + # Write to disk Rasters.write(f, filename, source, rast; - eltype, chunks, metadata, scale, offset, missingval, verbose, force, coerce, write, kw... + eltype, chunks, metadata, scale, offset, missingval=mv_inner, verbose, force, coerce, write, kw... ) do W - # write returns a variable, wrap it as a Raster + # `write`` returns a variable, wrap it as a Raster and run f f(rebuild(rast; data=W)) end # Don't pass in `missingval`, read it again from disk in case it changed - r = Raster(filename; source, lazy, metadata, dropband, coerce, missingval=mv_outer) - return r + return Raster(filename; source, lazy, metadata, dropband, coerce, missingval=mv_outer) end # Create on-disk RasterStack from filename, source, type and dims function create(filename::AbstractString, source::Source, layertypes::NamedTuple, dims::Tuple; diff --git a/src/lookup.jl b/src/lookup.jl index 07d1a68e0..75fd24660 100644 --- a/src/lookup.jl +++ b/src/lookup.jl @@ -95,15 +95,15 @@ dim(lookup::Projected) = lookup.dim end @inline function LA.selectindices(l::Projected, sel::LA.Selector{<:AbstractVector}; kw...) selval = reproject(mappedcrs(l), crs(l), dim(l), val(sel)) - _selectvec(l, rebuild(sel; val=selval)sel; kw...) + LA._selectvec(l, rebuild(sel; val=selval); kw...) end @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...) + LA._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...) + LA._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)) diff --git a/src/methods/burning/allocs.jl b/src/methods/burning/allocs.jl index 5d916cf76..7956e266f 100644 --- a/src/methods/burning/allocs.jl +++ b/src/methods/burning/allocs.jl @@ -11,22 +11,36 @@ function Allocs(buffer) return Allocs(buffer, edges, scratch, crossings) end -function _burning_allocs(x; nthreads=_nthreads(), threaded=true, kw...) - dims = commondims(x, DEFAULT_POINT_ORDER) +function _burning_allocs(x; + nthreads=_nthreads(), + threaded=true, + burncheck_metadata=Metadata(), + kw... +) if threaded - [Allocs(_init_bools(dims; metadata=Metadata())) for _ in 1:nthreads] + if isnothing(x) + [Allocs(nothing) for _ in 1:nthreads] + else + dims = commondims(x, DEFAULT_POINT_ORDER) + [Allocs(_init_bools(dims; metadata=deepcopy(burncheck_metadata))) for _ in 1:nthreads] + end else - Allocs(_init_bools(dims; metadata=Metadata())) + if isnothing(x) + Allocs() + else + dims = commondims(x, DEFAULT_POINT_ORDER) + Allocs(_init_bools(dims; metadata=burncheck_metadata)) + end end end _get_alloc(allocs::Vector{<:Allocs}) = _get_alloc(allocs[Threads.threadid()]) _get_alloc(allocs::Allocs) = allocs - # TODO include these in Allocs _alloc_burnchecks(n::Int) = fill(false, n) _alloc_burnchecks(x::AbstractArray) = _alloc_burnchecks(length(x)) + function _set_burnchecks(burnchecks, metadata::Metadata{<:Any,<:Dict}, verbose) metadata["missed_geometries"] = .!burnchecks verbose && _burncheck_info(burnchecks) diff --git a/src/methods/burning/array_init.jl b/src/methods/burning/array_init.jl index a700a6454..306a26e61 100644 --- a/src/methods/burning/array_init.jl +++ b/src/methods/burning/array_init.jl @@ -4,29 +4,22 @@ # TODO merge this with `create` somehow _init_bools(to; kw...) = _init_bools(to, BitArray; kw...) _init_bools(to, T::Type; kw...) = _init_bools(to, T, nothing; kw...) -_init_bools(to::AbstractRasterSeries, T::Type, data; kw...) = _init_bools(first(to), T, data; kw...) -_init_bools(to::AbstractRasterStack, T::Type, data; kw...) = _init_bools(first(to), T, data; kw...) -_init_bools(to::AbstractRaster, T::Type, data; kw...) = _init_bools(to, dims(to), T, data; kw...) -_init_bools(to::DimTuple, T::Type, data; kw...) = _init_bools(to, to, T, data; kw...) -function _init_bools(to::Nothing, T::Type, data; - geometrycolumn=nothing, - collapse=nokw, - res=nokw, - size=nokw, - kw... -) +_init_bools(to::AbstractRasterSeries, T::Type, data; kw...) = + _init_bools(dims(first(to)), T, data; kw...) +_init_bools(to::AbstractRasterStack, T::Type, data; kw...) = + _init_bools(dims(to), dims(to), T, data; kw...) +_init_bools(to::AbstractRaster, T::Type, data; kw...) = + _init_bools(to, dims(to), T, data; kw...) +_init_bools(to::DimTuple, T::Type, data; kw...) = + _init_bools(to, to, T, data; kw...) +function _init_bools(to::Nothing, T::Type, data; geometrycolumn=nothing, kw...) # Get the extent of the geometries ext = _extent(data; geometrycolumn) isnothing(ext) && throw(ArgumentError("no recognised dimensions, extent or geometry")) - return _init_bools(ext, T, data; collapse, res, size, kw...) -end -function _init_bools(to::Extents.Extent, T::Type, data; - collapse=nokw, size=nokw, res=nokw, sampling=nokw, kw... -) - # Convert the extent to dims (there must be `res` or `size` in `kw`) - dims = _extent2dims(to; size, res, sampling, kw...) - _init_bools(to, dims, T, data; collapse, kw...) + return _init_bools(ext, T, data; kw...) end +_init_bools(to::Extents.Extent, T::Type, data; kw...) = + _init_bools(to, _extent2dims(to; kw...), T, data; kw...) function _init_bools(to, dims::DimTuple, T::Type, data; collapse::Union{Bool,Nothing,NoKW}=nokw, kw... ) @@ -54,12 +47,12 @@ function _alloc_bools(to, dims::DimTuple, ::Type{<:Array{T}}; missingval=false, metadata=NoMetadata(), kw... ) where T # Use an Array - data = fill!(Raster{T}(undef, dims), missingval) + data = fill!(Raster{T}(undef, dims), missingval) return rebuild(data; missingval, metadata) end -function _prepare_for_burning(B, locus=Center()) - B1 = _forward_ordered(B) +function _prepare_for_burning(B; locus=Center(), order=ForwardOrdered()) + B1 = _maybe_lazy_reorder(order, B) start_dims = map(dims(B1, DEFAULT_POINT_ORDER)) do d # Shift lookup values to center of pixels d = DD.maybeshiftlocus(locus, d) @@ -69,9 +62,12 @@ function _prepare_for_burning(B, locus=Center()) end # Convert to Array if its not one already -_lookup_as_array(d::Dimension) = parent(lookup(d)) isa Array ? d : modify(Array, d) +_lookup_as_array(x) = setdims(x, _lookup_as_array(dims(x))) +_lookup_as_array(dims::Tuple) = map(_lookup_as_array, dims) +_lookup_as_array(d::Dimension) = parent(lookup(d)) isa Array ? d : modify(Array, d) -function _forward_ordered(B) +_maybe_lazy_reorder(::Nothing, B) = B +function _maybe_lazy_reorder(::ForwardOrdered, B) reduce(dims(B); init=B) do A, d if DD.order(d) isa ReverseOrdered A = view(A, rebuild(d, lastindex(d):-1:firstindex(d))) diff --git a/src/methods/burning/geometry.jl b/src/methods/burning/geometry.jl index fa61c1f24..2ee27211b 100644 --- a/src/methods/burning/geometry.jl +++ b/src/methods/burning/geometry.jl @@ -22,10 +22,11 @@ function _burn_geometry!(B::AbstractRaster, trait::Nothing, data; lock=Threads.SpinLock(), verbose=true, progress=true, - threaded=true, + threaded=false, fill=true, - allocs=_burning_allocs(B; threaded), geometrycolumn=nothing, + burncheck_metadata=Metadata(), + allocs=_burning_allocs(B; threaded, burncheck_metadata), kw... )::Bool geoms = _get_geometries(data, geometrycolumn) @@ -36,8 +37,8 @@ function _burn_geometry!(B::AbstractRaster, trait::Nothing, data; geom = _getgeom(geoms, i) ismissing(geom) && return nothing a = _get_alloc(allocs) - B1 = a.buffer - burnchecks[i] = _burn_geometry!(B1, geom; fill, allocs=a, lock, kw...) + buffer = a.buffer + burnchecks[i] = _burn_geometry!(buffer, geom; fill, allocs=a, lock, kw...) return nothing end if fill @@ -99,14 +100,17 @@ function _burn_geometry!(B::AbstractRaster, ::GI.AbstractGeometryTrait, geom; return false end # Take a view of the geometry extent - B1 = view(B, Touches(geomextent)) - buf1 = _init_bools(B1; missingval=false) + V = view(B, Touches(geomextent)) + buffer = _init_bools(V; missingval=false) # Burn the polygon into the buffer - allocs = isnothing(allocs) ? Allocs(B) : allocs - hasburned = _burn_polygon!(buf1, geom; shape, geomextent, allocs, boundary, kw...) - @inbounds for i in eachindex(B1) - if buf1[i] - B1[i] = fill + # We allocate a new bitarray for the view for performance + # and always fill with `true`. + allocs = isnothing(allocs) ? Allocs(nothing) : allocs + hasburned = _burn_polygon!(buffer, geom; shape, geomextent, allocs, boundary, kw...) + # We then transfer burned `true` values to B via the V view + @inbounds for i in eachindex(V) + if buffer[i] + V[i] = fill end end else diff --git a/src/methods/burning/line.jl b/src/methods/burning/line.jl index 9155479bb..24f924f6c 100644 --- a/src/methods/burning/line.jl +++ b/src/methods/burning/line.jl @@ -1,33 +1,53 @@ # _burn_lines! # Fill a raster with `fill` where pixels touch lines in a geom -# Separated for a type stability function barrier -function _burn_lines!(B::AbstractRaster, geom; fill=true, kw...) - _check_intervals(B) +# Usually `fill` is `true` of `false` +function _burn_lines!( + B::AbstractRaster, geom; fill=true, verbose=false, kw... +) + _check_intervals(B, verbose) B1 = _prepare_for_burning(B) - return _burn_lines!(B1, geom, fill) + + xdim, ydim = dims(B, DEFAULT_POINT_ORDER) + regular = map((xdim, ydim)) do d + # @assert (parent(lookup(d)) isa Array) + lookup(d) isa AbstractSampled && span(d) isa Regular + end + msg = """ + Can only fill lines where dimensions have `Regular` lookups. + Consider using `boundary=:center`, reprojecting the crs, + or make an issue in Rasters.jl on github if you need this to work. + """ + all(regular) || throw(ArgumentError(msg)) + + # Set indices of B as `fill` when a cell is found to burn. + _burn_lines!(identity, dims(B1), geom) do D + @inbounds B1[D] = fill + end end -_burn_lines!(B, geom, fill) = - _burn_lines!(B, GI.geomtrait(geom), geom, fill) -function _burn_lines!(B::AbstractArray, ::Union{GI.MultiLineStringTrait}, geom, fill) +_burn_lines!(f::F, c::C, dims::Tuple, geom) where {F<:Function,C<:Function} = + _burn_lines!(f, c, dims, GI.geomtrait(geom), geom) +function _burn_lines!( + f::F, c::C, dims::Tuple, ::Union{GI.MultiLineStringTrait}, geom +) where {F<:Function,C<:Function} n_on_line = 0 for linestring in GI.getlinestring(geom) - n_on_line += _burn_lines!(B, linestring, fill) + n_on_line += _burn_lines!(f, c, dims, linestring) end return n_on_line end function _burn_lines!( - B::AbstractArray, ::Union{GI.MultiPolygonTrait,GI.PolygonTrait}, geom, fill -) + f::F, c::C, dims::Tuple, ::Union{GI.MultiPolygonTrait,GI.PolygonTrait}, geom +) where {F<:Function,C<:Function} n_on_line = 0 for ring in GI.getring(geom) - n_on_line += _burn_lines!(B, ring, fill) + n_on_line += _burn_lines!(f, c, dims, ring) end return n_on_line end function _burn_lines!( - B::AbstractArray, ::GI.AbstractCurveTrait, linestring, fill -) + f::F, c::C, dims::Tuple, ::GI.AbstractCurveTrait, linestring +) where {F<:Function,C<:Function} isfirst = true local firstpoint, laststop n_on_line = 0 @@ -46,19 +66,19 @@ function _burn_lines!( stop=(x=GI.x(point), y=GI.y(point)), ) laststop = line.stop - n_on_line += _burn_line!(B, line, fill) + n_on_line += _burn_line!(f, c, dims, line) end return n_on_line end function _burn_lines!( - B::AbstractArray, t::GI.LineTrait, line, fill -) + f::F, c::C, dims::Tuple, t::GI.LineTrait, line +) where {F<:Function,C<:Function} p1, p2 = GI.getpoint(t, line) line1 = ( start=(x=GI.x(p1), y=GI.y(p1)), stop=(x=GI.x(p2), y=GI.y(p2)), ) - return _burn_line!(B, line1, fill) + return _burn_line!(f, c, dims, line1) end # _burn_line! @@ -66,26 +86,22 @@ end # Line-burning algorithm # Burns a single line into a raster with value where pixels touch a line # +# Function `f` does the actual work when passed a Dimension Tuple of a pixel to burn, +# and `c` is an initialisation callback that is passed the maximyum +# number of times `f` will be called. It may be called less than that. +# # TODO: generalise to Irregular spans? -function _burn_line!(A::AbstractRaster, line, fill) - - xdim, ydim = dims(A, DEFAULT_POINT_ORDER) - regular = map((xdim, ydim)) do d - @assert (parent(lookup(d)) isa Array) - lookup(d) isa AbstractSampled && span(d) isa Regular - end - msg = """ - Can only fill lines where dimensions have `Regular` lookups. - Consider using `boundary=:center`, reprojecting the crs, - or make an issue in Rasters.jl on github if you need this to work. - """ - all(regular) || throw(ArgumentError(msg)) +function _burn_line!(f::Function, c::Function, dims::Tuple, line::NamedTuple) + xdim, ydim = DD.dims(dims, DEFAULT_POINT_ORDER) + di = DimIndices(dims) + @assert xdim isa XDim + @assert ydim isa YDim @assert order(xdim) == order(ydim) == Lookups.ForwardOrdered() @assert locus(xdim) == locus(ydim) == Lookups.Center() - raster_x_step = abs(step(span(A, X))) - raster_y_step = abs(step(span(A, Y))) + raster_x_step = abs(step(span(xdim))) + raster_y_step = abs(step(span(ydim))) raster_x_offset = @inbounds xdim[1] - raster_x_step / 2 # Shift from center to start of pixel raster_y_offset = @inbounds ydim[1] - raster_y_step / 2 @@ -117,15 +133,18 @@ function _burn_line!(A::AbstractRaster, line, fill) # Int starting points for the line. +1 converts to julia indexing j, i = trunc(Int, relstart.x) + 1, trunc(Int, relstart.y) + 1 # Int - # For arbitrary dimension indexing - dimconstructors = map(DD.basetypeof, (xdim, ydim)) + n_on_line = 0 + + # Pass of number of runs of `f` to callback `c` + # This can help with e.g. allocating a vector + c(manhattan_distance + 1) if manhattan_distance == 0 - D = map((d, o) -> d(o), dimconstructors, (j, i)) - if checkbounds(Bool, A, D...) - @inbounds A[D...] = fill + D = map(rebuild, dims, (j, i)) + if checkbounds(Bool, di, D...) + f(D) + n_on_line += 1 end - n_on_line = 1 return n_on_line end @@ -153,18 +172,17 @@ function _burn_line!(A::AbstractRaster, line, fill) n_on_line = 0 countx = county = 0 - # Int steps to move allong the line step_j = signbit(diff_x) * -2 + 1 step_i = signbit(diff_y) * -2 + 1 # Travel one grid cell at a time. Start at zero for the current cell for _ in 0:manhattan_distance - D = map((d, o) -> d(o), dimconstructors, (j, i)) - if checkbounds(Bool, A, D...) - @inbounds A[D...] = fill + D = map(rebuild, dims, (j, i)) + if checkbounds(Bool, di, D...) + f(D) + n_on_line += 1 end - # Only move in either X or Y coordinates, not both. if abs(max_x) < abs(max_y) max_x += delta_x @@ -176,12 +194,6 @@ function _burn_line!(A::AbstractRaster, line, fill) county +=1 end end - return n_on_line -end -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. - """ - throw(ArgumentError(msg)) + + return n_on_line end diff --git a/src/methods/burning/polygon.jl b/src/methods/burning/polygon.jl index f6d4512e4..1ce6c48a9 100644 --- a/src/methods/burning/polygon.jl +++ b/src/methods/burning/polygon.jl @@ -24,14 +24,14 @@ function _burn_polygon!(B::AbstractDimArray, trait, geom; # Lines n_on_line = 0 if boundary !== :center - _check_intervals(B, boundary) + _check_intervals(B, boundary, verbose) if boundary === :touches - if _check_intervals(B, boundary) + if _check_intervals(B, boundary, verbose) # Add line pixels n_on_line = _burn_lines!(B, geom; fill)::Int end elseif boundary === :inside - if _check_intervals(B, boundary) + if _check_intervals(B, boundary, verbose) # Remove line pixels n_on_line = _burn_lines!(B, geom; fill=!fill)::Int end @@ -135,9 +135,9 @@ end const INTERVALS_INFO = "makes more sense on `Intervals` than `Points` and will have more correct results. You can construct dimensions with a `X(values; sampling=Intervals(Center()))` to achieve this" -@noinline _check_intervals(B) = - _chki(B) ? true : (@info "burning lines $INTERVALS_INFO"; false) -@noinline _check_intervals(B, boundary) = - _chki(B) ? true : (@info "`boundary=:$boundary` $INTERVALS_INFO"; false) +@noinline _check_intervals(B, verbose::Bool) = + _chki(B) ? true : (verbose && @info "burning lines $INTERVALS_INFO"; false) +@noinline _check_intervals(B, boundary, verbose::Bool) = + _chki(B) ? true : (verbose && @info "`boundary=:$boundary` $INTERVALS_INFO"; false) _chki(B) = all(map(s -> s isa Intervals, sampling(dims(B)))) diff --git a/src/methods/extract.jl b/src/methods/extract.jl index feb0b6001..a34cbc92e 100644 --- a/src/methods/extract.jl +++ b/src/methods/extract.jl @@ -1,9 +1,93 @@ +# Object to hold extract keywords +struct Extractor{T,P,K,N<:NamedTuple{K},ID,G,I,S,F} + names::N + id::ID + geometry::G + index::I + skipmissing::S + flatten::F +end +Base.@constprop :aggressive @inline function Extractor(x, data; + name::NTuple{<:Any,Symbol}, + skipmissing, + flatten, + geometry, + id, + index, + kw... +) + nt = NamedTuple{name}(name) + id = _booltype(id) + g = _booltype(geometry) + i = _booltype(index) + sm = _booltype(skipmissing) + f = _booltype(flatten) + Extractor(x, data, nt, id, g, i, sm, f) +end +function Extractor( + x, data, nt::N, id::ID, g::G, i::I, sm::S, f::F +) where {N<:NamedTuple{K},ID,G,I,S,F} where K + P = _proptype(x; skipmissing=sm, names=nt) + T = _geom_rowtype(x, data; id, geometry=g, index=i, names=nt, skipmissing=sm) + Extractor{T,P,K,N,ID,G,I,S,F}(nt, id, g, i, sm, f) +end + +Base.eltype(::Extractor{T}) where T = T + +# This mutable object is passed into closures as +# fields are type-stable in clusores but variables are not +mutable struct LineRefs{T} + i::Int + prev::CartesianIndex{2} + rows::Vector{T} + function LineRefs{T}() where T + i = 1 + prev = CartesianIndex((typemin(Int), typemin(Int))) + new{T}(i, prev, Vector{T}()) + end +end + +function _initialise!(lr::LineRefs{T}) where T + lr.i = 1 + lr.prev = CartesianIndex((typemin(Int), typemin(Int))) + lr.rows = Vector{T}() +end + +_geom_rowtype(A, geom; kw...) = + _geom_rowtype(A, GI.trait(geom), geom; kw...) +_geom_rowtype(A, ::Nothing, geoms; kw...) = + _geom_rowtype(A, first(skipmissing(geoms)); kw...) +_geom_rowtype(A, ::GI.AbstractGeometryTrait, geom; kw...) = + _geom_rowtype(A, first(GI.getpoint(geom)); kw...) +_geom_rowtype(A, ::Nothing, geoms::Missing; kw...) = + _rowtype(A, missing; kw...) +function _geom_rowtype(A, ::GI.AbstractPointTrait, p; kw...) + tuplepoint = if GI.is3d(p) + (GI.x(p), GI.y(p), GI.z(p)) + else + (GI.x(p), GI.y(p)) + end + return _rowtype(A, tuplepoint; kw...) +end + +# skipinvalid: can G and I be missing. skipmissing: can nametypes be missing +function _proptype(x; + skipmissing, names::NamedTuple{K}, kw... +) where K + NamedTuple{K,Tuple{_nametypes(x, names, skipmissing)...}} +end + """ - extract(x, data; kw...) + extract(x, geometries; kw...) + +Extracts the value of `Raster` or `RasterStack` for the passed in geometries, +returning an `Vector{NamedTuple}` with properties for `:geometry` and `Raster` +or `RasterStack` layer values. + +For lines, linestrings and linear rings points are extracted for each pixel that +the line touches. -Extracts the value of `Raster` or `RasterStack` at given points, returning -an iterable of `NamedTuple` with properties for `:geometry` and raster or -stack layer values. +For polygons, all cells witih centers covered by the polygon are returned. Note that if objects have more dimensions than the length of the point tuples, sliced arrays or stacks will be returned instead of single values. @@ -15,13 +99,22 @@ $DATA_ARGUMENT # Keywords -- `geometry`: include `:geometry` in returned `NamedTuple`, `true` by default. -- `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. +- `geometry`: include a `:geometry` field in rows, which will be a + tuple point. Either the original point for points or the pixel + center point for line and polygon extract. `true` by default. +- `index`: include `:index` field of extracted points in rows, `false` by default. +- `name`: a `Symbol` or `Tuple` of `Symbol` corresponding to layer/s + of a `RasterStack` to extract. All layers are extracted by default. - `skipmissing`: skip missing points automatically. +- `flatten`: flatten extracted points from multiple + geometries into a single vector. `true` by default. + Unmixed point geometries are always flattened. + Flattening is slow and single threaded, `flatten=false` may be a + large performance improvement in combination with `threaded=true`. - `atol`: a tolerance for floating point lookup values for when the `Lookup` contains `Points`. `atol` is ignored for `Intervals`. --$GEOMETRYCOLUMN_KEYWORD +$BOUNDARY_KEYWORD +$GEOMETRYCOLUMN_KEYWORD # Example @@ -50,205 +143,353 @@ extract(st, obs; skipmissing=true) (geometry = (0.32, 40.24), bio1 = 16.321388f0, bio3 = 41.659454f0, bio5 = 30.029825f0, bio7 = 25.544561f0, bio12 = 480.0f0) ``` -Note: passing in arrays, geometry collections or feature collections +Note: passing in arrays, geometry collections or feature collections containing a mix of points and other geometries has undefined results. """ function extract end -@inline function extract(x::RasterStackOrArray, data; - names=_names(x), name=names, skipmissing=false, geometry=true, index=false, geometrycolumn=nothing, kw... +Base.@constprop :aggressive @inline function extract(x::RasterStackOrArray, data; + geometrycolumn=nothing, + names::NTuple{<:Any,Symbol}=_names(x), + name::NTuple{<:Any,Symbol}=names, + skipmissing=false, + flatten=true, + id=false, + geometry=true, + index=false, + kw... ) n = DD._astuple(name) - geoms = _get_geometries(data, geometrycolumn) - _extract(x, geoms; - dims=DD.dims(x, DEFAULT_POINT_ORDER), - names=NamedTuple{n}(n), - # These keywords are converted to _True/_False for type stability later on - # The @inline above helps constant propagation of the Bools - geometry=_booltype(geometry), - index=_booltype(index), - skipmissing=_booltype(skipmissing), - kw... - ) + g, g1 = if GI.isgeometry(data) + data, data + elseif GI.isfeature(data) + g = GI.geometry(data) + g, g + else + gs = _get_geometries(data, geometrycolumn) + gs, (isempty(gs) || all(ismissing, gs)) ? missing : first(Base.skipmissing(gs)) + end + + xp = _prepare_for_burning(x) + e = Extractor(xp, g1; name, skipmissing, flatten, id, geometry, index) + id_init = 1 + return _extract(xp, e, id_init, g; kw...) end -function _extract(A::RasterStackOrArray, geom::Missing, names, kw...) - T = _rowtype(A, geom; names, kw...) - [_maybe_add_fields(T, map(_ -> missing, names), missing, missing)] +# TODO use a GeometryOpsCore method like `applyreduce` here? +function _extract(A::RasterStackOrArray, e::Extractor{T}, id::Int, geom::Missing; kw...) where T + return if istrue(e.skipmissing) + T[] + else + T[_maybe_add_fields(T, map(_ -> missing, e.names), id, missing, missing)] + end end -function _extract(A::RasterStackOrArray, geom; names, kw...) - _extract(A, GI.geomtrait(geom), geom; names, kw...) +function _extract(A::RasterStackOrArray, e::Extractor, id::Int, geom; kw...) + _extract(A, e, id, GI.geomtrait(geom), geom; kw...) end -function _extract(A::RasterStackOrArray, ::Nothing, geoms; - names, skipmissing, kw... -) - T = _rowtype(A, eltype(geoms); names, skipmissing, kw...) +function _extract(A::RasterStackOrArray, e::Extractor{T}, id::Int, ::Nothing, geoms::AbstractArray; + threaded=false, progress=true, kw... +) where T # Handle empty / all missing cases - (length(geoms) > 0 && any(!ismissing, geoms)) || return T[] - - geom1 = first(Base.skipmissing(geoms)) + isempty(geoms) && return T[] + geom1 = all(ismissing, geoms) ? missing : first(Base.skipmissing(geoms)) trait1 = GI.trait(geom1) - # We need to split out points from other geoms - # TODO this will fail with mixed point/geom vectors + # We split out points from other geoms for performance if trait1 isa GI.PointTrait - rows = Vector{T}(undef, length(geoms)) - if istrue(skipmissing) - j = 1 - for i in eachindex(geoms) - g = geoms[i] - ismissing(g) && continue - e = _extract_point(T, A, g, skipmissing; names, kw...) - if !ismissing(e) - rows[j] = e - j += 1 + allpoints = true + rows = Vector{T}(undef, length(geoms)) + if istrue(e.skipmissing) + if threaded + nonmissing = Vector{Bool}(undef, length(geoms)) + _run(1:length(geoms), threaded, progress, "Extracting points...") do i + g = geoms[i] + GI.geomtrait(g) isa GI.PointTrait || return nothing + loc_id = id + i - 1 + nonmissing[i] = _extract_point!(rows, A, e, loc_id, g, i; kw...)::Bool + return nothing + end + # This could use less memory if we reuse `rows` and shorten it + rows = rows[nonmissing] + else + j_ref = Ref(1) + # For non-threaded be more memory efficient + _run(1:length(geoms), false, progress, "Extracting points...") do i + g = geoms[i] + loc_id = id + i - 1 + if GI.geomtrait(g) isa GI.PointTrait + j_ref[] += _extract_point!(rows, A, e, loc_id, g, j_ref[]; kw...)::Bool + end + return nothing end - nothing + deleteat!(rows, j_ref[]:length(rows)) end - deleteat!(rows, j:length(rows)) else - for i in eachindex(geoms) + _run(1:length(geoms), threaded, progress, "Extracting points...") do i g = geoms[i] - rows[i] = _extract_point(T, A, g, skipmissing; names, kw...)::T - nothing + loc_id = id + i - 1 + _extract_point!(rows, A, e, loc_id, g, i; kw...) + return nothing end end - return rows + # If we found a non-point geometry, + # ignore these rows and start again generically + allpoints && return rows + end + + row_vecs = Vector{Vector{T}}(undef, length(geoms)) + if threaded + thread_line_refs = [LineRefs{T}() for _ in 1:Threads.nthreads()] + _run(1:length(geoms), threaded, progress, "Extracting geometries...") do i + line_refs = thread_line_refs[Threads.threadid()] + loc_id = id + i - 1 + row_vecs[i] = _extract(A, e, loc_id, geoms[i]; line_refs, kw...) + end else - # This will be a list of vectors, we need to flatten it into one - rows = Iterators.flatten(_extract(A, g; names, skipmissing, kw...) for g in geoms) - if istrue(skipmissing) - return collect(_skip_missing_rows(rows, _missingval_or_missing(A), names)) - else - return collect(rows) + line_refs = LineRefs{T}() + _run(1:length(geoms), threaded, progress, "Extracting geometries...") do i + loc_id = id + i - 1 + row_vecs[i] = _extract(A, e, loc_id, geoms[i]; line_refs, kw...) end end + + # TODO this is a bit slow and only on one thread + if istrue(e.flatten) + n = sum(map(length, row_vecs)) + out_rows = Vector{T}(undef, n) + k::Int = 1 + for row_vec in row_vecs + for row in row_vec + out_rows[k] = row + k += 1 + end + end + return out_rows + else + return row_vecs + end end -function _extract(A::RasterStackOrArray, ::GI.AbstractFeatureTrait, feature; kw...) - _extract(A, GI.geometry(feature); kw...) -end -function _extract(A::RasterStackOrArray, ::GI.FeatureCollectionTrait, fc; kw...) - # Fall back to the Array/iterator method for feature collections - _extract(A, [GI.geometry(f) for f in GI.getfeature(fc)]; kw...) -end -function _extract(A::RasterStackOrArray, ::GI.AbstractMultiPointTrait, geom; - skipmissing, kw... -) - T = _rowtype(A, GI.getpoint(geom, 1); names, skipmissing, kw...) - rows = (_extract_point(T, A, p, skipmissing; kw...) for p in GI.getpoint(geom)) - return skipmissing isa _True ? collect(_skip_missing_rows(rows, _missingval_or_missing(A), names)) : collect(rows) +function _extract( + A::RasterStackOrArray, e::Extractor{T}, id::Int, ::GI.AbstractMultiPointTrait, geom; kw... +)::Vector{T} where T + n = GI.npoint(geom) + rows = _init_rows(e, n) + i = 1 + for p in GI.getpoint(geom) + i += _extract_point!(rows, A, e, id, p, i; kw...)::Bool + end + if istrue(skipmissing) + # Remove excees rows where missing + deleteat!(rows, i:length(rows)) + end + return rows end -function _extract(A::RasterStackOrArray, ::GI.PointTrait, geom; - skipmissing, kw... -) - T = _rowtype(A, geom; names, skipmissing, kw...) - _extract_point(T, A, geom, skipmissing; kw...) +function _extract(A::RasterStackOrArray, e::Extractor{T}, id::Int, ::GI.PointTrait, p; kw...) where T + rows = _init_rows(e, 1) + _extract_point!(rows, A, e, id, p, 1; kw...) end -function _extract(A::RasterStackOrArray, t::GI.AbstractGeometryTrait, geom; - names, skipmissing, kw... -) - # Make a raster mask of the geometry - template = view(A, Touches(GI.extent(geom))) - ods = otherdims(A, DEFAULT_POINT_ORDER) - if length(ods) > 0 - template = view(template, map(d -> rebuild(d, firstindex(d)), ods)) +@noinline function _extract( + A::RasterStackOrArray, e::Extractor{T}, id::Int, ::GI.AbstractLineStringTrait, geom; + line_refs=LineRefs{T}(), + kw... +)::Vector{T} where T + _initialise!(line_refs) + # Subset/offst is not really needed for line buring + # But without it we can get different fp errors + # to polygons and eng up with lines in different + # places when they are right on the line. + dims, offset = _template(A, geom) + dp = DimPoints(A) + function _length_callback(n) + resize!(line_refs.rows, n) end - p = first(GI.getpoint(geom)) - tuplepoint = if GI.is3d(geom) - GI.x(p), GI.y(p), GI.z(p) - else - GI.x(p), GI.y(p) + + _burn_lines!(_length_callback, dims, geom) do D + I = CartesianIndex(map(val, D)) + # Avoid duplicates from adjacent line segments + line_refs.prev == I && return nothing + line_refs.prev = I + # Make sure we only hit this pixel once + # D is always inbounds + line_refs.i += _maybe_set_row!(line_refs.rows, e.skipmissing, e, id, A, dp, offset, I, line_refs.i) + return nothing end - T = _rowtype(A, tuplepoint; names, skipmissing, kw...) - B = boolmask(geom; to=template, kw...) - offset = CartesianIndex(map(x -> first(x) - 1, parentindices(parent(template)))) + deleteat!(line_refs.rows, line_refs.i:length(line_refs.rows)) + return line_refs.rows +end +@noinline function _extract( + A::RasterStackOrArray, e::Extractor{T}, id::Int, ::GI.AbstractGeometryTrait, geom; kw... +)::Vector{T} where T + # Make a raster mask of the geometry + dims, offset = _template(A, geom) + B = boolmask(geom; to=dims, burncheck_metadata=NoMetadata(), kw...) + n = count(B) + dp = DimPoints(A) + i = 1 + rows = _init_rows(e, n) # Add a row for each pixel that is `true` in the mask - rows = (_missing_or_fields(T, A, Tuple(I + offset), names, skipmissing) for I in CartesianIndices(B) if B[I]) - # Maybe skip missing rows - if istrue(skipmissing) - return collect(Base.skipmissing(rows)) + for I in CartesianIndices(B) + B[I] || continue + i += _maybe_set_row!(rows, e.skipmissing, e, id, A, dp, offset, I, i) + end + # Cleanup + deleteat!(rows, i:length(rows)) + return rows +end + +function _template(x, geom) + ods = otherdims(x, DEFAULT_POINT_ORDER) + # Build a dummy raster in case this is a stack + # views are just easier on an array than dims + t1 = Raster(FillArrays.Zeros(size(x)), dims(x)) + t2 = if length(ods) > 0 + view(t1, map(d -> rebuild(d, firstindex(d)), ods)) else - return collect(rows) + t1 end + I = dims2indices(dims(t2), Touches(GI.extent(geom))) + t3 = view(t2, I...) + offset = CartesianIndex(map(i -> first(i) - 1, I)) + return dims(t3), offset end -_extract(::Type{T}, A::RasterStackOrArray, trait::GI.PointTrait, point; skipmissing, kw...) where T = - _extract_point(T, A, point, skipmissing; kw...) -function _missing_or_fields(::Type{T}, A, I, names, skipmissing) where T - props = _prop_nt(A, I, names) - if istrue(skipmissing) && _ismissingval(A, props) - missing +Base.@assume_effects :foldable function _maybe_set_row!( + rows, skipmissing::_True, e, id, A, dp, offset, I, i; +) + props = _prop_nt(e, A, I, skipmissing) + return if ismissing(props) + 0 else - geom = DimPoints(A)[I...] - _maybe_add_fields(T, props, geom, I) + _maybe_set_row!(rows, _False(), e, id, A, dp, offset, I, i; props) end end +Base.@assume_effects :foldable function _maybe_set_row!( + rows::Vector{T}, skipmissing::_False, e, id, A, dp, offset, I, i; + props=_prop_nt(e, A, I, skipmissing) +) where T + Ioff = I + offset + geom = dp[Ioff] + rows[i] = _maybe_add_fields(T, props, id, geom, Tuple(Ioff)) + return 1 +end + +_init_rows(e::Extractor{T}, n=0) where T = Vector{T}(undef, n) -_ismissingval(A::Union{Raster,RasterStack}, props) = +Base.@assume_effects :foldable _ismissingval(A::Union{Raster,RasterStack}, props)::Bool = _ismissingval(missingval(A), props) -_ismissingval(A::Union{Raster,RasterStack}, props::NamedTuple) = +Base.@assume_effects :foldable _ismissingval(A::Union{Raster,RasterStack}, props::NamedTuple)::Bool = _ismissingval(missingval(A), props) -_ismissingval(mvs::NamedTuple, props::NamedTuple{K}) where K = - any(k -> ismissing(props[k]) || props[k] === mvs[k], K) -_ismissingval(mv, props::NamedTuple) = any(x -> ismissing(x) || x === mv, props) -_ismissingval(mv, prop) = (mv === prop) +Base.@assume_effects :foldable _ismissingval(mvs::NamedTuple{K}, props::NamedTuple{K}) where K = + any(DD.unrolled_map((p, mv) -> ismissing(p) || p === mv, props, mvs))::Bool +Base.@assume_effects :foldable _ismissingval(mvs::NamedTuple, props::NamedTuple{K}) where K = + _ismissingval(mvs[K], props)::Bool +Base.@assume_effects :foldable _ismissingval(mv, props::NamedTuple)::Bool = + any(DD.unrolled_map(p -> ismissing(p) || p === mv, props)) +Base.@assume_effects :foldable _ismissingval(mv, prop)::Bool = (mv === prop) -@inline _prop_nt(st::AbstractRasterStack, I, ::NamedTuple{K}) where K = st[I...][K] -@inline _prop_nt(A::AbstractRaster, I, ::NamedTuple{K}) where K = NamedTuple{K}((A[I...],)) +# We always return NamedTuple +Base.@assume_effects :foldable function _prop_nt( + ::Extractor{<:Any,P,K}, st::AbstractRasterStack, I, sm::_False +)::P where {P,K} + P(st[K][I]) +end +Base.@assume_effects :foldable function _prop_nt( + ::Extractor{<:Any,P,K}, st::AbstractRasterStack{K}, I, sm::_False +)::P where {P,K} + P(st[I]) +end +Base.@assume_effects :foldable function _prop_nt( + ::Extractor{<:Any,P,K}, A::AbstractRaster, I, sm::_False +)::P where {P,K} + P(NamedTuple{K}((A[I],))) +end +Base.@assume_effects :foldable function _prop_nt( + ::Extractor{<:Any,P,K}, st::AbstractRasterStack, I, sm::_True +)::Union{P,Missing} where {P,K} + x = st[K][I] + _ismissingval(A, x) ? missing : x::P +end +Base.@assume_effects :foldable function _prop_nt( + ::Extractor{<:Any,P,K}, st::AbstractRasterStack{K}, I, sm::_True +)::Union{P,Missing} where {P,K} + x = st[I] + _ismissingval(st, x) ? missing : x::P +end +Base.@assume_effects :foldable function _prop_nt( + ::Extractor{<:Any,P,K}, A::AbstractRaster, I, sm::_True +)::Union{P,Missing} where {P,K} + x = A[I] + _ismissingval(A, x) ? missing : NamedTuple{K}((x,))::P +end # Extract a single point -@inline function _extract_point(::Type, x::RasterStackOrArray, point::Missing, skipmissing::_True; kw...) - return missing +# Missing point to remove +@inline function _extract_point!(rows::Vector{T}, x::RasterStackOrArray, e::Extractor, id, point, i; kw...) where T + _extract_point!(rows, x, e, id, e.skipmissing, point, i; kw...) end -@inline function _extract_point(::Type{T}, x::RasterStackOrArray, point::Missing, skipmissing::_False; - names, kw... -) where T - props = map(_ -> missing, names) - geom = missing - I = missing - return _maybe_add_fields(T, props, geom, I) +@inline function _extract_point!( + rows::Vector{T}, x::RasterStackOrArray, e::Extractor{T,<:Any,K}, id, skipmissing::_True, point::Missing, i; + kw... +) where {T,K} + return false end -@inline function _extract_point(::Type{T}, x::RasterStackOrArray, point, skipmissing::_True; - dims, names::NamedTuple{K}, atol=nothing, kw... +# Normal point with missing / out of bounds data removed +@inline function _extract_point!( + rows::Vector{T}, x::RasterStackOrArray, e::Extractor{T,<:Any,K}, id, skipmissing::_True, point, i; + dims=DD.dims(x, DEFAULT_POINT_ORDER), + atol=nothing, + kw... ) where {T,K} # Get the actual dimensions available in the object - coords = map(DD.commondims(x, dims)) do d + coords = map(DD.dims(x, dims)) do d _dimcoord(d, point) end # If any are coordinates missing, also return missing for everything - if any(map(ismissing, coords)) - return missing - else + if !any(map(ismissing, coords)) selector_dims = map(dims, coords) do d, c _at_or_contains(d, c, atol) end selectors = map(val, DD.dims(selector_dims, DD.dims(x))) # Check the selector is in bounds / actually selectable I = DD.selectindices(DD.dims(x, dims), selectors; err=_False())::Union{Nothing,Tuple{Int,Int}} - if isnothing(I) - return missing - else + if !isnothing(I) D = map(rebuild, selector_dims, I) - if x isa Raster + if x isa Raster prop = x[D] - _ismissingval(missingval(x), prop) && return missing + _ismissingval(missingval(x), prop) && return false props = NamedTuple{K,Tuple{eltype(x)}}((prop,)) else props = x[D][K] - _ismissingval(missingval(x), props) && return missing + _ismissingval(missingval(x), props) && return false end - return _maybe_add_fields(T, props, point, I) + rows[i] = _maybe_add_fields(T, props, id, coords, I) + return true end end + return false end -@inline function _extract_point(::Type{T}, x::RasterStackOrArray, point, skipmissing::_False; - dims, names::NamedTuple{K}, atol=nothing, kw... -) where {T,K} +# Missing point to keep +@inline function _extract_point!( + rows::Vector{T}, x::RasterStackOrArray, e::Extractor{T,P,K}, id, skipmissing::_False, point::Missing, i; kw... +) where {T,P,K} + props = map(_ -> missing, e.names) + geom = missing + I = missing + rows[i] = _maybe_add_fields(T, props, id, geom, I) + return true +end +# Normal point with missing / out of bounds data kept with `missing` fields +@inline function _extract_point!( + rows::Vector{T}, x::RasterStackOrArray, e::Extractor{T,P,K}, id, skipmissing::_False, point, i; + dims=DD.dims(x, DEFAULT_POINT_ORDER), + atol=nothing, + kw... +) where {T,P,K} # Get the actual dimensions available in the object - coords = map(DD.commondims(x, dims)) do d + coords = map(dims) do d _dimcoord(d, point) end # If any are coordinates missing, also return missing for everything - if any(map(ismissing, coords)) - return _maybe_add_fields(T, map(_ -> missing, names), missing, missing) + rows[i] = if any(map(ismissing, coords)) + _maybe_add_fields(T, map(_ -> missing, names), id, missing, missing) else selector_dims = map(dims, coords) do d, c _at_or_contains(d, c, atol) @@ -257,40 +498,27 @@ end # Check the selector is in bounds / actually selectable I = DD.selectindices(DD.dims(x, dims), selectors; err=_False())::Union{Nothing,Tuple{Int,Int}} if isnothing(I) - return _maybe_add_fields(T, map(_ -> missing, names), point, missing) + _maybe_add_fields(T, map(_ -> missing, e.names), id, coords, missing) else D = map(rebuild, selector_dims, I) - props = if x isa Raster - NamedTuple{K,Tuple{eltype(x)}}((x[D],)) + props = if x isa Raster + P(NamedTuple{K,Tuple{eltype(x)}}((x[D],))) else - NamedTuple(x[D])[K] + P(NamedTuple(x[D])[K]) end - return _maybe_add_fields(T, props, point, I) + _maybe_add_fields(T, props, id, coords, I) end end + return 1 end # Maybe add optional fields -Base.@assume_effects :total function _maybe_add_fields(::Type{T}, props::NamedTuple, point, I)::T where {T<:NamedTuple{K}} where K - if :geometry in K - :index in K ? merge((; geometry=point, index=I), props) : merge((; geometry=point), props) - else - :index in K ? merge((; index=I), props) : props - end |> T -end - -@inline _skip_missing_rows(rows, ::Missing, names) = - Iterators.filter(row -> !any(ismissing, row), rows) -@inline _skip_missing_rows(rows, missingval, names) = - Iterators.filter(row -> !any(x -> ismissing(x) || x === missingval, row), rows) -@inline function _skip_missing_rows(rows, missingval::NamedTuple, names::NamedTuple{K}) where K - # first check if all fields are equal - if so just call with the first value - if Base.allequal(missingval) == 1 - return _skip_missing_rows(rows, first(missingval), names) - else - Iterators.filter(rows) do row - # rows may or may not contain a :geometry field, so map over keys instead - !any(key -> ismissing(row[key]) || row[key] === missingval[key], K) - end - end +# It is critically important for performance that this is type stable +Base.@assume_effects :total function _maybe_add_fields( + ::Type{T}, props::NamedTuple, id, point::Union{Tuple,Missing}, I +)::T where {T<:NamedTuple{K}} where K + row = :index in K ? merge((; index=I), props) : props + row = :geometry in K ? merge((; geometry=point), row) : row + row = :id in K ? merge((; id), row) : row + return T(row) end diff --git a/src/methods/mask.jl b/src/methods/mask.jl index 1401a990c..fc8d58da8 100644 --- a/src/methods/mask.jl +++ b/src/methods/mask.jl @@ -35,7 +35,8 @@ $SUFFIX_KEYWORD These can be used when `with` is a GeoInterface.jl compatible object: -$SHAPE_KEYWORDS +$SHAPE_KEYWORD +$BOUNDARY_KEYWORD $GEOMETRYCOLUMN_KEYWORD # Example @@ -332,7 +333,8 @@ function boolmask!(dest::AbstractRaster, data; lock=nothing, progress=true, threaded=false, - allocs=_burning_allocs(dest; threaded), + burncheck_metadata=Metadata(), + allocs=_burning_allocs(dest; threaded, burncheck_metadata), geometrycolumn=nothing, collapse=nokw, kw... diff --git a/src/methods/mosaic.jl b/src/methods/mosaic.jl index 2e9a3f5df..a823db4d2 100644 --- a/src/methods/mosaic.jl +++ b/src/methods/mosaic.jl @@ -1,3 +1,6 @@ +const RasterVecOrTuple = Union{Tuple{Vararg{AbstractRaster}},AbstractArray{<:AbstractRaster}} +const RasterStackVecOrTuple = Union{Tuple{Vararg{AbstractRasterStack}},AbstractArray{<:AbstractRasterStack}} + """ mosaic(f, regions...; missingval, atol) mosaic(f, regions; missingval, atol) @@ -28,42 +31,29 @@ If your mosaic has has apparent line errors, increase the `atol` value. Here we cut out Australia and Africa from a stack, and join them with `mosaic`. ```jldoctest -using Rasters, RasterDataSources, ArchGDAL, Plots -st = RasterStack(WorldClim{Climate}; month=1); - -africa = st[X(-20.0 .. 60.0), Y(-40.0 .. 35.0)] -a = plot(africa) - -aus = st[X(100.0 .. 160.0), Y(-50.0 .. -10.0)] -b = plot(aus) - -# Combine with mosaic -mos = mosaic(first, aus, africa) -c = plot(mos) +using Rasters, RasterDataSources, NaturalEarth, DataFrames +countries = naturalearth("admin_0_countries", 110) |> DataFrame +climate = RasterStack(WorldClim{Climate}, (:tmin, :tmax, :prec, :wind); month=July) +country_climates = map(("Norway", "Denmark", "Sweden")) do name + country = subset(countries, :NAME => ByRow(==("Norway"))) + trim(mask(climate; with=country); pad=10) +end +scandinavia_climate = trim(mosaic(first, country_climates)) +plot(scandinavia_climate) -savefig(a, "build/mosaic_example_africa.png") -savefig(b, "build/mosaic_example_aus.png") -savefig(c, "build/mosaic_example_combined.png") -nothing +savefig("build/mosaic_example_combined.png") # output ``` -### Individual continents - -![arica](mosaic_example_africa.png) - -![aus](mosaic_example_aus.png) - -### Mosaic of continents +### Mosaic of countries ![mosaic](mosaic_example_combined.png) $EXPERIMENTAL """ -function mosaic(f::Function, r1::RasterStackOrArray, rs::RasterStackOrArray...; kw...) +mosaic(f::Function, r1::RasterStackOrArray, rs::RasterStackOrArray...; kw...) = mosaic(f, (r1, rs...); kw...) -end mosaic(f::Function, regions; kw...) = _mosaic(f, first(regions), regions; kw...) function _mosaic(f::Function, A1::AbstractRaster, regions; missingval=nokw, @@ -74,23 +64,26 @@ function _mosaic(f::Function, A1::AbstractRaster, regions; force=false, kw... ) - isnothing(missingval) && throw(ArgumentError("missingval cannot be `nothing` for `mosaic`")) - missingval = if isnokw(missingval) + V = Vector{promote_type(map(Missings.nonmissingtype ∘ eltype, regions)...)} + T = Base.promote_op(f, V) + dims = _mosaic(Tuple(map(DD.dims, regions))) + l1 = first(regions) + + missingval = if isnothing(missingval) + throw(ArgumentError("missingval cannot be `nothing` for `mosaic`")) + elseif isnokw(missingval) mv = Rasters.missingval(first(regions)) isnokwornothing(mv) ? missing : mv else missingval end - missingval_pair = if !isnothing(filename) && (ismissing(missingval) || isnokwornothing(missingval)) - _type_missingval(eltype(A1)) => missing - elseif missingval isa Pair + missingval_pair = if missingval isa Pair missingval + elseif !isnothing(filename) && (ismissing(missingval) || isnokw(missingval)) + _type_missingval(eltype(A1)) => missing else missingval => missingval end - T = Base.promote_type(typeof(last(missingval_pair)), Base.promote_eltype(regions...)) - dims = _mosaic(Tuple(map(DD.dims, regions))) - l1 = first(regions) return create(filename, T, dims; name=name(l1), @@ -100,10 +93,10 @@ function _mosaic(f::Function, A1::AbstractRaster, regions; options, force ) do C - _mosaic!(f, C, regions; missingval, kw...) + mosaic!(f, C, regions; missingval, kw...) end end -function _mosaic(f::Function, ::AbstractRasterStack, regions; +function _mosaic(f::Function, r1::AbstractRasterStack, regions; filename=nothing, suffix=keys(first(regions)), kw... @@ -112,14 +105,14 @@ function _mosaic(f::Function, ::AbstractRasterStack, regions; layers = map(suffix, map(values, regions)...) do s, A... mosaic(f, A...; filename, suffix=s, kw...) end - return DD.rebuild_from_arrays(first(regions), Tuple(layers)) + return DD.rebuild_from_arrays(r1, layers) end """ mosaic!(f, x, regions...; missingval, atol) mosaic!(f, x, regions::Tuple; missingval, atol) -Combine `regions` in `x` using the function `f`. +Combine `regions` in `Raster` or `RasterStack` `x` using the function `f`. # Arguments @@ -143,18 +136,25 @@ Combine `regions` in `x` using the function `f`. # Example -Cut out Australia and Africa stacks, then combined them -into a single stack. +Cut out scandinavian countries and plot: ```jldoctest -using Rasters, RasterDataSources, ArchGDAL, Statistics, Plots -st = read(RasterStack(WorldClim{Climate}; month=1)) -aus = st[X=100.0 .. 160.0, Y=-50.0 .. -10.0] -africa = st[X=-20.0 .. 60.0, Y=-40.0 .. 35.0] -mosaic!(first, st, aus, africa) -plot(st) +using Rasters, RasterDataSources, NaturalEarth, DataFrames +# Get climate data form worldclim +climate = RasterStack(WorldClim{Climate}, (:tmin, :tmax, :prec, :wind); month=July) +# And country borders from natural earth +countries = naturalearth("admin_0_countries", 110) |> DataFrame +# Cut out each country +country_climates = map(("Norway", "Denmark", "Sweden")) do name + country = subset(countries, :NAME => ByRow(==(name))) + trim(mask(climate; with=country); pad=10) +end +# Mosaic together to a single raster +scandinavia_climate = mosaic(first, country_climates) +# And plot +plot(scandinavia_climate) + savefig("build/mosaic_bang_example.png") -nothing # output ``` @@ -164,49 +164,149 @@ nothing $EXPERIMENTAL """ mosaic!(f::Function, dest::RasterStackOrArray, regions::RasterStackOrArray...; kw...) = - _mosaic!(f, dest, regions; kw...) - -function _mosaic!(f::Function, A::AbstractRaster{T}, regions::Union{Tuple,AbstractArray}; + mosaic!(f, dest, regions; kw...) +function mosaic!( + f::Function, + dest::RasterStackOrArray, + regions::Union{Tuple,AbstractArray}; + op=_reduce_op(f, missingval(dest)), + kw... +) + # Centering avoids pixel edge floating point error + dest_centered = _prepare_for_burning(dest; order=nothing) + regions_centered = map(r -> _prepare_for_burning(r; order=nothing), regions) + _mosaic!(f, op, dest_centered, regions_centered; kw...) + return dest +end +function _mosaic!( + f::typeof(mean), op::Nothing, dest::Raster, regions::RasterVecOrTuple; + kw... +) + if length(regions) <= typemax(UInt8) + _mosaic_mean!(dest, UInt8, regions; kw...) + elseif length(regions) <= typemax(UInt16) + _mosaic_mean!(dest, UInt16, regions; kw...) + else + _mosaic_mean!(dest, UInt32, regions; kw...) + end +end +function _mosaic!( + f::typeof(length), op::Nothing, dest::AbstractRaster, regions::RasterVecOrTuple; + kw... +) + for region in regions + _count_region!(dest, region; kw...) + end + return dest +end +# Where there is a known reduction operator we can apply each region as a whole +function _mosaic!( + f::Function, op, dest::AbstractRaster, regions::RasterVecOrTuple; + kw... +) + for region in regions + _mosaic_region!(op, dest, region; kw...) + end + return dest +end +# Generic unknown functions +function _mosaic!( + f::Function, op::Nothing, A::AbstractRaster{T}, regions::RasterVecOrTuple; missingval=missingval(A), atol=nothing ) where T isnokwornothing(missingval) && throw(ArgumentError("destination array must have a `missingval`")) + R = promote_type(map(Missings.nonmissingtype ∘ eltype, regions)...) + buffer = Vector{R}(undef, length(regions)) _without_mapped_crs(A) do A1 broadcast!(A1, DimSelectors(A1; atol)) do ds # Get all the regions that have this point - ls = foldl(regions; init=()) do acc, l - if DD.hasselection(l, ds) - v = l[ds...] - (acc..., l) - else - acc - end - end - values = foldl(ls; init=()) do acc, l - v = l[ds...] - if isnothing(Rasters.missingval(l)) - (acc..., v) - elseif ismissing(Rasters.missingval(l)) - ismissing(v) ? acc : (acc..., v) - else - v === Rasters.missingval(l) ? acc : (acc..., v) + i = 0 + for r in regions + if DD.hasselection(r, ds) + x = r[ds...] + if x !== Rasters.missingval(r) + i += 1 + buffer[i] = x + end end end - if length(values) === 0 + if i === 0 missingval else - f(values) + f(view(buffer, 1:i)) end end end return A end -function _mosaic!(f::Function, st::AbstractRasterStack, regions::Union{Tuple,AbstractArray}; kw...) +function _mosaic!( + f::Function, + op, + st::AbstractRasterStack, + regions::RasterStackVecOrTuple; + kw... +) map(values(st), map(values, regions)...) do A, r... - mosaic!(f, A, r...; kw...) + _mosaic!(f, op, A, r; kw...) end return st end +function _mosaic_mean!(dest, ::Type{T}, regions; kw...) where T + # Note: sum and count are separate broadcasts because + # most disk formats don't support writing a tuple + + # Define a Raster to count into + counts = create(nothing, T, dest; missingval=zero(T)) + counts .= zero(T) + for region in regions + # Add region to dest + _mosaic_region!(Base.add_sum, dest, region; kw...) + # Count region + _count_region!(counts, region; kw...) + end + # Divide dest by counts + # Avoid divide by zero for missing values + dest .= ((d, c) -> d === missingval(dest) ? missingval(dest) : d / c).(dest, counts) + return dest +end + +function _mosaic_region!(op, dest, region; atol=nothing, kw...) + function skip_or_op(a, b) + if b === missingval(region) + a + elseif a === missingval(dest) + b + else + op(a, b) + end + end + ext = _maybe_pad_floats(extent(region), sampling(dest)) + selectors = map(sampling(dest)) do sa + ispoints(sa) ? At(; atol) : Contains() + end + ds = DimSelectors(view(dest, ext); selectors) + # `parent` needed to skip broadcast checks + dest[ext] .= skip_or_op.(parent(view(dest, ext)), parent(view(region, ds))) + return dest +end + +function _count_region!(count::AbstractRaster{T}, region::AbstractRaster; kw...) where T + function skip_or_count(a, b) + if b === missingval(region) + a + elseif a === missingval(count) + oneunit(Missings.nonmissingtype(T)) + else + a + oneunit(a) + end + end + ext = extent(region) + ds = DimSelectors(view(count, ext)) + view(count, ext) .= skip_or_count.(view(count, ext), view(region, ds)) + return count +end + _mosaic(alldims::Tuple{<:DimTuple,Vararg{DimTuple}}) = map(_mosaic, alldims...) function _mosaic(dims::Dimension...) map(dims) do d @@ -248,7 +348,6 @@ function _mosaic(span::Regular, lookup::AbstractSampled, lookups::LookupTuple) end return rebuild(lookup; data=newindex) end - function _mosaic(::Irregular, lookup::AbstractSampled, lookups::LookupTuple) newindex = sort(union(map(parent, lookups)...); order=LA.ordering(order(lookup))) return rebuild(lookup; data=newindex) @@ -264,3 +363,15 @@ function _mosaic(span::Explicit, lookup::AbstractSampled, lookups::LookupTuple) newbounds = vcat(permutedims(newlower), permutedims(newupper)) return rebuild(lookup; data=newindex, span=Explicit(newbounds)) end + +# Pad floats for intervals so that small floating point +# error doesn't exclude values in nealy matching lookups +function _maybe_pad_floats(ext::Extent{K}, sampling::Tuple) where K + map(values(Extents.bounds(ext)), sampling) do b, sa + if isintervals(sa) && eltype(first(b)) <: AbstractFloat + b[1] - 10eps(b[1]), b[2] + 10eps(b[2]) + else + b + end + end |> Extent{K} +end \ No newline at end of file diff --git a/src/methods/rasterize.jl b/src/methods/rasterize.jl index 2d9e278ec..1948a8e64 100644 --- a/src/methods/rasterize.jl +++ b/src/methods/rasterize.jl @@ -1,4 +1,4 @@ -struct _TakeFirst{MV} +struct _TakeFirst{MV} <: Function missingval::MV end (tf::_TakeFirst)(a, b) = a === tf.missingval ? b : a @@ -316,15 +316,10 @@ const RASTERIZE_KEYWORDS = """ - `op`: A reducing function that accepts two values and returns one, like `min` to `minimum`. For common methods this will be assigned for you, or is not required. But you can use it instead of a `reducer` as it will usually be faster. -- `shape`: force `data` to be treated as `:polygon`, `:line` or `:point`, where possible - Points can't be treated as lines or polygons, and lines may not work as polygons, but - an attempt will be made. -- `geometrycolumn`: `Symbol` to manually select the column the geometries are in - when `data` is a Tables.jl compatible table, or a tuple of `Symbol` for columns of - point coordinates. -- `progress`: show a progress bar, `true` by default, `false` to hide.. -- `verbose`: print information and warnings when there are problems with the rasterisation. - `true` by default. +$GEOM_KEYWORDS +$GEOMETRYCOLUMN_KEYWORD +$PROGRESS_KEYWORD +$VERBOSE_KEYWORD $THREADED_KEYWORD - `threadsafe`: specify that custom `reducer` and/or `op` functions are thread-safe, in that the order of operation or blocking does not matter. For example, @@ -360,8 +355,6 @@ $RASTERIZE_ARGUMENTS These are detected automatically from `data` where possible. -$GEOMETRYCOLUMN_KEYWORD -$GEOM_KEYWORDS $RASTERIZE_KEYWORDS $FILENAME_KEYWORD $SUFFIX_KEYWORD diff --git a/src/methods/reproject.jl b/src/methods/reproject.jl index e18910bb6..baad1cc04 100644 --- a/src/methods/reproject.jl +++ b/src/methods/reproject.jl @@ -1,4 +1,3 @@ - """ reproject(obj; crs) @@ -61,17 +60,14 @@ function _reproject(source::GeoFormat, target::GeoFormat, dim::Union{XDim,YDim}, # This is a dumb way to do this. But it save having to inspect crs, # and prevents reprojections that don't make sense from working. # A better method for this should be implemented in future. - return first(reproject(source, target, dim, [val])) + return first(_reproject(source, target, dim, [val])) end function _reproject(source::GeoFormat, target::GeoFormat, dim, vals::NTuple{N}) where N - reps = reproject(source, target, dim, [vals...]) + reps = _reproject(source, target, dim, [vals...]) return ntuple(x -> reps[x], N) end function _reproject(source::GeoFormat, target::GeoFormat, dim::Union{XDim,YDim}, vals::AbstractArray) - reshape(reproject(source, target, dim, vec(vals)), size(vals)) -end -function _reproject(source::Nothing, target::GeoFormat, dim::Union{XDim,YDim}, vals::AbstractArray) - reshape(reproject(source, target, dim, vec(vals)), size(vals)) + reshape(_reproject(source, target, dim, vec(vals)), size(vals)) end function _reproject(source::GeoFormat, target::GeoFormat, dim::Union{XDim,YDim}, vals::AbstractVector) throw_extension_error(reproject, "Proj", :RastersProjExt, (source, target, dim, vals)) diff --git a/src/methods/shared_docstrings.jl b/src/methods/shared_docstrings.jl index 88bb5437d..0790f50ab 100644 --- a/src/methods/shared_docstrings.jl +++ b/src/methods/shared_docstrings.jl @@ -30,15 +30,15 @@ const CRS_KEYWORD = """ - `crs`: a `crs` which will be attached to the resulting raster when `to` not passed or is an `Extent`. Otherwise the crs from `to` is used. """ - -const SHAPE_KEYWORDS = """ -- `shape`: Force `data` to be treated as `:polygon`, `:line` or `:point` geometries. - using points or lines as polygons may have unexpected results. +const BOUNDARY_KEYWORD = """ - `boundary`: for polygons, include pixels where the `:center` is inside the polygon, where the polygon `:touches` the pixel, or that are completely `:inside` the polygon. The default is `:center`. """ - +const SHAPE_KEYWORD = """ +- `shape`: Force `data` to be treated as `:polygon`, `:line` or `:point` geometries. + using points or lines as polygons may have unexpected results. +""" const THREADED_KEYWORD = """ - `threaded`: run operations in parallel, `false` by default. In some circumstances `threaded` can give large speedups over single-threaded operation. This can be true for complicated @@ -53,7 +53,8 @@ $TO_KEYWORD $RES_KEYWORD $SIZE_KEYWORD $CRS_KEYWORD -$SHAPE_KEYWORDS +$BOUNDARY_KEYWORD +$SHAPE_KEYWORD """ const DATA_ARGUMENT = """ diff --git a/src/methods/zonal.jl b/src/methods/zonal.jl index c08b62cc5..41fb4c5fb 100644 --- a/src/methods/zonal.jl +++ b/src/methods/zonal.jl @@ -72,16 +72,17 @@ insertcols!(january_stats, 1, :country => first.(split.(countries.ADMIN, r"[^A-Z 3 columns and 243 rows omitted ``` """ -zonal(f, x::RasterStackOrArray; of, kw...) = _zonal(f, x, of; kw...) +zonal(f, x::RasterStackOrArray; of, skipmissing=true, kw...) = + _zonal(f, _prepare_for_burning(x), of; skipmissing, kw...) _zonal(f, x::RasterStackOrArray, of::RasterStackOrArray; kw...) = _zonal(f, x, Extents.extent(of); kw...) _zonal(f, x::RasterStackOrArray, of::DimTuple; kw...) = _zonal(f, x, Extents.extent(of); kw...) # We don't need to `mask` with an extent, it's square so `crop` will do enough. -_zonal(f, x::Raster, of::Extents.Extent; skipmissing=true) = +_zonal(f, x::Raster, of::Extents.Extent; skipmissing) = _maybe_skipmissing_call(f, crop(x; to=of, touches=true), skipmissing) -function _zonal(f, x::RasterStack, ext::Extents.Extent; skipmissing=true) +function _zonal(f, x::RasterStack, ext::Extents.Extent; skipmissing) cropped = crop(x; to=ext, touches=true) length(cropped) > 0 || return missing return maplayers(cropped) do A @@ -96,7 +97,7 @@ _zonal(f, x, ::GI.AbstractFeatureCollectionTrait, fc; kw...) = _zonal(f, x::RasterStackOrArray, ::GI.AbstractFeatureTrait, feature; kw...) = _zonal(f, x, GI.geometry(feature); kw...) function _zonal(f, x::AbstractRaster, ::GI.AbstractGeometryTrait, geom; - skipmissing=true, kw... + skipmissing, kw... ) cropped = crop(x; to=geom, touches=true) length(cropped) > 0 || return missing @@ -104,7 +105,7 @@ function _zonal(f, x::AbstractRaster, ::GI.AbstractGeometryTrait, geom; return _maybe_skipmissing_call(f, masked, skipmissing) end function _zonal(f, st::AbstractRasterStack, ::GI.AbstractGeometryTrait, geom; - skipmissing=true, kw... + skipmissing, kw... ) cropped = crop(st; to=geom, touches=true) length(cropped) > 0 || return map(_ -> missing, st) diff --git a/src/sources/grd.jl b/src/sources/grd.jl index b7f7556ba..9f03e9148 100644 --- a/src/sources/grd.jl +++ b/src/sources/grd.jl @@ -218,7 +218,7 @@ function Base.write(filename::String, ::GRDsource, A::AbstractRaster; if write _mmapgrd(filename, source_eltype(mod), size(A); write=true) do M - f(rebuild(A, _maybe_modify(M, mod))) + f(_maybe_modify(M, mod)) end end diff --git a/src/utils.jl b/src/utils.jl index 41c11f609..84878461c 100644 --- a/src/utils.jl +++ b/src/utils.jl @@ -284,7 +284,6 @@ end # to distinguish between objects returned by _get_geometries and other objects struct IterableOfGeometries end - # Chunking # NoKW means true @@ -434,7 +433,7 @@ end # Other -_progress(args...; kw...) = ProgressMeter.Progress(args...; color=:blue, barlen=50, kw...) +_progress(args...; kw...) = ProgressMeter.Progress(args...; dt=0.1, color=:blue, barlen=50, kw...) # Function barrier for splatted vector broadcast @noinline _do_broadcast!(f, x, args...) = broadcast!(f, x, args...) @@ -525,47 +524,48 @@ _names(A::AbstractRaster) = (Symbol(name(A)),) _names(A::AbstractRasterStack) = keys(A) using DimensionalData.Lookups: _True, _False + +_booltype(::_True) = _True() +_booltype(::_False) = _False() _booltype(x) = x ? _True() : _False() + istrue(::_True) = true istrue(::_False) = false +istrue(::Type{_True}) = true +istrue(::Type{_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... + id=_False(), 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) + keys = _rowkeys(id, geometry, index, names) + types = _rowtypes(x, _G, _I, id, 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} +@generated function _rowtypes( + x, ::Type{G}, ::Type{I}, id, geometry, index, 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)...} + ts = Expr(:tuple) + istrue(id) && push!(ts.args, Int) + if istrue(skipmissing) + istrue(geometry) && push!(ts.args, G) + istrue(index) && push!(ts.args, I) + else + istrue(geometry) && push!(ts.args, Union{Missing,G}) + istrue(index) && push!(ts.args, Union{Missing,I}) + end + :(Tuple{$ts...,_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},) +@inline _nametypes(::Raster{T}, ::NamedTuple{Names}, sm::_True) where {T,Names} = + (nonmissingtype(T),) +@inline _nametypes(::Raster{T}, ::NamedTuple{Names}, sm::_False) where {T,Names} = + (Union{Missing,T},) # This only compiles away when generated @generated function _nametypes( ::RasterStack{<:Any,T}, ::NamedTuple{PropNames}, skipmissing::_True @@ -580,7 +580,11 @@ end 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...) +_rowkeys(id::_True, geometry::_False, index::_False, names::NamedTuple{Names}) where Names = (:id, Names...) +_rowkeys(id::_True, geometry::_True, index::_False, names::NamedTuple{Names}) where Names = (:id, :geometry, Names...) +_rowkeys(id::_True, geometry::_True, index::_True, names::NamedTuple{Names}) where Names = (:id, :geometry, :index, Names...) +_rowkeys(id::_True, geometry::_False, index::_True, names::NamedTuple{Names}) where Names = (:id, :index, Names...) +_rowkeys(id::_False, geometry::_False, index::_False, names::NamedTuple{Names}) where Names = Names +_rowkeys(id::_False, geometry::_True, index::_False, names::NamedTuple{Names}) where Names = (:geometry, Names...) +_rowkeys(id::_False, geometry::_True, index::_True, names::NamedTuple{Names}) where Names = (:geometry, :index, Names...) +_rowkeys(id::_False, geometry::_False, index::_True, names::NamedTuple{Names}) where Names = (:index, Names...) \ No newline at end of file diff --git a/test/extract.jl b/test/extract.jl new file mode 100644 index 000000000..9247c3754 --- /dev/null +++ b/test/extract.jl @@ -0,0 +1,506 @@ +using Rasters, Test, DataFrames, Extents +import GeoInterface as GI +using Rasters.Lookups, Rasters.Dimensions + +include(joinpath(dirname(pathof(Rasters)), "../test/test_utils.jl")) + +dimz = (X(9.0:1.0:10.0), Y(0.1:0.1:0.2)) +rast = Raster(Union{Int,Missing}[1 2; 3 4], dimz; name=:test, missingval=missing) +rast2 = Raster([5 6; 7 8], dimz; name=:test2, missingval=5) +rast_m = Raster([1 2; 3 missing], dimz; name=:test, missingval=missing) +st = RasterStack(rast, rast2) + +pts = [missing, (9.0, 0.1), (9.0, 0.2), (10.0, 0.3), (10.0, 0.2)] +poly = GI.Polygon([[(8.0, 0.0), (11.0, 0.0), (11.0, 0.4), (8.0, 0.0)]]) +linestring = GI.LineString([(8.0, 0.0), (9.5, 0.0), (10.0, 0.4)]) +line = GI.Line([(8.0, 0.0), (12.0, 0.4)]) +table = (geometry=pts, foo=zeros(4)) + +@testset "Points" begin + @testset "From Raster" begin + @testset "skipmissing=false" begin + ex = extract(rast, pts; skipmissing=false) + T = @NamedTuple{geometry::Union{Missing,Tuple{Float64,Float64}},test::Union{Missing,Int64}} + @test eltype(ex) == T + @test all(ex .=== T[ + (geometry = missing, test = missing) + (geometry = (9.0, 0.1), test=1) + (geometry = (9.0, 0.2), test=2) + (geometry = (10.0, 0.3), test=missing) + (geometry = (10.0, 0.2), test=4) + ]) + end + @testset "skipmissing=true" begin + ex = extract(rast_m, pts; skipmissing=true); + T = @NamedTuple{geometry::Tuple{Float64, Float64}, test::Int64} + @test eltype(ex) == T + @test all(ex .=== T[(geometry = (9.0, 0.1), test = 1), (geometry = (9.0, 0.2), test = 2)]) + end + + @testset "skipmissing=true, geometry=false" begin + ex = extract(rast_m, pts; skipmissing=true, geometry=false) + T = @NamedTuple{test::Int64} + @test eltype(ex) == T + @test all(ex .=== T[(test = 1,), (test = 2,)]) + @test all(extract(rast_m, pts; skipmissing=true, geometry=false, index=true) .=== [ + (index = (1, 1), test = 1,) + (index = (1, 2), test = 2,) + ]) + end + @testset "reverse points" begin + # NamedTuple (reversed) points - tests a Table that iterates over points + T = @NamedTuple{geometry::Union{Missing,Tuple{Float64,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 = (9.0, 0.1), test = 1) + (geometry = (10.0, 0.2), test = 4) + (geometry = (10.0, 0.3), test = missing) + ]) + end + @testset "Vector points" begin + @test extract(rast, [[9.0, 0.1], [10.0, 0.2]]) == [ + (geometry = (9.0, 0.1), test = 1) + (geometry = (10.0, 0.2), test = 4) + ] + end + end + + @testset "From RasterStack" begin + T = @NamedTuple{geometry::Union{Missing,Tuple{Float64,Float64}},test::Union{Missing,Int64},test2::Union{Missing,Int64}} + @test all(extract(st, [missing, (9.0, 0.1), (10.0, 0.2), (10.0, 0.3)]) .=== T[ + (geometry = missing, test = missing, test2 = missing) + (geometry = (9.0, 0.1), test = 1, test2 = 5) + (geometry = (10.0, 0.2), test = 4, test2 = 8) + (geometry = (10.0, 0.3), test = missing, test2 = missing) + ]) + @test extract(st, [missing, (9.0, 0.1), (10.0, 0.2), (10.0, 0.3)]; skipmissing=true) == [ + (geometry = (10.0, 0.2), test = 4, test2 = 8) + ] + @test extract(st, [missing, (9.0, 0.1), (10.0, 0.2), (10.0, 0.3)]; skipmissing=true, geometry=false) == [ + (test = 4, test2 = 8) + ] + T = @NamedTuple{index::Union{Missing, Tuple{Int,Int}}, test::Union{Missing, Int64}, test2::Union{Missing, Int64}} + @test extract(st, [missing, (9.0, 0.1), (10.0, 0.2), (10.0, 0.3)]; skipmissing=true, geometry=false, index=true) == T[ + (index = (2, 2), test = 4, test2 = 8) + ] + # Subset with `names` + T = @NamedTuple{geometry::Union{Missing, Tuple{Float64, Float64}}, test2::Union{Missing, Int64}} + @test all(extract(st, [missing, (9.0, 0.1), (10.0, 0.2), (10.0, 0.3)]; name=(:test2,)) .=== T[ + (geometry = missing, test2 = missing) + (geometry = (9.0, 0.1), test2 = 5) + (geometry = (10.0, 0.2), test2 = 8) + (geometry = (10.0, 0.3), test2 = missing) + ]) + # Subset with `names` and `skipmissing` with mixed missingvals + @test extract(st, [missing, (9.0, 0.1), (10.0, 0.2), (10.0, 0.3)]; name=(:test2,), skipmissing=true) == [ + (geometry = (10.0, 0.2), test2 = 8) + ] + @test extract(st, [missing, (9.0, 0.1), (10.0, 0.2), (10.0, 0.3)]; name=(:test,), skipmissing=true) == [ + (geometry = (9.0, 0.1), test = 1) + (geometry = (10.0, 0.2), test = 4) + ] + end + + @testset "Errors" begin + # Missing coord errors + @test_throws ArgumentError extract(rast, [(0.0, missing), (9.0, 0.1), (9.0, 0.2), (10.0, 0.3)]) + @test_throws ArgumentError extract(rast, [(9.0, 0.1), (0.0, missing), (9.0, 0.2), (10.0, 0.3)]) + @test_throws ArgumentError extract(rast, [(X=0.0, Y=missing), (9.0, 0.1), (9.0, 0.2), (10.0, 0.3)]) + end +end + +@testset "Polygons" begin + # Extract a polygon + T = @NamedTuple{geometry::Union{Missing,Tuple{Float64,Float64}},test::Union{Missing,Int64}} + @test all(extract(rast_m, poly) .=== T[ + (geometry = (9.0, 0.1), test = 1) + (geometry = (10.0, 0.1), test = 3) + (geometry = (10.0, 0.2), test = missing) + ]) + # Test all the keyword combinations + @test all(extract(rast_m, poly) .=== T[ + (geometry = (9.0, 0.1), test = 1) + (geometry = (10.0, 0.1), test = 3) + (geometry = (10.0, 0.2), test = missing) + ]) + T = @NamedTuple{test::Union{Missing,Int64}} + @test all(extract(rast_m, poly; geometry=false) .=== T[ + (test = 1,) + (test = 3,) + (test = missing,) + ]) + T = @NamedTuple{index::Union{Missing,Tuple{Int,Int}},test::Union{Missing,Int64}} + @test all(extract(rast_m, poly; geometry=false, index=true) .=== T[ + (index = (1, 1), test = 1) + (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}} + @test all(extract(rast_m, poly; index=true) .=== T[ + (geometry = (9.0, 0.1), index = (1, 1), test = 1) + (geometry = (10.0, 0.1), index = (2, 1), test = 3) + (geometry = (10.0, 0.2), index = (2, 2), test = missing) + ]) + @test extract(rast_m, poly; skipmissing=true) == [ + (geometry = (9.0, 0.1), test = 1) + (geometry = (10.0, 0.1), test = 3) + ] + @test extract(rast_m, poly; skipmissing=true, geometry=false) == [ + (test = 1,) + (test = 3,) + ] + @test extract(rast_m, poly; skipmissing=true, geometry=false, index=true) == [ + (index = (1, 1), test = 1) + (index = (2, 1), test = 3) + ] + @test extract(rast_m, poly; skipmissing=true, index=true) == [ + (geometry = (9.0, 0.1), index = (1, 1), test = 1) + (geometry = (10.0, 0.1), index = (2, 1), test = 3) + ] + @test extract(rast2, poly; skipmissing=true) == [ + (geometry = (10.0, 0.1), test2 = 7) + (geometry = (10.0, 0.2), test2 = 8) + ] + T = @NamedTuple{geometry::Union{Missing,Tuple{Float64,Float64}},test::Union{Missing,Int64}} + @test all(extract(rast_m, poly) .=== T[ + (geometry = (9.0, 0.1), test = 1) + (geometry = (10.0, 0.1), test = 3) + (geometry = (10.0, 0.2), test = missing) + ]) + @test extract(rast_m, poly; skipmissing=true, geometry=false, id=true) == [ + (id=1, test = 1,) + (id=1, test = 3,) + ] + @test extract(rast_m, [poly, poly]; skipmissing=true, geometry=false, id=true) == [ + (id=1, test = 1,) + (id=1, test = 3,) + (id=2, test = 1,) + (id=2, test = 3,) + ] + @test extract(rast_m, [poly, poly]; skipmissing=true, geometry=false, id=true, flatten=false) == [ + [(id=1, test = 1,), (id=1, test = 3,)], + [(id=2, test = 1,), (id=2, test = 3,)], + ] + + @testset "Vector of polygons" begin + ex = extract(rast_m, [poly, poly]) + @test eltype(ex) == T + @test all(ex .=== T[ + (geometry = (9.0, 0.1), test = 1) + (geometry = (10.0, 0.1), test = 3) + (geometry = (10.0, 0.2), test = missing) + (geometry = (9.0, 0.1), test = 1) + (geometry = (10.0, 0.1), test = 3) + (geometry = (10.0, 0.2), test = missing) + ]) + end +end + +@testset "Extract a linestring" begin + T = @NamedTuple{geometry::Union{Missing,Tuple{Float64,Float64}},test::Union{Missing,Int64}} + Tsm = @NamedTuple{geometry::Tuple{Float64,Float64},test::Int64} + linestrings = [linestring, linestring, linestring] + fc = GI.FeatureCollection(map(GI.Feature, linestrings)) + + # Single LineString + @test extract(rast, linestring) isa Vector{T} + @test all(extract(rast_m, linestring) .=== T[ + (geometry = (9.0, 0.1), test = 1) + (geometry = (10.0, 0.1), test = 3) + (geometry = (10.0, 0.2), test = missing) + ]) + @test all(extract(rast_m, linestring; skipmissing=true) .=== Tsm[ + (geometry = (9.0, 0.1), test = 1) + (geometry = (10.0, 0.1), test = 3) + ]) + + # Multiple linstrings + # Test all combinations of skipmissing, flatten and threaded + @test extract(rast_m, linestrings; skipmissing=true) isa Vector{Tsm} + @test extract(rast_m, fc; skipmissing=true) == + extract(rast_m, fc; skipmissing=true, threaded=true) == + extract(rast_m, linestrings; skipmissing=true) == + extract(rast_m, linestrings; skipmissing=true, threaded=true) == Tsm[ + (geometry = (9.0, 0.1), test = 1) + (geometry = (10.0, 0.1), test = 3) + (geometry = (9.0, 0.1), test = 1) + (geometry = (10.0, 0.1), test = 3) + (geometry = (9.0, 0.1), test = 1) + (geometry = (10.0, 0.1), test = 3) + ] + + @test extract(rast_m, linestrings; skipmissing=false) isa Vector{T} + @test all( + extract(rast_m, fc; skipmissing=false) .=== + extract(rast_m, fc; skipmissing=false, threaded=true) .=== + extract(rast_m, linestrings; skipmissing=false) .=== + extract(rast_m, linestrings; skipmissing=false, threaded=true) .=== T[ + (geometry = (9.0, 0.1), test = 1) + (geometry = (10.0, 0.1), test = 3) + (geometry = (10.0, 0.2), test = missing) + (geometry = (9.0, 0.1), test = 1) + (geometry = (10.0, 0.1), test = 3) + (geometry = (10.0, 0.2), test = missing) + (geometry = (9.0, 0.1), test = 1) + (geometry = (10.0, 0.1), test = 3) + (geometry = (10.0, 0.2), test = missing) + ]) + + @test extract(rast_m, linestrings; skipmissing=true, flatten=false) isa Vector{Vector{Tsm}} + @test extract(rast_m, linestrings; skipmissing=true, flatten=false) == + extract(rast_m, linestrings; skipmissing=true, flatten=false, threaded=true) == Vector{Tsm}[ + [(geometry = (9.0, 0.1), test = 1), (geometry = (10.0, 0.1), test = 3)], + [(geometry = (9.0, 0.1), test = 1), (geometry = (10.0, 0.1), test = 3)], + [(geometry = (9.0, 0.1), test = 1), (geometry = (10.0, 0.1), test = 3)], + ] + + # Nested Vector holding missing needs special handling to check equality + @test extract(rast_m, linestrings; skipmissing=false, flatten=false) isa Vector{Vector{T}} + ref = Vector{T}[ + [(geometry = (9.0, 0.1), test = 1), (geometry = (10.0, 0.1), test = 3), (geometry = (10.0, 0.2), test = missing)] , + [(geometry = (9.0, 0.1), test = 1), (geometry = (10.0, 0.1), test = 3), (geometry = (10.0, 0.2), test = missing)], + [(geometry = (9.0, 0.1), test = 1), (geometry = (10.0, 0.1), test = 3), (geometry = (10.0, 0.2), test = missing)], + ] + matching(a, b) = all(map(===, a, b)) + @test all(map(matching, extract(rast_m, linestrings; skipmissing=false, flatten=false, threaded=false), ref)) + @test all(map(matching, extract(rast_m, linestrings; skipmissing=false, flatten=false, threaded=true), ref)) +end + +@testset "Extract a line" begin + T = @NamedTuple{geometry::Union{Missing,Tuple{Float64,Float64}},test::Union{Missing,Int64}} + Tsm = @NamedTuple{geometry::Tuple{Float64,Float64},test::Int64} + lines = [line, line] + fc = GI.FeatureCollection(map(GI.Feature, lines)) + + # Single LineString + @test extract(rast, line) isa Vector{T} + @test all(extract(rast_m, line) .=== T[ + (geometry = (9.0, 0.1), test = 1) + (geometry = (9.0, 0.2), test = 2) + (geometry = (10.0, 0.2), test = missing) + ]) + @test all(extract(rast_m, line; skipmissing=true) .=== Tsm[ + (geometry = (9.0, 0.1), test = 1) + (geometry = (9.0, 0.2), test = 2) + ]) + + # Multiple linstrings + # Test all combinations of skipmissing, flatten and threaded + @test extract(rast_m, lines; skipmissing=true) isa Vector{Tsm} + @test extract(rast_m, fc; skipmissing=true) == + extract(rast_m, fc; skipmissing=true, threaded=true) == + extract(rast_m, lines; skipmissing=true) == + extract(rast_m, lines; skipmissing=true, threaded=true) == Tsm[ + (geometry = (9.0, 0.1), test = 1) + (geometry = (9.0, 0.2), test = 2) + (geometry = (9.0, 0.1), test = 1) + (geometry = (9.0, 0.2), test = 2) + ] + + @test extract(rast_m, lines; skipmissing=false) isa Vector{T} + @test all( + extract(rast_m, fc; skipmissing=false) .=== + extract(rast_m, fc; skipmissing=false, threaded=true) .=== + extract(rast_m, lines; skipmissing=false) .=== + extract(rast_m, lines; skipmissing=false, threaded=true) .=== T[ + (geometry = (9.0, 0.1), test = 1) + (geometry = (9.0, 0.2), test = 2) + (geometry = (10.0, 0.2), test = missing) + (geometry = (9.0, 0.1), test = 1) + (geometry = (9.0, 0.2), test = 2) + (geometry = (10.0, 0.2), test = missing) + ]) + + Tsm_i = @NamedTuple{id::Int,geometry::Tuple{Float64,Float64},test::Int64} + @test extract(rast, lines; skipmissing=true, id=true) isa Vector{Tsm_i} + @test all(extract(rast_m, fc; skipmissing=true, id=true) .=== + extract(rast_m, fc; skipmissing=true, threaded=true, id=true) .=== + extract(rast_m, lines; skipmissing=true, id=true) .=== + extract(rast_m, lines; skipmissing=true, threaded=true, id=true) .=== Tsm_i[ + (id=1, geometry = (9.0, 0.1), test = 1) + (id=1, geometry = (9.0, 0.2), test = 2) + (id=2, geometry = (9.0, 0.1), test = 1) + (id=2, geometry = (9.0, 0.2), test = 2) + ]) + + @test extract(rast_m, lines; skipmissing=true, flatten=false) isa Vector{Vector{Tsm}} + @test extract(rast_m, lines; skipmissing=true, flatten=false) == + extract(rast_m, lines; skipmissing=true, flatten=false, threaded=true) == Vector{Tsm}[ + [(geometry = (9.0, 0.1), test = 1), (geometry = (9.0, 0.2), test = 2)], + [(geometry = (9.0, 0.1), test = 1), (geometry = (9.0, 0.2), test = 2)], + ] + + # Nested Vector holding missing needs special handling to check equality + @test extract(rast_m, lines; skipmissing=false, flatten=false) isa Vector{Vector{T}} + ref = Vector{T}[ + [(geometry = (9.0, 0.1), test = 1), (geometry = (9.0, 0.2), test = 2), (geometry = (10.0, 0.2), test = missing)] , + [(geometry = (9.0, 0.1), test = 1), (geometry = (9.0, 0.2), test = 2), (geometry = (10.0, 0.2), test = missing)], + ] + matching(a, b) = all(map(===, a, b)) + @test all(map(matching, extract(rast_m, lines; skipmissing=false, flatten=false, threaded=false), ref)) + @test all(map(matching, extract(rast_m, lines; skipmissing=false, flatten=false, threaded=true), ref)) +end + +@testset "Table" begin + T = @NamedTuple{geometry::Union{Missing,Tuple{Float64, Float64}}, test::Union{Missing, Int64}} + @test all(extract(rast, table) .=== T[ + (geometry = missing, test = missing) + (geometry = (9.0, 0.1), test = 1) + (geometry = (9.0, 0.2), test = 2) + (geometry = (10.0, 0.3), test = missing) + (geometry = (10.0, 0.2), test = 4) + ]) + @test extract(rast, table; skipmissing=true) == [ + (geometry = (9.0, 0.1), test = 1) + (geometry = (9.0, 0.2), test = 2) + (geometry = (10.0, 0.2), test = 4) + ] + @test extract(rast, table; skipmissing=true, geometry=false) == [ + (test = 1,) + (test = 2,) + (test = 4,) + ] + @test extract(rast, table; skipmissing=true, geometry=false, index=true) == [ + (index = (1, 1), test = 1,) + (index = (1, 2), test = 2,) + (index = (2, 2), test = 4,) + ] + @test extract(rast, table; skipmissing=true, geometry=false, id=true) == [ + (id=2, test = 1,) + (id=3, test = 2,) + (id=5, test = 4,) + ] + T = @NamedTuple{id::Int, test::Union{Missing, Int64}} + @test all(extract(rast, table; skipmissing=false, geometry=false, id=true) .=== T[ + (id=1, test = missing,) + (id=2, test = 1,) + (id=3, test = 2,) + (id=4, test = missing,) + (id=5, test = 4,) + ]) + @test_throws ArgumentError extract(rast, (foo = zeros(4),)) +end + +@testset "Empty geoms" begin + @test all(extract(rast, []) .=== @NamedTuple{geometry::Missing,test::Union{Missing,Int}}[]) + @test all(extract(rast, []; geometry=false) .=== @NamedTuple{test::Tuple{Missing}}[]) + @test all(extract(rast, [missing, missing]; geometry=false) .=== @NamedTuple{test::Union{Missing,Int}}[ + (test = missing,) + (test = missing,) + ]) + @test typeof(extract(rast, []; geometry=false, skipmissing=true)) == Vector{@NamedTuple{test::Int}} + @test typeof(extract(rast, [missing, missing]; geometry=false, skipmissing=true)) == Vector{@NamedTuple{test::Int}} +end + +#= +Some benchmark and plotting code + +using BenchmarkTools +using ProfileView +using Chairmarks +using Cthulhu +using Rasters +using DataFrames +import GeoInterface as GI +using Rasters.Lookups, Rasters.Dimensions + +dimz = (X(5.0:0.1:15.0; sampling=Intervals(Start())), Y(-0.1:0.01:0.5; sampling=Intervals(Start()))) +rast_m = Raster(rand([missing, 1, 2], dimz); name=:test, missingval=missing) +rast = Raster(rand(Int, dimz); name=:test2, missingval=5) +st = RasterStack((a=rast, b=rast, c=Float32.(rast))) +ks = ntuple(x -> gensym(), 100) +st100 = RasterStack(NamedTuple{ks}(ntuple(x -> rast, length(ks)))) + +poly = GI.Polygon([[(8.0, 0.0), (11.0, 0.0), (11.0, 0.4), (8.0, 0.0)]]) +linestring = GI.LineString([(8.0, 0.0), (11.0, 0.0), (11.0, 0.4), (8.0, 0.0)]) +line = GI.Line([(8.0, 0.0), (12.0, 4.0)]) +polys = [poly for i in 1:10000] +lines = [line for i in 1:100000] +linestrings = [linestring for i in 1:100000] + +extract(rast, linestring) +extract(rast, polys; threaded=false) +@time extract(rast, polys; geometry=false, threaded=true, flatten=true) |> length +@time extract(rast, lines; geometry=false, threaded=true, flatten=true) |> length +@time extract(st, polys; geometry=false, threaded=true, flatten=true) |> length +@time extract(st, lines; geometry=false, threaded=false, flatten=true) |> length + +@time extract(rast_m, lines; geometry=false, threaded=true, flatten=false) |> length +@time extract(rast, lines; geometry=false, threaded=true, flatten=false) |> length +@time extract(rast_m, lines; skipmissing=true, geometry=false, threaded=true, flatten=false) |> length +@time extract(rast, lines; skipmissing=true, geometry=false, threaded=true, flatten=false) |> length +@time extract(st, lines; skipmissing=true, geometry=false, threaded=true, flatten=false) |> length +@time extract(st100, lines; geometry=false, threaded=false, flatten=false) |> length +@time extract(st100, lines; geometry=false, threaded=true, flatten=false) |> length +@time extract(rast, lines; geometry=false, threaded=true, flatten=true) |> length +@time extract(rast, lines; threaded=true, flatten=true, progress=false) |> length +@time extract(rast, lines; threaded=false, flatten=false, progress=false) |> length +@benchmark extract(rast, lines; threaded=false, flatten=true) +@benchmark extract(rast, lines; threaded=true, flatten=true) +@benchmark extract(rast, lines; threaded=false, flatten=false) +@benchmark extract(rast, lines; threaded=true, flatten=false) +@time extract(rast, polys; flatten=false); +extract(rast, polys; flatten=true); +extract(st100, linestring) +extract(st100, poly) + +@profview extract(rast, lines) +@profview extract(rast, lines; threaded=true) +@profview extract(rast, linestrings) +@profview extract(rast_m, linestrings) +@profview extract(rast, polys) +@profview extract(rast, polys; flatten=false) +@profview extract(st, linestring) +@profview extract(st, poly) + +# Profile running one function a lot. +# People will do this +f(rast, ls, n; skipmissing=true) = for _ in 1:n extract(rast, ls; skipmissing) end +@profview f(rast, polies, 10; skipmissing=false) +@profview f(rast, poly, 10000; skipmissing=false) +@profview f(rast_m, poly, 10000; skipmissing=false) +@profview f(rast, linestring, 10000; skipmissing=false) +@profview f(rast_m, linestring, 10000; skipmissing=false) +@profview f(st, linestring, 10000; skipmissing=false) +@profview f(rast, linestring, 10000) +@profview f(st, linestring, 100000; skipmissing=false) +@profview f(st, linestring, 100000; skipmissing=true) +@profview f(st100, linestring, 1000) +@profview f(st100, line, 1000) +@profview f(st100, poly, 10000; skipmissing=true) + +@profview f(st100, linestring, 1000; skipmissing=true) +@profview f(st100, linestring, 1000; skipmissing=false) +@profview f(st100, linestring, 1000; names=skipmissing=false) +keys(st100)[1:50] + + + +# Demo plots +using GLMakie +using GeoInterfaceMakie + +Makie.plot(rast; colormap=:Reds) +Makie.plot!(rast[Touches(GI.extent(first(polys)))] ; colormap=:Greens) +Makie.plot!(polys; alpha=0.2) +ps = getindex.(extract(rast, poly; index=true, flatten=true), :geometry) +Makie.scatter!(ps; color=:yellow) + +Makie.plot(rast; colormap=:Reds) +Makie.plot!(rast[Touches(GI.extent(first(polys)))] ; colormap=:Greens) +Makie.plot!(poly; alpha=0.7) +ps = getindex.(extract(rast, poly; index=true, flatten=true, boundary=:touches), :geometry); +Makie.scatter!(ps; color=:pink) + +Makie.plot(rast; colormap=:Reds) +Makie.plot!(rast[Touches(GI.extent(first(polys)))] ; colormap=:Greens) +Makie.plot!(polys; alpha=0.7) +ps = getindex.(extract(rast, poly; index=true, flatten=true, boundary=:inside), :geometry) +Makie.scatter!(ps; color=:pink) + +Makie.plot(rast; colormap=:Reds) +Makie.plot!(rast[Touches(GI.extent(first(polys)))] ; colormap=:Greens) +Makie.plot!(linestring; color=:violet, linewidth=5) +ps = getindex.(extract(rast, linestring; index=true), :geometry); +Makie.scatter!(ps; color=:yellow) + +=# diff --git a/test/methods.jl b/test/methods.jl index c2e15124f..5ebfe1fc2 100644 --- a/test/methods.jl +++ b/test/methods.jl @@ -17,6 +17,8 @@ pointvec = [(-20.0, 30.0), (0.0, 10.0), (0.0, 30.0), (-20.0, 30.0)] + + vals = [1, 2, 3, 4, 5] polygon = ArchGDAL.createpolygon(pointvec) multi_polygon = ArchGDAL.createmultipolygon([[pointvec]]) @@ -409,190 +411,6 @@ end createpoint(args...) = ArchGDAL.createpoint(args...) -@testset "extract" begin - dimz = (X(9.0:1.0:10.0), Y(0.1:0.1:0.2)) - rast = Raster(Union{Int,Missing}[1 2; 3 4], dimz; name=:test, missingval=missing) - rast2 = Raster([5 6; 7 8], dimz; name=:test2, missingval=5) - rast_m = Raster([1 2; 3 missing], dimz; name=:test, missingval=missing) - mypoints = [missing, (9.0, 0.1), (9.0, 0.2), (10.0, 0.3), (10.0, 0.2)] - table = (geometry=mypoints, foo=zeros(4)) - st = RasterStack(rast, rast2) - @testset "from Raster" begin - # Tuple points - ex = extract(rast, mypoints) - T = @NamedTuple{geometry::Union{Missing,Tuple{Float64,Float64}},test::Union{Missing,Int64}} - @test eltype(ex) == T - @test all(ex .=== T[ - (geometry = missing, test = missing) - (geometry = (9.0, 0.1), test=1) - (geometry = (9.0, 0.2), test=2) - (geometry = (10.0, 0.3), test=missing) - (geometry = (10.0, 0.2), test=4) - ]) - ex = extract(rast_m, mypoints; skipmissing=true) - T = @NamedTuple{geometry::Tuple{Float64, Float64}, test::Int64} - @test eltype(ex) == T - @test all(ex .=== T[(geometry = (9.0, 0.1), test = 1), (geometry = (9.0, 0.2), test = 2)]) - ex = extract(rast_m, mypoints; skipmissing=true, geometry=false) - T = @NamedTuple{test::Int64} - @test eltype(ex) == T - @test all(ex .=== T[(test = 1,), (test = 2,)]) - @test all(extract(rast_m, mypoints; skipmissing=true, geometry=false, index=true) .=== [ - (index = (1, 1), test = 1,) - (index = (1, 2), test = 2,) - ]) - # NamedTuple (reversed) points - tests a Table that iterates over points - 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) - (geometry = (Y = 0.3, X = 10.0), test = missing) - ]) - # Vector points - @test all(extract(rast, [[9.0, 0.1], [10.0, 0.2]]) .== [ - (geometry = [9.0, 0.1], test = 1) - (geometry = [10.0, 0.2], test = 4) - ]) - # 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{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) - (geometry = (10.0, 0.2), test = missing) - ]) - # Extract a vector of polygons - ex = extract(rast_m, [p, p]) - @test eltype(ex) == T - @test all(ex .=== T[ - (geometry = (9.0, 0.1), test = 1) - (geometry = (10.0, 0.1), test = 3) - (geometry = (10.0, 0.2), test = missing) - (geometry = (9.0, 0.1), test = 1) - (geometry = (10.0, 0.1), test = 3) - (geometry = (10.0, 0.2), test = missing) - ]) - # Test all the keyword combinations - @test all(extract(rast_m, p) .=== T[ - (geometry = (9.0, 0.1), test = 1) - (geometry = (10.0, 0.1), test = 3) - (geometry = (10.0, 0.2), test = missing) - ]) - T = @NamedTuple{test::Union{Missing,Int64}} - @test all(extract(rast_m, p; geometry=false) .=== T[ - (test = 1,) - (test = 3,) - (test = missing,) - ]) - T = @NamedTuple{index::Union{Missing,Tuple{Int,Int}},test::Union{Missing,Int64}} - @test all(extract(rast_m, p; geometry=false, index=true) .=== T[ - (index = (1, 1), test = 1) - (index = (2, 1), test = 3) - (index = (2, 2), test = missing) - ]) - 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) - (geometry = (10.0, 0.2), index = (2, 2), test = missing) - ]) - @test extract(rast_m, p; skipmissing=true) == [ - (geometry = (9.0, 0.1), test = 1) - (geometry = (10.0, 0.1), test = 3) - ] - @test extract(rast_m, p; skipmissing=true, geometry=false) == [ - (test = 1,) - (test = 3,) - ] - @test extract(rast_m, p; skipmissing=true, geometry=false, index=true) == [ - (index = (1, 1), test = 1) - (index = (2, 1), test = 3) - ] - @test extract(rast_m, p; skipmissing=true, index=true) == [ - (geometry = (9.0, 0.1), index = (1, 1), test = 1) - (geometry = (10.0, 0.1), index = (2, 1), test = 3) - ] - @test extract(rast2, p; skipmissing=true) == [ - (geometry = (10.0, 0.1), test2 = 7) - (geometry = (10.0, 0.2), test2 = 8) - ] - # Empty geoms - @test extract(rast, []) == NamedTuple{(:geometry, :test),Tuple{Missing,Missing}}[] - @test extract(rast, []; geometry=false) == NamedTuple{(:test,),Tuple{Missing}}[] - # Missing coord errors - @test_throws ArgumentError extract(rast, [(0.0, missing), (9.0, 0.1), (9.0, 0.2), (10.0, 0.3)]) - @test_throws ArgumentError extract(rast, [(9.0, 0.1), (0.0, missing), (9.0, 0.2), (10.0, 0.3)]) - @test_throws ArgumentError extract(rast, [(X=0.0, Y=missing), (9.0, 0.1), (9.0, 0.2), (10.0, 0.3)]) - end - - @testset "with table" begin - T = @NamedTuple{geometry::Union{Missing, Tuple{Float64, Float64}}, test::Union{Missing, Int64}} - @test all(extract(rast, table) .=== T[ - (geometry = missing, test = missing) - (geometry = (9.0, 0.1), test = 1) - (geometry = (9.0, 0.2), test = 2) - (geometry = (10.0, 0.3), test = missing) - (geometry = (10.0, 0.2), test = 4) - ]) - @test extract(rast, table; skipmissing=true) == [ - (geometry = (9.0, 0.1), test = 1) - (geometry = (9.0, 0.2), test = 2) - (geometry = (10.0, 0.2), test = 4) - ] - @test extract(rast, table; skipmissing=true, geometry=false) == [ - (test = 1,) - (test = 2,) - (test = 4,) - ] - @test extract(rast, table; skipmissing=true, geometry=false, index=true) == [ - (index = (1, 1), test = 1,) - (index = (1, 2), test = 2,) - (index = (2, 2), test = 4,) - ] - - @test_throws ArgumentError extract(rast, (foo = zeros(4),)) - end - - @testset "from stack" begin - T = @NamedTuple{geometry::Union{Missing,Tuple{Float64,Float64}},test::Union{Missing,Int64},test2::Union{Missing,Int64}} - @test all(extract(st, [missing, (9.0, 0.1), (10.0, 0.2), (10.0, 0.3)]) .=== T[ - (geometry = missing, test = missing, test2 = missing) - (geometry = (9.0, 0.1), test = 1, test2 = 5) - (geometry = (10.0, 0.2), test = 4, test2 = 8) - (geometry = (10.0, 0.3), test = missing, test2 = missing) - ]) - @test extract(st, [missing, (9.0, 0.1), (10.0, 0.2), (10.0, 0.3)]; skipmissing=true) == [ - (geometry = (10.0, 0.2), test = 4, test2 = 8) - ] - @test extract(st2, [missing, (2, 2), (2, 1)]; skipmissing=true) == [ - (geometry = (2, 1), a = 7.0, b = 2.0) - ] - @test extract(st, [missing, (9.0, 0.1), (10.0, 0.2), (10.0, 0.3)]; skipmissing=true, geometry=false) == [ - (test = 4, test2 = 8) - ] - T = @NamedTuple{index::Union{Missing, Tuple{Int,Int}}, test::Union{Missing, Int64}, test2::Union{Missing, Int64}} - @test extract(st, [missing, (9.0, 0.1), (10.0, 0.2), (10.0, 0.3)]; skipmissing=true, geometry=false, index=true) == T[ - (index = (2, 2), test = 4, test2 = 8) - ] - # Subset with `names` - T = @NamedTuple{geometry::Union{Missing, Tuple{Float64, Float64}}, test2::Union{Missing, Int64}} - @test all(extract(st, [missing, (9.0, 0.1), (10.0, 0.2), (10.0, 0.3)]; name=(:test2,)) .=== T[ - (geometry = missing, test2 = missing) - (geometry = (9.0, 0.1), test2 = 5) - (geometry = (10.0, 0.2), test2 = 8) - (geometry = (10.0, 0.3), test2 = missing) - ]) - # Subset with `names` and `skipmissing` with mixed missingvals - @test extract(st, [missing, (9.0, 0.1), (10.0, 0.2), (10.0, 0.3)]; name=(:test2,), skipmissing=true) == [ - (geometry = (10.0, 0.2), test2 = 8) - ] - @test extract(st, [missing, (9.0, 0.1), (10.0, 0.2), (10.0, 0.3)]; name=(:test,), skipmissing=true) == [ - (geometry = (9.0, 0.1), test = 1) - (geometry = (10.0, 0.2), test = 4) - ] - end -end - @testset "trim, crop, extend" begin A = [missing missing missing missing 2.0 0.5 @@ -676,51 +494,6 @@ end end -@testset "mosaic" begin - reg1 = Raster([0.1 0.2; 0.3 0.4], (X(2.0:1.0:3.0), Y(5.0:1.0:6.0))) - reg2 = Raster([1.1 1.2; 1.3 1.4], (X(3.0:1.0:4.0), Y(6.0:1.0:7.0))) - irreg1 = Raster([0.1 0.2; 0.3 0.4], (X([2.0, 3.0]), Y([5.0, 6.0]))) - irreg2 = Raster([1.1 1.2; 1.3 1.4], (X([3.0, 4.0]), Y([6.0, 7.0]))) - - span_x1 = Explicit(vcat((1.5:1.0:2.5)', (2.5:1.0:3.5)')) - span_x2 = Explicit(vcat((2.5:1.0:3.5)', (3.5:1.0:4.5)')) - exp1 = Raster([0.1 0.2; 0.3 0.4], (X(Sampled([2.0, 3.0]; span=span_x1)), Y([5.0, 6.0]))) - exp2 = Raster([1.1 1.2; 1.3 1.4], (X(Sampled([3.0, 4.0]; span=span_x2)), Y([6.0, 7.0]))) - @test val(span(mosaic(first, exp1, exp2), X)) == [1.5 2.5 3.5; 2.5 3.5 4.5] - @test all(mosaic(first, [reg1, reg2]) .=== - mosaic(first, irreg1, irreg2) .=== - mosaic(first, (irreg1, irreg2)) .=== - [0.1 0.2 missing; - 0.3 0.4 1.2; - missing 1.3 1.4]) - @test all(mosaic(last, reg1, reg2) .=== - mosaic(last, irreg1, irreg2) .=== - mosaic(last, exp1, exp2) .=== [0.1 0.2 missing; - 0.3 1.1 1.2; - missing 1.3 1.4]) - - @test all(mosaic(first, [reverse(reg2; dims=Y), reverse(reg1; dims=Y)]) .=== - [missing 0.2 0.1; - 1.2 1.1 0.3; - 1.4 1.3 missing] - ) - - # 3 dimensions - A1 = Raster(ones(2, 2, 2), (X(2.0:-1.0:1.0), Y(5.0:1.0:6.0), Ti(DateTime(2001):Year(1):DateTime(2002)))) - A2 = Raster(zeros(2, 2, 2), (X(3.0:-1.0:2.0), Y(4.0:1.0:5.0), Ti(DateTime(2002):Year(1):DateTime(2003)))) - @test all(mosaic(mean, A1, A2) |> parent .=== - mosaic(mean, RasterStack(A1), RasterStack(A2)).layer1 .=== - cat([missing missing missing - missing 1.0 1.0 - missing 1.0 1.0 ], - [0.0 0.0 missing - 0.0 0.5 1.0 - missing 1.0 1.0 ], - [0.0 0.0 missing - 0.0 0.0 missing - missing missing missing], dims=3)) -end - using StableRNGs, StatsBase test = rebuild(ga; name = :test) @testset "sample" begin @@ -754,9 +527,9 @@ test = rebuild(ga; name = :test) 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)) + # extract and sample return the same type + @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) .=== diff --git a/test/mosaic.jl b/test/mosaic.jl new file mode 100644 index 000000000..b9242a976 --- /dev/null +++ b/test/mosaic.jl @@ -0,0 +1,67 @@ +using Rasters, Dates, Statistics, Test +using Rasters.Lookups, Rasters.Dimensions + +@testset "mosaic" begin + reg1 = Raster([0.1 0.2; 0.3 0.4], (X(2.0:1.0:3.0), Y(5.0:1.0:6.0))) + reg2 = Raster([1.1 1.2; 1.3 1.4], (X(3.0:1.0:4.0), Y(6.0:1.0:7.0))) + irreg1 = Raster([0.1 0.2; 0.3 0.4], (X([2.0, 3.0]), Y([5.0, 6.0]))) + irreg2 = Raster([1.1 1.2; 1.3 1.4], (X([3.0, 4.0]), Y([6.0, 7.0]))) + + span_x1 = Explicit(vcat((1.5:1.0:2.5)', (2.5:1.0:3.5)')) + span_x2 = Explicit(vcat((2.5:1.0:3.5)', (3.5:1.0:4.5)')) + exp1 = Raster([0.1 0.2; 0.3 0.4], (X(Sampled([2.0, 3.0]; span=span_x1)), Y([5.0, 6.0]))) + exp2 = Raster([1.1 1.2; 1.3 1.4], (X(Sampled([3.0, 4.0]; span=span_x2)), Y([6.0, 7.0]))) + @test val(span(mosaic(first, exp1, exp2), X)) == [1.5 2.5 3.5; 2.5 3.5 4.5] + @test all(mosaic(first, [reg1, reg2]) .=== + mosaic(first, irreg1, irreg2) .=== + mosaic(first, (irreg1, irreg2)) .=== + [0.1 0.2 missing; + 0.3 0.4 1.2; + missing 1.3 1.4]) + @test all(mosaic(last, reg1, reg2) .=== + mosaic(last, irreg1, irreg2) .=== + mosaic(last, exp1, exp2) .=== [0.1 0.2 missing; + 0.3 1.1 1.2; + missing 1.3 1.4]) + + @test all(mosaic(first, [reverse(reg2; dims=Y), reverse(reg1; dims=Y)]) .=== + [missing 0.2 0.1; + 1.2 1.1 0.3; + 1.4 1.3 missing] + ) + + @testset "Generic functions" begin + @test all(mosaic(xs -> count(x -> x > 0, xs), reg1, reg2) .=== + [1 1 missing + 1 2 1 + missing 1 1] + ) + end + + @testset "3 dimensions" begin + A1 = Raster(ones(2, 2, 2), (X(2.0:-1.0:1.0), Y(5.0:1.0:6.0), Ti(DateTime(2001):Year(1):DateTime(2002)))) + A2 = Raster(zeros(2, 2, 2), (X(3.0:-1.0:2.0), Y(4.0:1.0:5.0), Ti(DateTime(2002):Year(1):DateTime(2003)))) + mean_mos = cat([missing missing missing + missing 1.0 1.0 + missing 1.0 1.0 ], + [0.0 0.0 missing + 0.0 0.5 1.0 + missing 1.0 1.0 ], + [0.0 0.0 missing + 0.0 0.0 missing + missing missing missing], dims=3) + @test all(mosaic(mean, A1, A2) .=== + mosaic(mean, RasterStack(A1), RasterStack(A2)).layer1 .=== + mean_mos) + @test mosaic(length, A1, A2; missingval=0) == cat([0 0 0 + 0 1 1 + 0 1 1], + [1 1 0 + 1 2 1 + 0 1 1], + [1 1 0 + 1 1 0 + 0 0 0], dims=3) + end + +end diff --git a/test/runtests.jl b/test/runtests.jl index 052ba3137..a9ed7ea04 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -10,16 +10,18 @@ using Rasters, Test, Aqua, SafeTestsets 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 "adapt" begin include("adapt.jl") end +@time @safetestset "methods" begin include("methods.jl") end @time @safetestset "aggregate" begin include("aggregate.jl") end +@time @safetestset "mosaic" begin include("mosaic.jl") end @time @safetestset "rasterize" begin include("rasterize.jl") end -@time @safetestset "adapt" begin include("adapt.jl") end +@time @safetestset "extract" begin include("extract.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 @@ -35,7 +37,7 @@ if !Sys.iswindows() @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 +# Only test RasterDataSources locally for now, because CI downloads keep breaking if !haskey(ENV, "CI") @time @safetestset "rasterdatasources" begin include("sources/rasterdatasources.jl") end end diff --git a/test/sources/gdal.jl b/test/sources/gdal.jl index fc058d85b..aa4737723 100644 --- a/test/sources/gdal.jl +++ b/test/sources/gdal.jl @@ -187,7 +187,7 @@ gdalpath = maybedownload(url) @test gdalarray[Y(4.224e6..4.226e6), Band(1)] isa Raster end - @testset "methods" begin + @testset "methods" begin @testset "mean" begin @test all(mean(gdalarray; dims=Y) .=== mean(parent(gdalarray); dims=2)) end @@ -289,20 +289,26 @@ gdalpath = maybedownload(url) @testset "mosaic" begin @time gdalarray = Raster(gdalpath; name=:test) A1 = gdalarray[X(1:300), Y(1:200)] - A2 = gdalarray[X(57:500), Y(101:301)] + A2 = gdalarray[X(57:End()), Y(101:End())] tempfile1 = tempname() * ".tif" tempfile2 = tempname() * ".tif" tempfile3 = tempname() * ".tif" - Afile = mosaic(first, A1, A2; missingval=0xff, atol=1e-8, filename=tempfile1) + view(gdalarray, extent(A2)) + Afile = mosaic(first, A1, A2; missingval=0xff, atol=1e-1, filename=tempfile1) + @test missingval(Afile) === 0xff Afile2 = mosaic(first, A1, A2; atol=1e-8, filename=tempfile2) - collect(Afile) - collect(Afile2) @test missingval(Afile2) === missing - Amem = mosaic(first, A1, A2; missingval=0xff, atol=1e-8) - Atest = gdalarray[X(1:500), Y(1:301)] - Atest[X(1:56), Y(201:301)] .= 0xff - Atest[X(301:500), Y(1:100)] .= 0xff + Amem = mosaic(first, A1, A2; missingval=0xff, atol=1e-5) + + Atest = rebuild(gdalarray[Extents.union(extent(A1), extent(A2))]; missingval=0xff) + Atest .= 0xff + + Atest[DimSelectors(Atest[extent(A1)]; selectors=Contains())] .= A1 + Atest[DimSelectors(Atest[extent(A2)]; selectors=Contains())] .= A2 + + @test size(Amem) == size(gdalarray) @test all(Atest .=== Amem .=== Afile .=== replace_missing(Afile2, 0xff)) + filter(x -> !Bool(first(x)), tuple.((Atest .=== Amem), Atest, Amem)) end end # methods diff --git a/test/sources/grd.jl b/test/sources/grd.jl index ca3c72425..dc6a87b79 100644 --- a/test/sources/grd.jl +++ b/test/sources/grd.jl @@ -171,7 +171,6 @@ grdpath = stem * ".gri" A = Raster(tempgrd) @test count(==(100.0f0), A) + count(==(255.0f0), A) == length(A) end - @testset "mosaic" begin @time grdarray = Raster(grdpath) A1 = grdarray[X(1:40), Y(1:30)] diff --git a/test/sources/ncdatasets.jl b/test/sources/ncdatasets.jl index 4c49ad607..bf9337429 100644 --- a/test/sources/ncdatasets.jl +++ b/test/sources/ncdatasets.jl @@ -210,12 +210,13 @@ end A1 = ncarray[X(1:80), Y(1:100)] A2 = ncarray[X(50:150), Y(90:150)] tempfile = tempname() * ".nc" - Afile = mosaic(first, read(A1), read(A2); atol=1e-7, filename=tempfile, force=true) + Afile = mosaic(first, A1, A2; + atol=1e-7, filename=tempfile, force=true + ) |> read Amem = mosaic(first, A1, A2; atol=1e-7) Atest = ncarray[X(1:150), Y(1:150)] Atest[X(1:49), Y(101:150)] .= missing Atest[X(81:150), Y(1:89)] .= missing - read(Afile) @test all(Atest .=== Afile .=== Amem) end @testset "slice" begin @@ -400,7 +401,7 @@ end B = rebuild(boolmask(ncarray) .* 1; missingval=nothing) write(filename, B) nomissing = Raster(filename) - @test missingval(nomissing) == nothing + @test missingval(nomissing) === nothing @test eltype(nomissing) == Int64 end