From 88af714ec6194c321c6d1188e2c38c5e318318cf Mon Sep 17 00:00:00 2001 From: Rafael Schouten Date: Sun, 22 Oct 2023 00:11:06 +0200 Subject: [PATCH] use BitArray for burning buffer (#533) * optimise rasterize * threading is `false` by default * simplify locking * update threaded docstring * remove lock arg * bugfix rasterize and test threading --- ext/RastersArchGDALExt/warp.jl | 4 +- src/Rasters.jl | 1 - src/methods/burning/array_init.jl | 14 +- src/methods/burning/geometry.jl | 2 +- src/methods/burning/point.jl | 2 +- src/methods/coverage.jl | 31 ++-- src/methods/mask.jl | 14 +- src/methods/rasterize.jl | 20 +-- src/methods/shared_docstrings.jl | 9 ++ src/sectorlock.jl | 83 ---------- test/methods.jl | 6 +- test/rasterize.jl | 241 ++++++++++++++++-------------- 12 files changed, 184 insertions(+), 243 deletions(-) delete mode 100644 src/sectorlock.jl diff --git a/ext/RastersArchGDALExt/warp.jl b/ext/RastersArchGDALExt/warp.jl index ef38806f4..1559cf474 100644 --- a/ext/RastersArchGDALExt/warp.jl +++ b/ext/RastersArchGDALExt/warp.jl @@ -72,7 +72,7 @@ end function _warp(A::AbstractRaster, flags::Dict; filename=nothing, suffix="", kw...) filename = RA._maybe_add_suffix(filename, suffix) - flagvect = reduce([flags...]; init=[]) do acc, (key, val) + flagvect = reduce([flags...]; init=String[]) do acc, (key, val) append!(acc, String[_asflag(key), _stringvect(val)...]) end tempfile = isnothing(filename) ? nothing : tempname() * ".tif" @@ -84,7 +84,7 @@ function _warp(A::AbstractRaster, flags::Dict; filename=nothing, suffix="", kw.. # Either read the MEM dataset, or get the filename as a FileArray # And permute the dimensions back to what they were in A p_raster = _maybe_permute_from_gdal(raster, dims(A)) - # Either read the MEM dataset to an Array, or keep a filename base raster lazy + # Either read the MEM dataset to an Array, or keep the raster lazy return isnothing(filename) ? read(p_raster) : p_raster end end diff --git a/src/Rasters.jl b/src/Rasters.jl index 2cb7c45a7..64bc8131c 100644 --- a/src/Rasters.jl +++ b/src/Rasters.jl @@ -109,7 +109,6 @@ include("read.jl") include("write.jl") include("show.jl") include("plotrecipes.jl") -include("sectorlock.jl") include("methods/burning/edges.jl") include("methods/burning/allocs.jl") diff --git a/src/methods/burning/array_init.jl b/src/methods/burning/array_init.jl index 2be68f4d1..16ee70e00 100644 --- a/src/methods/burning/array_init.jl +++ b/src/methods/burning/array_init.jl @@ -2,7 +2,7 @@ # Like `create` but without disk writes, mostly for Bool/Union{Missing,Boo}, # and uses `similar` where possible # TODO merge this with `create` somehow -_init_bools(to; kw...) = _init_bools(to, Bool; kw...) +_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...) @@ -31,15 +31,11 @@ function _init_bools(to, dims::DimTuple, T::Type, data; collapse::Union{Bool,Not end end -function _alloc_bools(to, dims::DimTuple, ::Type{Bool}; missingval=false, metadata=NoMetadata(), kw...) - if length(dims) > 2 - # Use a BitArray - return Raster(falses(size(dims)), dims; missingval, metadata) # Use a BitArray - else - return Raster(zeros(Bool, size(dims)), dims; missingval, metadata) # Use a BitArray - end +function _alloc_bools(to, dims::DimTuple, ::Type{BitArray}; missingval=false, metadata=NoMetadata(), kw...) + # Use a BitArray + return Raster(falses(size(dims)), dims; missingval, metadata) end -function _alloc_bools(to, dims::DimTuple, ::Type{T}; missingval=false, metadata=NoMetadata(), kw...) where T +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) return rebuild(data; missingval, metadata) diff --git a/src/methods/burning/geometry.jl b/src/methods/burning/geometry.jl index 5afae6204..4263ead93 100644 --- a/src/methods/burning/geometry.jl +++ b/src/methods/burning/geometry.jl @@ -23,7 +23,7 @@ function _burn_geometry!(B::AbstractRaster, ::GI.AbstractFeatureCollectionTrait, end # Where geoms is an iterator function _burn_geometry!(B::AbstractRaster, trait::Nothing, geoms; - collapse::Union{Bool,Nothing}=nothing, lock=SectorLocks(), verbose=true, progress=true, threaded=true, + collapse::Union{Bool,Nothing}=nothing, lock=Threads.SpinLock(), verbose=true, progress=true, threaded=true, allocs=_burning_allocs(B; threaded), kw... )::Bool range = _geomindices(geoms) diff --git a/src/methods/burning/point.jl b/src/methods/burning/point.jl index bc7b0e499..935481cd5 100644 --- a/src/methods/burning/point.jl +++ b/src/methods/burning/point.jl @@ -34,7 +34,7 @@ end _fill_index!(x, fill, I) else sector = CartesianIndices(map(i -> i:i, I)) - Base.lock(lock, sector) + Base.lock(lock) _fill_index!(x, fill, I) Base.unlock(lock) end diff --git a/src/methods/coverage.jl b/src/methods/coverage.jl index 4895f9cd3..39f903725 100644 --- a/src/methods/coverage.jl +++ b/src/methods/coverage.jl @@ -24,6 +24,7 @@ const COVERAGE_KEYWORDS = """ 10 x 10 or 100 points that contribute to coverage. Using `100` means 10,000 points contribute. Performance will decline as `scale` increases. Memory use will grow by `scale^2` when `mode=:union`. +$THREADED_KEYWORD $PROGRESS_KEYWORD $VERBOSE_KEYWORD """ @@ -92,11 +93,12 @@ function _coverage!(A::AbstractRaster, ::GI.AbstractGeometryTrait, geom, r; scal return A end function _coverage!(A::AbstractRaster, ::Nothing, geoms, r::Rasterizer; mode, scale) - n = _nthreads() + n = r.threaded ? _nthreads() : 1 buffers = ( - allocs = _burning_allocs(A; threaded=r.threaded), - linebuffer = [_init_bools(A, Bool; missingval=false) for _ in 1:n], - centerbuffer = [_init_bools(A, Bool; missingval=false) for _ in 1:n], + # We always need a vector of threads, but just + allocs = _burning_allocs(A; nthreads=n, threaded=true), + linebuffer = [_init_bools(A, BitArray; missingval=false) for _ in 1:n], + centerbuffer = [_init_bools(A, BitArray; missingval=false) for _ in 1:n], block_crossings = [[Vector{Float64}(undef, 0) for _ in 1:scale] for _ in 1:n], burnstatus=[fill(BurnStatus(), scale) for _ in 1:n], subbuffer = [fill!(Array{Bool}(undef, scale, scale), false) for _ in 1:n], @@ -116,11 +118,11 @@ end # Combines coverage at the sub-pixel level for a final value 0-1 function _union_coverage!(A::AbstractRaster, geoms, buffers; - scale, subpixel_dims, progress=true, threaded=true + scale, subpixel_dims, progress=true, threaded=false ) n = _nthreads() - centeracc = [_init_bools(A, Bool; missingval=false) for _ in 1:n] - lineacc = [_init_bools(A, Bool; missingval=false) for _ in 1:n] + centeracc = [_init_bools(A, BitArray; missingval=false) for _ in 1:n] + lineacc = [_init_bools(A, BitArray; missingval=false) for _ in 1:n] subpixel_buffer = [falses(size(A) .* scale) for _ in 1:n] allbuffers = merge(buffers, (; centeracc, lineacc, subpixel_buffer)) @@ -129,6 +131,7 @@ function _union_coverage!(A::AbstractRaster, geoms, buffers; _run(range, threaded, progress, "Calculating coverage buffers...") do i geom = _getgeom(geoms, i) idx = Threads.threadid() + # Get buffers for each thread as a NamedTuple thread_buffers = map(b -> b[idx], allbuffers) _union_coverage!(A, geom; scale, subpixel_buffer, thread_buffers...) fill!(thread_buffers.linebuffer, false) @@ -170,12 +173,12 @@ function _union_coverage!(A::AbstractRaster, geoms, buffers; end function _union_coverage!(A::AbstractRaster, geom; scale, - linebuffer=_init_bools(A, Bool; missingval=false), - centerbuffer=_init_bools(A, Bool; missingval=false), + linebuffer=_init_bools(A, BitArray; missingval=false), + centerbuffer=_init_bools(A, BitArray; missingval=false), allocs=Allocs(linebuffer), subpixel_buffer=falses(size(A) .* scale), - centeracc=_init_bools(A, Bool; missingval=false), - lineacc=_init_bools(A, Bool; missingval=false), + centeracc=_init_bools(A, BitArray; missingval=false), + lineacc=_init_bools(A, BitArray; missingval=false), burnstatus=[BurnStatus() for _ in 1:scale], block_crossings=[Vector{Float64}(undef, 0) for _ in 1:scale], subbuffer=falses(scale, scale), @@ -251,7 +254,7 @@ end # Sums all coverage function _sum_coverage!(A::AbstractRaster, geoms, buffers; - scale, subpixel_dims, verbose=true, progress=true, threaded=true + scale, subpixel_dims, verbose=true, progress=true, threaded=false ) n = _nthreads() coveragebuffers = [fill!(similar(A), 0.0) for _ in 1:n] @@ -276,8 +279,8 @@ function _sum_coverage!(A::AbstractRaster, geoms, buffers; end function _sum_coverage!(A::AbstractRaster, geom; scale, - linebuffer=_init_bools(A, Bool; missingval=false), - centerbuffer=_init_bools(A, Bool; missingval=false), + linebuffer=_init_bools(A, BitArray; missingval=false), + centerbuffer=_init_bools(A, BitArray; missingval=false), allocs=Allocs(linebuffer), block_crossings=[Vector{Float64}(undef, 0) for _ in 1:scale], subbuffer=falses(scale, scale), diff --git a/src/methods/mask.jl b/src/methods/mask.jl index a08d2c5c6..25adf8787 100644 --- a/src/methods/mask.jl +++ b/src/methods/mask.jl @@ -208,9 +208,11 @@ The array returned from calling `boolmask` on a `AbstractRaster` is a - `missingval`: The missing value of the source array, with default `missingval(raster)`. -# Geometry keywords +# Keywords $GEOM_KEYWORDS +$THREADED_KEYWORD +$PROGRESS_KEYWORD And specifically for `shape=:polygon`: @@ -245,14 +247,14 @@ function boolmask end boolmask(series::AbstractRasterSeries; kw...) = boolmask(first(series); kw...) boolmask(stack::AbstractRasterStack; kw...) = boolmask(first(stack); kw...) function boolmask(source::AbstractRaster; kw...) - dest = _init_bools(source, Bool, nothing; kw..., missingval=false) + dest = _init_bools(source, BitArray, nothing; kw..., missingval=false) return boolmask!(dest, source; kw...) end function boolmask(x; to=nothing, kw...) if to isa Union{AbstractDimArray,AbstractDimStack,DimTuple} to = dims(to, DEFAULT_POINT_ORDER) end - A = _init_bools(to, Bool, x; kw..., missingval=false) + A = _init_bools(to, BitArray, x; kw..., missingval=false) return boolmask!(A, x; kw...) end @@ -264,7 +266,7 @@ function boolmask!(dest::AbstractRaster, src::AbstractRaster; end end function boolmask!(dest::AbstractRaster, geoms; - allocs=nothing, lock=nothing, progress=true, threaded=true, kw... + allocs=nothing, lock=nothing, progress=true, threaded=false, kw... ) if isnothing(allocs) allocs = _burning_allocs(dest; threaded) @@ -322,11 +324,11 @@ $EXPERIMENTAL missingmask(series::AbstractRasterSeries; kw...) = missingmask(first(series); kw...) missingmask(stack::AbstractRasterStack; kw...) = missingmask(first(stack); kw...) function missingmask(source::AbstractRaster; kw...) - dest = _init_bools(source, Union{Missing,Bool}, nothing; kw..., missingval=missing) + dest = _init_bools(source, Array{Union{Missing,Bool}}, nothing; kw..., missingval=missing) return missingmask!(dest, source; kw...) end function missingmask(x; to=nothing, kw...) - B = _init_bools(to, Union{Missing,Bool}, x; kw..., missingval=missing) + B = _init_bools(to, Array{Union{Missing,Bool}}, x; kw..., missingval=missing) return missingmask!(B, x; kw...) end diff --git a/src/methods/rasterize.jl b/src/methods/rasterize.jl index 7401be41c..207ae319c 100644 --- a/src/methods/rasterize.jl +++ b/src/methods/rasterize.jl @@ -75,7 +75,7 @@ struct Rasterizer{T,G,F,R,O,I,M} op::O init::I missingval::M - lock::Union{SectorLocks,Nothing} + lock::Union{Threads.SpinLock,Nothing} shape::Symbol boundary::Symbol verbose::Bool @@ -93,7 +93,7 @@ function Rasterizer(geom, fill, fillitr; filename=nothing, verbose=true, progress=true, - threaded=true, + threaded=false, kw... ) # A single geometry does not need a reducing function @@ -138,7 +138,7 @@ function Rasterizer(geom, fill, fillitr; @warn "currently `:points` rasterization of multiple non-`PointTrait` geometries may be innaccurate for `reducer` methods besides $stable_reductions. Make a Rasters.jl github issue if you need this to work" end eltype, missingval = get_eltype_missingval(eltype, missingval, fillitr, init, filename, op, reducer) - lock = threaded ? SectorLocks() : nothing + lock = threaded ? Threads.SpinLock() : nothing return Rasterizer(eltype, geom, fillitr, reducer, op, init, missingval, lock, shape, boundary, verbose, progress, threaded) end @@ -323,7 +323,7 @@ const RASTERIZE_KEYWORDS = """ - `progress`: show a progress bar, `true` by default, `false` to hide.. - `verbose`: print information and warnings whne there are problems with the rasterisation. `true` by default. -- `threaded`: run operations in parallel. `true` by default. +$THREADED_KEYWORD """ const RASTERIZE_ARGUMENTS = """ @@ -569,12 +569,12 @@ function _rasterize!(A, ::GI.AbstractGeometryTrait, geom, fill, r::Rasterizer; a V = view(A, Touches(ext)) length(V) > 0 || return false - bools = _init_bools(commondims(V, DEFAULT_POINT_ORDER), Bool; metadata=metadata(A)) + bools = _init_bools(commondims(V, DEFAULT_POINT_ORDER), BitArray; metadata=metadata(A)) boolmask!(bools, geom; allocs, lock, shape, boundary, verbose, progress) hasburned = any(bools) if hasburned - # Avoid race conditions with a SectorLock - isnothing(lock) || Base.lock(lock, V) + # Avoid race conditions + isnothing(lock) || Base.lock(lock) _fill!(V, bools, fill, op, init, missingval) isnothing(lock) || Base.unlock(lock) end @@ -584,11 +584,7 @@ end # Fill points function _rasterize!(A, trait::GI.AbstractPointTrait, point, fill, r::Rasterizer; allocs=nothing) # Avoid race conditions whern Point is in a mixed set of Geometries - # isnothing(r.lock) || Base.lock(r.lock, A) - hasburned = _fill_point!(A, trait, point; fill, r.lock) - # isnothing(r.lock) || Base.unlock(r.lock) - # for all points we avoid parallel rasterization completely - this method should not be hit often - return hasburned + return _fill_point!(A, trait, point; fill, r.lock) end function _rasterize!(A, trait::Nothing, geoms, fill, r::Rasterizer; allocs=nothing) if r.shape === :point diff --git a/src/methods/shared_docstrings.jl b/src/methods/shared_docstrings.jl index 152f73aa1..f184534c4 100644 --- a/src/methods/shared_docstrings.jl +++ b/src/methods/shared_docstrings.jl @@ -27,6 +27,15 @@ const SHAPE_KEYWORDS = """ The default is `:center`. """ +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 + geometries written into low-resolution rasters, but may not be for simple geometries with + high-resolution rasters. With very large rasters threading may be counter productive due + to excessing memory use. Caution should also be used: `threaded` should not be used in in-place + functions wrinting to `BitArray` or other arrays where race conditions can occur. +""" + const GEOM_KEYWORDS = """ $TO_KEYWORD $RES_KEYWORD diff --git a/src/sectorlock.jl b/src/sectorlock.jl deleted file mode 100644 index 28cc3637f..000000000 --- a/src/sectorlock.jl +++ /dev/null @@ -1,83 +0,0 @@ -############################## -# SectorLocks - -# These lets us lock part of the array so we can work in parallel -# without always blocking other threads. - -# Usually most polygons cover a small subset of a raster. -mutable struct SectorLock - lock::Threads.SpinLock - sector::CartesianIndices{2,Tuple{UnitRange{Int},UnitRange{Int}}} -end -SectorLock() = SectorLock(Threads.SpinLock(), CartesianIndices((1:0, 1:0))) - -Base.lock(sl::SectorLock) = lock(sl.lock) -Base.unlock(sl::SectorLock) = unlock(sl.lock) -Base.islocked(sl::SectorLock) = islocked(sl.lock) - -# One lock for each thread - this is only for threads -struct SectorLocks - seclocks::Vector{SectorLock} - spinlock::Threads.SpinLock -end -SectorLocks(n::Int) = SectorLocks([SectorLock() for _ in 1:n], Threads.SpinLock()) -SectorLocks() = SectorLocks(_nthreads()) - -Base.length(sl::SectorLocks) = length(sl.seclocks) -Base.getindex(sl::SectorLocks, i::Int) = sl.seclocks[i] -Base.iterate(sl::SectorLocks, args...) = iterate(sl.seclocks, args...) -Base.eachindex(sl::SectorLocks) = 1:length(sl) - -Base.lock(sl::SectorLocks, sector::RasterStack) = Base.lock(sl, first(sector)) -function Base.lock(sl::SectorLocks, A::Raster) - o = otherdims(A, DEFAULT_POINT_ORDER) - if length(o) > 0 - slice = first(eachslice(A; dims=o)) - return Base.lock(sl, parent(slice)) - else - return Base.lock(sl, parent(A)) - end -end -# Base.lock(sl::SectorLocks, sector::SubArray) = Base.lock(sl, CartesianIndices(_unitranges(sector.indices...))) -# Base.lock(sl::SectorLocks, sector::DiskArrays.SubDiskArray) = Base.lock(sl, sector.v) -Base.lock(sl::SectorLocks, sector::Tuple) = Base.lock(sl, CartesianIndices(sector)) -Base.lock(sl::SectorLocks, sector::AbstractArray) = Base.lock(sl, CartesianIndices(sector)) -function Base.lock(seclocks::SectorLocks, sector::CartesianIndices) - idx = Threads.threadid() - thread_lock = seclocks.seclocks[idx] - thread_lock.sector = sector - # Lock the main lock so no other sectors can be locked - lock(seclocks.spinlock) - # Check for any other lock that intersects our sector - for i in eachindex(seclocks) - # Ignore a lock from this thread - i == idx && continue - seclock = seclocks[i] - # Most of the time the lock will not intersect - if islocked(seclock) && _intersects(seclock, sector) - # If it does, lock the sector and wait - # All other threads are also waiting, - # but if this doens't happen often its fine. - lock(seclock) - unlock(seclock) - end - end - # Lock this sector - lock(thread_lock) - # Unlock the main spinlock so other sectors can be locked - unlock(seclocks.spinlock) - # Work is now done in the locked sector, knowing no other thread - # can write to it. `unlock` is run when it is finished. -end -Base.unlock(sl::SectorLocks) = begin - # println("unlocking thread", idx) - unlock(sl[Threads.threadid()]) -end - -_intersects(seclock::SectorLock, sector) = _intersects(seclock.sector, sector) -_intersects(sector1, sector2) = length(intersect(sector1, sector2)) > 0 - -_unitranges(x1, xs...) = (_unitranges(x1)..., _unitranges(xs...)...) -_unitranges(x1::UnitRange) = (x1,) -_unitranges(x1::AbstractRange{Int}) = (UnitRange(min(first(x1), last(x1)), max(first(x1), last(x1))),) -_unitranges(x1) = () diff --git a/test/methods.jl b/test/methods.jl index 2db44ad8e..aef73bde2 100644 --- a/test/methods.jl +++ b/test/methods.jl @@ -65,13 +65,13 @@ end @testset "boolmask" begin @test boolmask(ga) == [false true; true false] - @test parent(boolmask(ga)) isa Matrix{Bool} + @test parent(boolmask(ga)) isa BitMatrix @test boolmask(ga99) == [false true; true false] @test boolmask(gaNaN) == [false true; true false] @test dims(boolmask(ga)) == (X(NoLookup(Base.OneTo(2))), Y(NoLookup(Base.OneTo(2)))) x = boolmask(polygon; res=1.0) @test x == trues(X(Projected(-20:1.0:-1.0; crs=nothing)), Y(Projected(10.0:1.0:29.0; crs=nothing))) - @test parent(x) isa Matrix{Bool} + @test parent(x) isa BitMatrix # With a :geometry axis x = boolmask([polygon, polygon]; collapse=false, res=1.0) @test eltype(x) == Bool @@ -81,7 +81,7 @@ end x = boolmask([polygon, polygon]; collapse=true, res=1.0) @test size(x) == (20, 20) @test sum(x) == 400 - @test parent(x) isa Matrix{Bool} + @test parent(x) isa BitMatrix end @testset "missingmask" begin diff --git a/test/rasterize.jl b/test/rasterize.jl index 8c166524e..c11c06cc9 100644 --- a/test/rasterize.jl +++ b/test/rasterize.jl @@ -49,57 +49,66 @@ st = RasterStack((A1, copy(A1))) geom = line_collection geom = poly_collection - for A in (A1, A2), geom in (pointvec, pointfc, multi_point, linestring, multi_linestring, linearring, polygon, multi_polygon, table, line_collection, poly_collection) + for A in (A1, A2), + geom in (pointvec, pointfc, multi_point, linestring, multi_linestring, linearring, polygon, multi_polygon, table, line_collection, poly_collection), + threaded in (true, false) + fill!(A, 0) - rasterize!(sum, A, geom; shape=:point, fill=1); + rasterize!(sum, A, geom; shape=:point, fill=1, threaded); @test sum(A) == 5.0 - @test sum(rasterize(sum, geom; to=A, shape=:point, fill=1, missingval=0)) == 5.0 - rasterize!(last, A, geom; shape=:point, fill=1); + @test sum(rasterize(sum, geom; to=A, shape=:point, fill=1, missingval=0, threaded)) == 5.0 + rasterize!(last, A, geom; shape=:point, fill=1, threaded); @test sum(A) == 4.0 - @test sum(rasterize(last, geom; to=A, shape=:point, fill=1, missingval=0)) == 4.0 + @test sum(rasterize(last, geom; to=A, shape=:point, fill=1, missingval=0, threaded)) == 4.0 fill!(A, 0) if !Tables.istable(geom) - rasterize!(count, A, [geom, geom]; shape=:point) + rasterize!(count, A, [geom, geom]; shape=:point, threaded) @test sum(A) == 10.0 fill!(A, 0) end end geom = multi_point - for A in (A1, A2), geom in (table, pointvec, pointfc, multi_point, linestring, multi_linestring, linearring, polygon, multi_polygon, line_collection, poly_collection) - rasterize!(sum, st, geom; shape=:point, fill=(layer1=2, layer2=3)) + for A in (A1, A2), + geom in (table, pointvec, pointfc, multi_point, linestring, multi_linestring, linearring, polygon, multi_polygon, line_collection, poly_collection), + threaded in (true, false) + + rasterize!(sum, st, geom; shape=:point, fill=(layer1=2, layer2=3), threaded) st.layer1 .= st.layer2 .= 0 - rasterize!(sum, st, geom; shape=:point, fill=(layer1=2, layer2=3)) + rasterize!(sum, st, geom; shape=:point, fill=(layer1=2, layer2=3), threaded) @test sum(st.layer1) == 10 @test sum(st.layer2) == 15 @test parent(st.layer1) isa Array{Float64,2} st.layer1 .= 0; st.layer2 .= 0 - rasterize!(last, st, geom; shape=:point, fill=(layer1=2, layer2=3)) + rasterize!(last, st, geom; shape=:point, fill=(layer1=2, layer2=3), threaded) @test sum(st.layer1) == 8 @test sum(st.layer2) == 12 st.layer1 .= 0; st.layer2 .= 0 - rasterize!(first, st, geom; shape=:point, fill=(layer1=2, layer2=3)) + rasterize!(first, st, geom; shape=:point, fill=(layer1=2, layer2=3), threaded) @test sum(st.layer1) == 8 @test sum(st.layer2) == 12 st.layer1 .= 0; st.layer2 .= 0 end - for A in (A1, A2), geom in (table, pointvec, pointfc, multi_point, linestring, multi_linestring, linearring, polygon, multi_polygon) + for A in (A1, A2), + geom in (table, pointvec, pointfc, multi_point, linestring, multi_linestring, linearring, polygon, multi_polygon), + threaded in (true, false) + st[:layer1] .= 0; st[:layer2] .= 0 - rasterize!(sum, st, geom; shape=:point, fill=(layer1=1:5, layer2=6:10)) + rasterize!(sum, st, geom; shape=:point, fill=(layer1=1:5, layer2=6:10), threaded) @test sum(st[:layer1]) == sum(1:5) @test sum(st[:layer2]) == sum(6:10) @test parent(st[:layer1]) isa Array{Float64,2} st[:layer1] .= 0; st[:layer2] .= 0 - rasterize!(last, st, geom; shape=:point, fill=(layer1=1:5, layer2=6:10)) + rasterize!(last, st, geom; shape=:point, fill=(layer1=1:5, layer2=6:10), threaded) @test sum(st[:layer1]) == sum(2:5) @test sum(st[:layer2]) == sum(7:10) st[:layer1] .= 0; st[:layer2] .= 0 - rasterize!(first, st, geom; shape=:point, fill=(layer1=1:5, layer2=6:10)) + rasterize!(first, st, geom; shape=:point, fill=(layer1=1:5, layer2=6:10), threaded) @test sum(st[:layer1]) == sum(1:4) @test sum(st[:layer2]) == sum(6:9) st[:layer1] .= 0; st[:layer2] .= 0 - @test_nowarn rasterize!(sum, A, geom; shape=:point, fill=1) + @test_nowarn rasterize!(sum, A, geom; shape=:point, fill=1, threaded) end end @@ -107,18 +116,21 @@ end A = A1 geom = linestring geom = line_collection - for A in (A1, A2), geom in (linestring, multi_linestring, linearring, polygon, multi_polygon, line_collection, poly_collection) + for A in (A1, A2), + geom in (linestring, multi_linestring, linearring, polygon, multi_polygon, line_collection, poly_collection), + threaded in (true, false) + A .= 0 - rasterize!(sum, A, geom; shape=:line, fill=1) + rasterize!(sum, A, geom; shape=:line, fill=1, threaded) @test sum(A) == 20 + 20 + 20 + 20 - @test sum(rasterize(sum, geom; to=A, shape=:line, fill=1, missingval=0)) == 80 + @test sum(rasterize(sum, geom; to=A, shape=:line, fill=1, missingval=0, threaded)) == 80 end @testset ":line is detected for line geometries" begin - for A in (A1, A2), geom in (linestring, multi_linestring) + for A in (A1, A2), geom in (linestring, multi_linestring), threaded in (true, false) A .= 0 - rasterize!(A, geom; fill=1) + rasterize!(A, geom; fill=1, threaded) @test sum(A) == 20 + 20 + 20 + 20 - @test sum(rasterize(geom; to=A, fill=1, missingval=0)) == 80 + @test sum(rasterize(geom; to=A, fill=1, missingval=0, threaded)) == 80 end end end @@ -127,38 +139,38 @@ end A = A1 poly = polygon poly = poly_collection - for A in (A1, A2), poly in (polygon, multi_polygon, poly_collection) + for A in (A1, A2), poly in (polygon, multi_polygon, poly_collection), threaded in (true, false) A .= 0 - ra = rasterize(last, poly; to=A, missingval=0, shape=:polygon, fill=1, boundary=:center) - ra_res = rasterize(last, poly; res=map(step, span(A)), missingval=0, shape=:polygon, fill=1, boundary=:center) + ra = rasterize(last, poly; to=A, missingval=0, shape=:polygon, fill=1, boundary=:center, threaded) + ra_res = rasterize(last, poly; res=map(step, span(A)), missingval=0, shape=:polygon, fill=1, boundary=:center, threaded) @test parent(ra) isa Matrix{Int} @test sum(ra) == sum(ra_res) === 20 * 20 - ra = rasterize(last, poly; to=A, shape=:polygon, fill=1, boundary=:touches) + ra = rasterize(last, poly; to=A, shape=:polygon, fill=1, boundary=:touches, threaded) @test parent(ra) isa Array{Union{Missing,Int},2} @test sum(skipmissing(ra)) === 21 * 21 - rasterize!(last, A, poly; shape=:polygon, fill=1, boundary=:inside) + rasterize!(last, A, poly; shape=:polygon, fill=1, boundary=:inside, threaded) @test sum(A) === 19.0 * 19.0 A .= 0 @testset "polygon is detected for polygon geometries" begin A = Raster(zeros(X(-20:5; sampling=Intervals()), Y(0:30; sampling=Intervals()))) - R = rasterize(last, poly; to=A, fill=1) + R = rasterize(last, poly; to=A, fill=1, threaded) @test parent(R) isa Array{} @test sum(skipmissing(R)) == 20 * 20 - @test_throws ArgumentError rasterize!(A, poly; shape=:notashape, fill=1) - @test_throws ArgumentError rasterize!(A, poly; shape=:polygon, fill=1, boundary=:notaboundary) + @test_throws ArgumentError rasterize!(A, poly; shape=:notashape, fill=1, threaded) + @test_throws ArgumentError rasterize!(A, poly; shape=:polygon, fill=1, boundary=:notaboundary, threaded) end - st1 = rasterize(last, poly; fill=(layer1=1, layer2=2), to=st) + st1 = rasterize(last, poly; fill=(layer1=1, layer2=2), to=st, threaded) @test sum(skipmissing(st1[:layer1])) == 400 # The last value overwrites the first @test sum(skipmissing(st1[:layer2])) == 800 # Missing size / res - @test_throws ArgumentError rasterize(poly; fill=1) + @test_throws ArgumentError rasterize(poly; fill=1, threaded) # Both size + res - @test_throws ArgumentError rasterize(poly; res=0.1, size=200, fill=1) - @test_throws ArgumentError rasterize(poly; res=(0.1, 0.2), size=200, fill=1) - @test_throws ArgumentError rasterize(poly; res=0.1, size=(200, 200), fill=1) - @test_throws ArgumentError rasterize(poly; res=(0.1, 0.2), size=(200, 200), fill=1) + @test_throws ArgumentError rasterize(poly; res=0.1, size=200, fill=1, threaded) + @test_throws ArgumentError rasterize(poly; res=(0.1, 0.2), size=200, fill=1, threaded) + @test_throws ArgumentError rasterize(poly; res=0.1, size=(200, 200), fill=1, threaded) + @test_throws ArgumentError rasterize(poly; res=(0.1, 0.2), size=(200, 200), fill=1, threaded) end end @@ -210,14 +222,14 @@ end @testset "feature collection, table from fill of Symbol keys" begin data = pointfc Tables.istable(data) - for data in (pointfc, DataFrame(pointfc)) + for data in (pointfc, DataFrame(pointfc)), threaded in (true, false) @testset "NTuple of Symbol fill makes an stack" begin - rst = rasterize(sum, data; to=A, fill=(:val1, :val2)) + rst = rasterize(sum, data; to=A, fill=(:val1, :val2), threaded) @test parent(rst.val1) isa Array{Union{Missing,Int},2} @test parent(rst.val2) isa Array{Union{Missing,Float32},2} @test keys(rst) == (:val1, :val2) @test map(sum ∘ skipmissing, rst) === (val1=15, val2=30.0f0) - @test_throws ArgumentError rasterize(data; to=A, fill=(:val1, :not_a_column)) + @test_throws ArgumentError rasterize(data; to=A, fill=(:val1, :not_a_column), threaded) end @testset "Symbol fill makes an array" begin ra = rasterize(sum, data; to=A, fill=:val1) @@ -352,77 +364,79 @@ end polygon = ArchGDAL.createpolygon(pointvec) polygons = ArchGDAL.createpolygon.([[pointvec1], [pointvec2], [pointvec3], [pointvec4]]) # With fill of 1 these are all the same thing - for f in (last, first, mean, median, maximum, minimum) - r = rasterize(f, polygons; res=5, fill=1, boundary=:center, threaded=false, crs=EPSG(4326)) - @test parent(r) isa Array{<:Union{Missing,<:Real},2} - @test sum(skipmissing(r)) == 12 + 12 + 12 + 16 - @test crs(r) == EPSG(4326) - end - for f in (last, maximum) - r = rasterize(last, polygons; res=5, fill=1:4, boundary=:center) - @test parent(r) isa Array{Union{Missing,Int},2} - @test sum(skipmissing(r)) == 12 * 1 + 12 * 2 + 12 * 3 + 16 * 4 - end - for f in (first, minimum) - r = rasterize(f, polygons; res=5, fill=1:4, boundary=:center) - @test parent(r) isa Array{Union{Missing,Int},2} - @test sum(skipmissing(r)) == 16 * 1 + 12 * 2 + 12 * 3 + 12 * 4 - end - for f in (mean, median) - r = rasterize(f, polygons; res=5, fill=1:4, boundary=:center) - @test parent(r) isa Array{Union{Missing,Float64},2} - @test sum(skipmissing(r)) == - (12 * 1 + 8 * 2 + 8 * 3 + 12 * 4) + (4 * 1.5 + 4 * 2.5 + 4 * 3.5) + for threaded in (true, false) + for f in (last, first, mean, median, maximum, minimum) + r = rasterize(f, polygons; res=5, fill=1, boundary=:center, threaded, crs=EPSG(4326)) + @test parent(r) isa Array{<:Union{Missing,<:Real},2} + @test sum(skipmissing(r)) == 12 + 12 + 12 + 16 + @test crs(r) == EPSG(4326) + end + for f in (last, maximum) + r = rasterize(last, polygons; res=5, fill=1:4, boundary=:center, threaded) + @test parent(r) isa Array{Union{Missing,Int},2} + @test sum(skipmissing(r)) == 12 * 1 + 12 * 2 + 12 * 3 + 16 * 4 + end + for f in (first, minimum) + r = rasterize(f, polygons; res=5, fill=1:4, boundary=:center, threaded) + @test parent(r) isa Array{Union{Missing,Int},2} + @test sum(skipmissing(r)) == 16 * 1 + 12 * 2 + 12 * 3 + 12 * 4 + end + for f in (mean, median) + r = rasterize(f, polygons; res=5, fill=1:4, boundary=:center) + @test parent(r) isa Array{Union{Missing,Float64},2} + @test sum(skipmissing(r)) == + (12 * 1 + 8 * 2 + 8 * 3 + 12 * 4) + (4 * 1.5 + 4 * 2.5 + 4 * 3.5) + end + prod_r = rasterize(prod, polygons; res=5, fill=1:4, boundary=:center, filename="test.tif", threaded) + prod_r = rasterize(prod, polygons; res=5, fill=1:4, boundary=:center, threaded) + @test sum(skipmissing(prod_r)) == + (12 * 1 + 8 * 2 + 8 * 3 + 12 * 4) + (4 * 1 * 2 + 4 * 2 * 3 + 4 * 3 * 4) + + prod_st = rasterize(prod, polygons; res=5, fill=(a=1:4, b=4:-1:1), missingval=missing, boundary=:center, threaded) + @test all(prod_st.a .=== rot180(prod_st.b)) + @test all(prod_r .=== prod_st.a) + prod_r_m = rasterize(prod, polygons; res=5, fill=1:4, missingval=-1, boundary=:center, threaded) + prod_st_m = rasterize(prod, polygons; res=5, fill=(a=1:4, b=4.0:-1.0:1.0), missingval=(a=-1, b=-1.0), boundary=:center, threaded) + @test all(prod_st_m.a .=== prod_r_m) + @test all( prod_st_m.b .=== rot180(Float64.(prod_r_m))) + + r = rasterize(last, polygons; res=5, fill=(a=1, b=2), boundary=:center, threaded) + @test all(r.a .* 2 .=== r.b) + + reduced_raster_sum_center = rasterize(sum, polygons; res=5, fill=1, boundary=:center, threaded) + reduced_raster_count_center = rasterize(count, polygons; res=5, fill=1, boundary=:center, threaded) + @test name(reduced_raster_sum_center) == :sum + @test name(reduced_raster_count_center) == :count + @test sum(skipmissing(reduced_raster_sum_center)) == + sum(skipmissing(reduced_raster_count_center)) == 16 * 4 + reduced_raster_sum_touches = rasterize(sum, polygons; res=5, fill=1, boundary=:touches, threaded) + reduced_raster_count_touches = rasterize(count, polygons; res=5, fill=1, boundary=:touches, threaded) + @test name(reduced_raster_sum_touches) == :sum + @test name(reduced_raster_count_touches) == :count + # plot(reduced_raster_count_touches) + # plot(reduced_raster_sum_touches) + # This is broken because the raster area isn't big enough + @test_broken sum(skipmissing(reduced_raster_sum_touches)) == + sum(skipmissing(reduced_raster_count_touches)) == 25 * 4 + @test sum(skipmissing(reduced_raster_sum_touches)) == + sum(skipmissing(reduced_raster_count_touches)) == 25 * 4 - 9 + # The outlines of these plots should exactly mactch, + # with three values of 2 on the diagonal + # using Plots + # Plots.plot(reduced_raster; clims=(0, 3)) + # Plots.plot!(polygons; opacity=0.3, fillcolor=:black) + reduced_center = rasterize(sum, polygons; res=5, fill=1, boundary=:center, threaded) + reduced_touches = rasterize(sum, polygons; res=5, fill=1, boundary=:touches, threaded) + reduced_inside = rasterize(sum, polygons; res=5, fill=1, boundary=:inside, threaded) + # Plots.plot(reduced_inside; clims=(0, 3)) + # Plots.plot(reduced_center; clims=(0, 3)) + # Plots.plot(reduced_touches; clims=(0, 3)) + # Plots.plot!(polygons; opacity=0.3, fillcolor=:black) + # Its not clear what the results here should be - there + # are differences between different implementations. + # Soon we will define the pixel intervals so we don't need + # arbitrary choices of which lines are touched for :touches end - prod_r = rasterize(prod, polygons; res=5, fill=1:4, boundary=:center, filename="test.tif") - prod_r = rasterize(prod, polygons; res=5, fill=1:4, boundary=:center) - @test sum(skipmissing(prod_r)) == - (12 * 1 + 8 * 2 + 8 * 3 + 12 * 4) + (4 * 1 * 2 + 4 * 2 * 3 + 4 * 3 * 4) - - prod_st = rasterize(prod, polygons; res=5, fill=(a=1:4, b=4:-1:1), missingval=missing, boundary=:center) - @test all(prod_st.a .=== rot180(prod_st.b)) - @test all(prod_r .=== prod_st.a) - prod_r_m = rasterize(prod, polygons; res=5, fill=1:4, missingval=-1, boundary=:center) - prod_st_m = rasterize(prod, polygons; res=5, fill=(a=1:4, b=4.0:-1.0:1.0), missingval=(a=-1, b=-1.0), boundary=:center) - @test all(prod_st_m.a .=== prod_r_m) - @test all( prod_st_m.b .=== rot180(Float64.(prod_r_m))) - - r = rasterize(last, polygons; res=5, fill=(a=1, b=2), boundary=:center) - @test all(r.a .* 2 .=== r.b) - - reduced_raster_sum_center = rasterize(sum, polygons; res=5, fill=1, boundary=:center) - reduced_raster_count_center = rasterize(count, polygons; res=5, fill=1, boundary=:center) - @test name(reduced_raster_sum_center) == :sum - @test name(reduced_raster_count_center) == :count - @test sum(skipmissing(reduced_raster_sum_center)) == - sum(skipmissing(reduced_raster_count_center)) == 16 * 4 - reduced_raster_sum_touches = rasterize(sum, polygons; res=5, fill=1, boundary=:touches) - reduced_raster_count_touches = rasterize(count, polygons; res=5, fill=1, boundary=:touches) - @test name(reduced_raster_sum_touches) == :sum - @test name(reduced_raster_count_touches) == :count - # plot(reduced_raster_count_touches) - # plot(reduced_raster_sum_touches) - # This is broken because the raster area isn't big enough - @test_broken sum(skipmissing(reduced_raster_sum_touches)) == - sum(skipmissing(reduced_raster_count_touches)) == 25 * 4 - @test sum(skipmissing(reduced_raster_sum_touches)) == - sum(skipmissing(reduced_raster_count_touches)) == 25 * 4 - 9 - # The outlines of these plots should exactly mactch, - # with three values of 2 on the diagonal - # using Plots - # Plots.plot(reduced_raster; clims=(0, 3)) - # Plots.plot!(polygons; opacity=0.3, fillcolor=:black) - reduced_center = rasterize(sum, polygons; res=5, fill=1, boundary=:center) - reduced_touches = rasterize(sum, polygons; res=5, fill=1, boundary=:touches) - reduced_inside = rasterize(sum, polygons; res=5, fill=1, boundary=:inside) - # Plots.plot(reduced_inside; clims=(0, 3)) - # Plots.plot(reduced_center; clims=(0, 3)) - # Plots.plot(reduced_touches; clims=(0, 3)) - # Plots.plot!(polygons; opacity=0.3, fillcolor=:black) - # Its not clear what the results here should be - there - # are differences between different implementations. - # Soon we will define the pixel intervals so we don't need - # arbitrary choices of which lines are touched for :touches end @testset "Rasterizing with extra dimensions" begin @@ -446,8 +460,13 @@ end end @testset "coverage" begin - @time covsum = coverage(sum, shphandle.shapes; res=1, scale=10) - @time covunion = coverage(union, shphandle.shapes; res=1, scale=10) + @time covsum = coverage(sum, shphandle.shapes; threaded=false, res=1, scale=10) + @time covunion = coverage(union, shphandle.shapes; threaded=false, res=1, scale=10) + @test parent(covsum) isa Array{Float64,2} + @test parent(covunion) isa Array{Float64,2} + + @time covsum = coverage(sum, shphandle.shapes; threaded=true, res=1, scale=10) + @time covunion = coverage(union, shphandle.shapes; threaded=true, res=1, scale=10) @test parent(covsum) isa Array{Float64,2} @test parent(covunion) isa Array{Float64,2} # using Plots