From 79c984e98310958199631779341eb209a28402a9 Mon Sep 17 00:00:00 2001 From: rafaqz Date: Tue, 3 Oct 2023 19:43:43 +0200 Subject: [PATCH] fix tests --- src/Rasters.jl | 2 + src/methods/burning/point.jl | 9 + src/methods/rasterize.jl | 29 +- src/polygon_ops.jl | 778 +++++++++++++++++++++++++++++++++++ src/utils.jl | 98 ----- test/rasterize.jl | 62 +-- 6 files changed, 845 insertions(+), 133 deletions(-) create mode 100644 src/polygon_ops.jl diff --git a/src/Rasters.jl b/src/Rasters.jl index 9c24ca8c0..46795200b 100644 --- a/src/Rasters.jl +++ b/src/Rasters.jl @@ -119,6 +119,8 @@ include("methods/burning/geometry.jl") include("methods/burning/point.jl") include("methods/burning/line.jl") include("methods/burning/polygon.jl") +include("methods/burning/extents.jl") +include("methods/burning/utils.jl") include("methods/shared_docstrings.jl") include("methods/mask.jl") diff --git a/src/methods/burning/point.jl b/src/methods/burning/point.jl index fecad73fa..bc7b0e499 100644 --- a/src/methods/burning/point.jl +++ b/src/methods/burning/point.jl @@ -1,7 +1,16 @@ +# TODO use rasterize_points for this instead # _fill_point! # Fill a raster with `fill` where points are inside raster pixels @noinline _fill_point!(x::RasterStackOrArray, geom; kw...) = _fill_point!(x, GI.geomtrait(geom), geom; kw...) +@noinline function _fill_point!(x::RasterStackOrArray, ::GI.GeometryCollectionTrait, geom; kw...) + _without_mapped_crs(x) do x1 + for geom in GI.getgeom(geom) + _fill_point!(x, geom; kw...) + end + end + return true +end @noinline function _fill_point!(x::RasterStackOrArray, ::GI.AbstractGeometryTrait, geom; kw...) # Just find which pixels contain the points, and set them to true _without_mapped_crs(x) do x1 diff --git a/src/methods/rasterize.jl b/src/methods/rasterize.jl index 39c468104..7401be41c 100644 --- a/src/methods/rasterize.jl +++ b/src/methods/rasterize.jl @@ -632,11 +632,29 @@ function _rasterize_points!(A, ::GI.AbstractGeometryTrait, geom, fill, r::Raster fill1 =_iterable_fill(points, fill) _rasterize_points!(A, nothing, points, fill1, r) end +function _rasterize_points!(A, ::GI.GeometryCollectionTrait, collection, fill, r::Rasterizer) + # TODO How to handle fill when there is another level of nesting + hasburned = false + for geom in _getgeom(collection) + hasburned |= _rasterize_points!(A, geom, fill, r) + end + return hasburned +end function _rasterize_points!(A, ::Nothing, geoms, fillitr, r::Rasterizer) (; reducer, op, missingval, init) = r t1 = GI.trait(first(skipmissing(geoms))) hasburned = false - if !(t1 isa GI.PointTrait) + if t1 isa GI.PointTrait + # Get extent information to properly shift the points + # to the region of the array during rounding + ext = Extents.extent(A) + xrange = ext.X[2] - ext.X[1] + yrange = ext.Y[2] - ext.Y[1] + xsize = size(A, X) + ysize = size(A, Y) + s = (; ext, xrange, yrange, xsize, ysize) + return _rasterize_points_inner!(A, geoms, fillitr, s, reducer, op, missingval, init) + else # Recurse down until we hit points if r.fillitr isa Function for geom in _getgeom(geoms) @@ -656,15 +674,6 @@ function _rasterize_points!(A, ::Nothing, geoms, fillitr, r::Rasterizer) end return hasburned end - # Get extent information to properly shift the points - # to the region of the array during rounding - ext = Extents.extent(A) - xrange = ext.X[2] - ext.X[1] - yrange = ext.Y[2] - ext.Y[1] - xsize = size(A, X) - ysize = size(A, Y) - s = (; ext, xrange, yrange, xsize, ysize) - _rasterize_points_inner!(A, geoms, fillitr, s, reducer, op, missingval, init) end @noinline function _rasterize_points_inner!(A, geoms, fillitr::F, s, reducer::R, op::O, missingval, init)::Bool where {F,O,R} diff --git a/src/polygon_ops.jl b/src/polygon_ops.jl new file mode 100644 index 000000000..56b540afe --- /dev/null +++ b/src/polygon_ops.jl @@ -0,0 +1,778 @@ +const DEFAULT_POINT_ORDER = (X(), Y()) +const DEFAULT_TABLE_DIM_KEYS = (:X, :Y) + +# Tracks the burning status for each column +struct BurnStatus + ic::Int + burn::Bool + hasburned::Bool +end +BurnStatus() = BurnStatus(1, false, false) + +# Simple point positions offset relative to the raster +struct Position + offset::Tuple{Float64,Float64} + yind::Int32 +end +function Position(point::Tuple, start::Tuple, step::Tuple) + (x, y) = point + xoff, yoff = (x - start[1]), (y - start[2]) + offset = (xoff / step[1] + 1.0, yoff / step[2] + 1.0) + yind = trunc(Int32, offset[2]) + return Position(offset, yind) +end + +# Simple edge offset with (x, y) start relative to the raster +# gradient of line and integer start/stop for y +struct Edge + start::Tuple{Float64,Float64} + gradient::Float64 + iystart::Int32 + iystop::Int32 + function Edge(start::Position, stop::Position) + if start.offset[2] > stop.offset[2] + stop, start = start, stop + end + gradient = (stop.offset[1] - start.offset[1]) / (stop.offset[2] - start.offset[2]) + new(start.offset, gradient, start.yind + 1, stop.yind) + end +end + +Base.isless(e1::Edge, e2::Edge) = isless(e1.iystart, e2.iystart) +Base.isless(e::Edge, x::Real) = isless(e.iystart, x) +Base.isless(x::Real, e::Edge) = isless(x, e.iystart) + +x_at_y(e::Edge, y) = (y - e.start[2]) * e.gradient + e.start[1] + + + +function can_skip(prevpos::Position, nextpos::Position, xlookup, ylookup) + # ignore edges between grid lines on y axis + (nextpos.yind == prevpos.yind) && + (prevpos.offset[2] != prevpos.yind) && + (nextpos.offset[2] != nextpos.yind) && return true + # ignore edges outside the grid on the y axis + (prevpos.offset[2] < 0) && (nextpos.offset[2] < 0) && return true + (prevpos.offset[2] > lastindex(ylookup)) && (nextpos.offset[2] > lastindex(ylookup)) && return true + # ignore horizontal edges + (prevpos.offset[2] == nextpos.offset[2]) && return true + return false +end + +struct Allocs{B} + buffer::B + edges::Vector{Edge} + scratch::Vector{Edge} + crossings::Vector{Float64} +end +function Allocs(buffer) + edges = Vector{Edge}(undef, 0) + scratch = Vector{Edge}(undef, 0) + crossings = Vector{Float64}(undef, 0) + return Allocs(buffer, edges, scratch, crossings) +end + +function _burning_allocs(x; nthreads=_nthreads(), threaded=true, kw...) + dims = commondims(x, DEFAULT_POINT_ORDER) + if threaded + [Allocs(_init_bools(dims; metadata=Metadata())) for _ in 1:nthreads] + else + Allocs(_init_bools(dims; metadata=Metadata())) + end +end + +_get_alloc(allocs::Vector{<:Allocs}) = + _get_alloc(allocs[Threads.threadid()]) +_get_alloc(allocs::Allocs) = allocs + + +struct Edges <: AbstractVector{Edge} + edges::Vector{Edge} + max_ylen::Int + edge_count::Int +end +Edges(geom, dims; kw...) = Edges(GI.geomtrait(geom), geom, dims; kw...) +function Edges( + tr::Union{GI.AbstractCurveTrait,GI.AbstractPolygonTrait,GI.AbstractMultiPolygonTrait}, + geom, dims; + allocs::Union{Allocs,Vector{Allocs}}, + kw... +) + (; edges, scratch) = _get_alloc(allocs) + + # TODO fix bug that requires this to be redefined + edges = Vector{Edge}(undef, 0) + local edge_count = max_ylen = 0 + if tr isa GI.AbstractCurveTrait + edge_count, max_ylen = _to_edges!(edges, geom, dims, edge_count) + else + for ring in GI.getring(geom) + edge_count, ring_max_ylen = _to_edges!(edges, ring, dims, edge_count) + max_ylen = max(max_ylen, ring_max_ylen) + end + end + + # We may have allocated too much + edges1 = view(edges, 1:edge_count) + @static if VERSION < v"1.9-alpha1" + sort!(edges1) + else + sort!(edges1; scratch) + end + + return Edges(edges, max_ylen, edge_count) +end + +Base.parent(edges::Edges) = edges.edges +Base.length(edges::Edges) = edges.edge_count +Base.axes(edges::Edges) = axes(parent(edges)) +Base.getindex(edges::Edges, I...) = getindex(parent(edges), I...) +Base.setindex(edges::Edges, x, I...) = setindex!(parent(edges), x, I...) + +@noinline function _to_edges!(edges, geom, dims, edge_count) + # Dummy Initialisation + local firstpos = prevpos = nextpos = Position((0.0, 0.0), 0) + isfirst = true + local max_ylen = 0 + + GI.npoint(geom) > 0 || return edge_count, max_ylen + xlookup, ylookup = lookup(dims, (X(), Y())) + (length(xlookup) > 0 && length(ylookup) > 0) || return edge_count, max_ylen + + # Raster properties + starts = (Float64(first(xlookup)), Float64(first(ylookup))) + steps = (Float64(Base.step(xlookup)), Float64(Base.step(ylookup))) + local prevpoint = (0.0, 0.0) + + # Loop over points to generate edges + for point in GI.getpoint(geom) + p = (Float64(GI.x(point)), Float64(GI.y(point))) + + # For the first point just set variables + if isfirst + prevpos = firstpos = Position(p, starts, steps) + prevpoint = p + isfirst = false + continue + end + + # Get the next offsets and indices + nextpos = Position(p, starts, steps) + + # Check if we need an edge between these offsets + # This is the performance-critical step that reduces the size of the edge list + if can_skip(prevpos, nextpos, xlookup, ylookup) + prevpos = nextpos + continue + end + + # Add the edge to our `edges` vector + edge_count += 1 + edge = Edge(prevpos, nextpos) + add_edge!(edges, edge, edge_count) + max_ylen = max(max_ylen, edge.iystop - edge.iystart) + prevpos = nextpos + prevpoint = p + end + # Check in case the polygon is not closed + if prevpos != firstpos + edge_count += 1 + edge = Edge(prevpos, firstpos) + max_ylen = max(max_ylen, edge.iystop - edge.iystart) + add_edge!(edges, edge, edge_count) + end + + return edge_count, max_ylen +end + +function add_edge!(edges, edge, edge_count) + if edge_count <= lastindex(edges) + @inbounds edges[edge_count] = edge + else + push!(edges, edge) + end +end + + +# _burn_geometry! +# Fill a raster with `fill` where it interacts with a geometry. +# This is used in `boolmask` TODO move to mask.jl ? +# +# _istable keyword is a hack so we know not to pay the +# price of calling `istable` which calls `hasmethod` +function burn_geometry!(B::AbstractRaster, data::T; kw...) where T + if Tables.istable(T) + geomcolname = first(GI.geometrycolumns(data))::Symbol + geoms = Tables.getcolumn(data, geomcolname) + _burn_geometry!(B, nothing, geoms; kw...) + else + _burn_geometry!(B, GI.trait(data), data; kw...) + end + return B +end + +# This feature filling is simplistic in that it does not use any feature properties. +# This is suitable for masking. See `rasterize` for a version using properties. +_burn_geometry!(B, obj; kw...) = _burn_geometry!(B, GI.trait(obj), obj; kw...)::Bool +function _burn_geometry!(B::AbstractRaster, ::GI.AbstractFeatureTrait, feature; kw...)::Bool + _burn_geometry!(B, GI.geometry(feature); kw...) +end +function _burn_geometry!(B::AbstractRaster, ::GI.AbstractFeatureCollectionTrait, fc; kw...)::Bool + geoms = (GI.geometry(f) for f in GI.getfeature(fc)) + _burn_geometry!(B, nothing, geoms; kw...) +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, + allocs=_burning_allocs(B; threaded), kw... +)::Bool + range = _geomindices(geoms) + burnchecks = _alloc_burnchecks(range) + if isnothing(collapse) || collapse + _run(range, threaded, progress, "") do i + geom = _getgeom(geoms, i) + ismissing(geom) && return nothing + a = _get_alloc(allocs) + B1 = a.buffer + burnchecks[i] = _burn_geometry!(B1, geom; allocs=a, lock, kw...) + return nothing + end + if allocs isa Allocs + _do_broadcast!(|, B, allocs.buffer) + else + buffers = map(a -> a.buffer, allocs) + _do_broadcast!(|, B, buffers...) + end + else + _run(range, threaded, progress, "") do i + geom = _getgeom(geoms, i) + ismissing(geom) && return nothing + B1 = view(B, Dim{:geometry}(i)) + a = _get_alloc(allocs) + burnchecks[i] = _burn_geometry!(B1, geom; allocs=a, lock, kw...) + return nothing + end + end + + _set_burnchecks(burnchecks, metadata(B), verbose) + return false +end + +function _burn_geometry!(B::AbstractRaster, ::GI.AbstractGeometryTrait, geom; + shape=nothing, verbose=true, boundary=:center, allocs=nothing, kw... +)::Bool + hasburned = false + GI.npoint(geom) > 0 || return false + # Use the specified shape or detect it + shape = shape isa Symbol ? shape : _geom_shape(geom) + if shape === :point + hasburned = _fill_point!(B, geom; fill=true, shape, kw...) + elseif shape === :line + n_on_line = _burn_lines!(B, geom; shape, kw...) + hasburned = n_on_line > 0 + elseif shape === :polygon + # Get the extents of the geometry and array + geomextent = _extent(geom) + arrayextent = Extents.extent(B, DEFAULT_POINT_ORDER) + # Only fill if the gemoetry bounding box overlaps the array bounding box + if !Extents.intersects(geomextent, arrayextent) + verbose && _verbose_extent_info(geomextent, arrayextent) + return false + end + # Take a view of the geometry extent + B1 = view(B, Touches(geomextent)) + buf1 = _init_bools(B1) + # 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] = true + end + end + else + _shape_error(shape) + end + return hasburned +end + +@noinline _shape_error(shape) = + throw(ArgumentError("`shape` is $shape, must be `:point`, `:line`, `:polygon` or `nothing`")) + +@noinline _verbose_extent_info(geomextent, arrayextent) = + @info "A geometry was ignored at $geomextent as it was outside of the supplied extent $arrayextent" + +# _burn_polygon! +# Burn `true` values into a raster +# `boundary` determines how edges are handled +function _burn_polygon!(B::AbstractDimArray, geom; kw...)::Bool + B1 = _prepare_for_burning(B) + _burn_polygon!(B1::AbstractDimArray, GI.geomtrait(geom), geom; kw...) +end +function _burn_polygon!(B::AbstractDimArray, trait, geom; + fill=true, boundary=:center, geomextent, verbose=false, allocs=Allocs(B), kw... +)::Bool + allocs = _get_alloc(allocs) + edges = Edges(geom, dims(B); allocs) + + hasburned::Bool = _burn_polygon!(B, edges, allocs.crossings) + + # Lines + n_on_line = 0 + if boundary !== :center + _check_intervals(B, boundary) + if boundary === :touches + if _check_intervals(B, boundary) + # Add line pixels + n_on_line = _burn_lines!(B, geom; fill)::Int + end + elseif boundary === :inside + if _check_intervals(B, boundary) + # Remove line pixels + n_on_line = _burn_lines!(B, geom; fill=!fill)::Int + end + else + throw(ArgumentError("`boundary` can be :touches, :inside, or :center, got :$boundary")) + end + if verbose + (n_on_line > 0) || @info "$n_on_line pixels were on lines" + end + end + + hasburned |= (n_on_line > 0) + + return hasburned +end +function _burn_polygon!(A::AbstractDimArray, edges::Edges, crossings::Vector{Float64}; + offset=nothing, verbose=true +)::Bool + local prev_ypos = 0 + hasburned = false + # Loop over each index of the y axis + for iy in axes(A, YDim) + # Calculate where on the x axis iy is crossed + ncrossings, prev_ypos = _set_crossings!(crossings, edges, iy, prev_ypos) + # Burn between alternate crossings + status = _burn_crossings!(A, crossings, ncrossings, iy) + hasburned |= status.hasburned + end + return hasburned +end + +function _set_crossings!(crossings::Vector{Float64}, edges::Edges, iy::Int, prev_ypos::Int) + # max_ylen tells us how big the largest y edge is. + # We can use this to jump back from the last y position + # rather than iterating from the start of the edges + ypos = max(1, prev_ypos - edges.max_ylen - 1) + ncrossings = 0 + # We know the maximum size on y, so we can start from ypos + start_ypos = searchsortedfirst(edges, ypos) + prev_ypos = start_ypos + for i in start_ypos:lastindex(edges) + e = @inbounds edges[i] + # Edges are sorted on y, so we can skip + # some at the end once they are larger than iy + if iy < e.iystart + prev_ypos = iy + break + end + if iy <= e.iystop + ncrossings += 1 + if ncrossings <= length(crossings) + @inbounds crossings[ncrossings] = Rasters.x_at_y(e, iy) + else + push!(crossings, Rasters.x_at_y(e, iy)) + end + end + end + # For some reason this is much faster than `partialsort!` + sort!(view(crossings, 1:ncrossings)) + return ncrossings, prev_ypos +end + +function _burn_crossings!(A, crossings, ncrossings, iy; + status::BurnStatus=BurnStatus() +) + stop = false + # Start burning loop from outside any rings + (; ic, burn) = status + ix = firstindex(A, X()) + hasburned = false + while ic <= ncrossings + crossing = crossings[ic] + # Burn/skip until we hit the next edge crossing + while ix < crossing + if ix > lastindex(A, X()) + stop = true + break + end + if burn + @inbounds A[X(ix), Y(iy)] = true + hasburned = true + end + ix += 1 + end + if stop + break + else + # Alternate burning/skipping with each edge crossing + burn = !burn + ic += 1 + end + end + # Maybe fill in the end of the row + if burn + for x in ix:lastindex(A, X()) + @inbounds A[X(ix), Y(iy)] = true + end + end + return BurnStatus(ic, burn, hasburned) +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 acheive 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) + +_chki(B) = all(map(s -> s isa Intervals, sampling(dims(B)))) + +function _prepare_for_burning(B, locus=Center()) + B1 = _forward_ordered(B) + start_dims = map(dims(B1, DEFAULT_POINT_ORDER)) do d + # Shift lookup values to center of pixels + d = DD.maybeshiftlocus(locus, d) + _lookup_as_array(d) + end + return setdims(B1, start_dims) +end + +# Convert to Array if its not one already +_lookup_as_array(d::Dimension) = parent(lookup(d)) isa Array ? d : modify(Array, d) + +function _forward_ordered(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))) + set(A, d => reverse(d)) + else + A + end + end +end + + +# _fill_point! +# Fill a raster with `fill` where points are inside raster pixels +@noinline _fill_point!(x::RasterStackOrArray, geom; kw...) = _fill_point!(x, GI.geomtrait(geom), geom; kw...) +@noinline function _fill_point!(x::RasterStackOrArray, ::GI.AbstractGeometryTrait, geom; kw...) + # Just find which pixels contain the points, and set them to true + _without_mapped_crs(x) do x1 + for point in GI.getpoint(geom) + _fill_point!(x, point; kw...) + end + end + return true +end +@noinline function _fill_point!(x::RasterStackOrArray, ::GI.AbstractPointTrait, point; + fill, atol=nothing, lock=nothing, kw... +) + dims1 = commondims(x, DEFAULT_POINT_ORDER) + selectors = map(dims1) do d + _at_or_contains(d, _dimcoord(d, point), atol) + end + # TODO make a check in dimensionaldata that returns the index if it is inbounds + if hasselection(x, selectors) + I = dims2indices(dims1, selectors) + if isnothing(lock) + _fill_index!(x, fill, I) + else + sector = CartesianIndices(map(i -> i:i, I)) + Base.lock(lock, sector) + _fill_index!(x, fill, I) + Base.unlock(lock) + end + return true + else + return false + end +end + +# Fill Int indices directly +_fill_index!(st::AbstractRasterStack, fill::NamedTuple, I::NTuple{<:Any,Int}) = st[I...] = fill +_fill_index!(A::AbstractRaster, fill, I::NTuple{<:Any,Int}) = A[I...] = fill +_fill_index!(A::AbstractRaster, fill::Function, I::NTuple{<:Any,Int}) = A[I...] = fill(A[I...]) + +_fill_index!(st::AbstractRasterStack, fill::NamedTuple, I) = + map((A, f) -> A[I...] .= Ref(f), st, fill) +_fill_index!(A::AbstractRaster, fill, I) = A[I...] .= Ref(fill) +_fill_index!(A::AbstractRaster, fill::Function, I) = A[I...] .= fill.(view(A, I...)) + +# _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) + B1 = _prepare_for_burning(B) + return _burn_lines!(B1, geom, fill) +end + +_burn_lines!(B, geom, fill) = + _burn_lines!(B, GI.geomtrait(geom), geom, fill) +function _burn_lines!(B::AbstractArray, ::Union{GI.MultiLineStringTrait}, geom, fill) + n_on_line = 0 + for linestring in GI.getlinestring(geom) + n_on_line += _burn_lines!(B, linestring, fill) + end + return n_on_line +end +function _burn_lines!( + B::AbstractArray, ::Union{GI.MultiPolygonTrait,GI.PolygonTrait}, geom, fill +) + n_on_line = 0 + for ring in GI.getring(geom) + n_on_line += _burn_lines!(B, ring, fill) + end + return n_on_line +end +function _burn_lines!( + B::AbstractArray, ::GI.AbstractCurveTrait, linestring, fill +) + isfirst = true + local firstpoint, laststop + n_on_line = 0 + for point in GI.getpoint(linestring) + if isfirst + isfirst = false + firstpoint = point + laststop = (x=GI.x(point), y=GI.y(point)) + continue + end + if point == firstpoint + isfirst = true + end + line = ( + start=laststop, + stop=(x=GI.x(point), y=GI.y(point)), + ) + laststop = line.stop + n_on_line += _burn_line!(B, line, fill) + end + return n_on_line +end +function _burn_lines!( + B::AbstractArray, t::GI.LineTrait, line, fill +) + 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) +end + +# _burn_line! +# +# Line-burning algorithm +# Burns a single line into a raster with value where pixels touch a line +# +# 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)) + + @assert order(xdim) == order(ydim) == LookupArrays.ForwardOrdered() + @assert locus(xdim) == locus(ydim) == LookupArrays.Center() + + raster_x_step = abs(step(span(A, X))) + raster_y_step = abs(step(span(A, Y))) + 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 + + # TODO merge this with Edge generation + # Converted lookup to array axis values (still floating) + relstart = (x=(line.start.x - raster_x_offset) / raster_x_step, + y=(line.start.y - raster_y_offset) / raster_y_step) + relstop = (x=(line.stop.x - raster_x_offset) / raster_x_step, + y=(line.stop.y - raster_y_offset) / raster_y_step) + + # Ray/Slope calculations + # Straight distance to the first vertical/horizontal grid boundaries + if relstop.x > relstart.x + xoffset = trunc(relstart.x) - relstart.x + 1 + xmoves = trunc(Int, relstop.x) - trunc(Int, relstart.x) + else + xoffset = relstart.x - trunc(relstart.x) + xmoves = trunc(Int, relstart.x) - trunc(Int, relstop.x) + end + if relstop.y > relstart.y + yoffset = trunc(relstart.y) - relstart.y + 1 + ymoves = trunc(Int, relstop.y) - trunc(Int, relstart.y) + else + yoffset = relstart.y - trunc(relstart.y) + ymoves = trunc(Int, relstart.y) - trunc(Int, relstop.y) + end + manhattan_distance = xmoves + ymoves + + # 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)) + + if manhattan_distance == 0 + D = map((d, o) -> d(o), dimconstructors, (j, i)) + if checkbounds(Bool, A, D...) + @inbounds A[D...] = fill + end + n_on_line = 1 + return n_on_line + end + + diff_x = relstop.x - relstart.x + diff_y = relstop.y - relstart.y + + # Angle of ray/slope. + # max: How far to move along the ray to cross the first cell boundary. + # delta: How far to move along the ray to move 1 grid cell. + hyp = @fastmath sqrt(diff_y^2 + diff_x^2) + cs = diff_x / hyp + si = -diff_y / hyp + + delta_x, max_x = if isapprox(cs, zero(cs); atol=1e-10) + -Inf, Inf + else + 1.0 / cs, xoffset / cs + end + delta_y, max_y = if isapprox(si, zero(si); atol=1e-10) + -Inf, Inf + else + 1.0 / si, yoffset / si + end + # Count how many exactly hit lines + 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 + end + + # Only move in either X or Y coordinates, not both. + if abs(max_x) < abs(max_y) + max_x += delta_x + j += step_j + countx +=1 + else + max_y += delta_y + i += step_i + 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)) +end + +# Get the GeoInterface coord from a point for a specific Dimension +_dimcoord(::XDim, point) = GI.x(point) +_dimcoord(::YDim, point) = GI.y(point) +_dimcoord(::ZDim, point) = GI.z(point) + +# Get the shape category for a geometry +@inline _geom_shape(geom) = _geom_shape(GI.geomtrait(geom), geom) +@inline _geom_shape(::Union{<:GI.PointTrait,<:GI.MultiPointTrait}, geom) = :point +@inline _geom_shape(::Union{<:GI.LineTrait,<:GI.LineStringTrait,<:GI.MultiLineStringTrait}, geom) = :line +@inline _geom_shape(::Union{<:GI.LinearRingTrait,<:GI.PolygonTrait,<:GI.MultiPolygonTrait}, geom) = :polygon +@inline _geom_shape(x, geom) = throw(ArgumentError("Geometry trait $x cannot be rasterized"), geom) +@inline _geom_shape(::Nothing, geom) = throw(ArgumentError("Object is not a GeoInterface.jl compatible geometry: $geom"), geom) + + +# 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, 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::Extents.Extent, T::Type, data; kw...) = _init_bools(to, _extent2dims(to; kw...), 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; kw...) + # Get the extent of the geometries + ext = _extent(data) + isnothing(ext) && throw(ArgumentError("no recognised dimensions, extent or geometry")) + # Convert the extent to dims (there must be `res` or `size` in `kw`) + dims = _extent2dims(ext; kw...) + return _init_bools(to, dims, T, data; kw...) +end +function _init_bools(to, dims::DimTuple, T::Type, data; collapse::Union{Bool,Nothing}=nothing, kw...) + if isnothing(data) || isnothing(collapse) || collapse + _alloc_bools(to, dims, T; kw...) + else + n = if Base.IteratorSize(data) isa Base.HasShape + length(data) + else + count(_ -> true, data) + end + geomdim = Dim{:geometry}(1:n) + _alloc_bools(to, (dims..., geomdim), T; kw...) + 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 +end +function _alloc_bools(to, dims::DimTuple, ::Type{T}; missingval=false, metadata=NoMetadata(), kw...) where T + # Use an `Array` + data = fill!(Raster{T}(undef, dims), missingval) + return rebuild(data; missingval, metadata) +end + +_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) +end +_set_burnchecks(burnchecks, metadata, verbose) = verbose && _burncheck_info(burnchecks) +function _burncheck_info(burnchecks) + nburned = sum(burnchecks) + nmissed = length(burnchecks) - nburned + nmissed > 0 && @info "$nmissed geometries did not affect any pixels. See `metadata(raster)[\"missed_geometries\"]` for a vector of misses" +end + +_nthreads() = Threads.nthreads() + +function _at_or_contains(d, v, atol) + selector = sampling(d) isa Intervals ? Contains(v) : At(v; atol=atol) + DD.basetypeof(d)(selector) +end diff --git a/src/utils.jl b/src/utils.jl index b91671166..f180331e0 100644 --- a/src/utils.jl +++ b/src/utils.jl @@ -183,27 +183,6 @@ function _as_intervals(ds::Tuple) return setdims(ds, interval_dims) end -_geomindices(geoms) = _geomindices(GI.trait(geoms), geoms) -_geomindices(::Nothing, geoms) = eachindex(geoms) -_geomindices(::GI.FeatureCollectionTrait, geoms) = 1:GI.nfeature(geoms) -_geomindices(::GI.FeatureTrait, geoms) = _geomindices(GI.geometry(geoms)) -_geomindices(::GI.AbstractGeometryTrait, geoms) = 1:GI.ngeom(geoms) - -_getgeom(geoms, i::Integer) = _getgeom(GI.trait(geoms), geoms, i) -_getgeom(::GI.FeatureCollectionTrait, geoms, i::Integer) = GI.geometry(GI.getfeature(geoms, i)) -_getgeom(::GI.FeatureTrait, geoms, i::Integer) = GI.getgeom(GI.geometry(geoms), i) -_getgeom(::GI.AbstractGeometryTrait, geoms, i::Integer) = GI.getgeom(geom, i) -_getgeom(::GI.PointTrait, geom, i::Integer) = error("PointTrait should not be reached") -_getgeom(::Nothing, geoms, i::Integer) = geoms[i] # Otherwise we can probably just index? - -_getgeom(geoms) = _getgeom(GI.trait(geoms), geoms) -_getgeom(::GI.FeatureCollectionTrait, geoms) = (GI.geometry(f) for f in GI.getfeature(geoms)) -_getgeom(::GI.FeatureTrait, geoms) = GI.getgeom(GI.geometry(geoms)) -_getgeom(::GI.AbstractGeometryTrait, geoms) = GI.getgeom(geoms) -_getgeom(::GI.PointTrait, geom) = error("PointTrait should not be reached") -_getgeom(::Nothing, geoms) = geoms - - _warn_disk() = @warn "Disk-based objects may be very slow here. User `read` first." _filenotfound_error(filename) = throw(ArgumentError("file \"$filename\" not found")) @@ -242,80 +221,3 @@ function _run(f, range::OrdinalRange, threaded::Bool, progress::Bool, desc::Stri end end end - - -const XYExtent = Extents.Extent{(:X,:Y),Tuple{Tuple{Float64,Float64},Tuple{Float64,Float64}}} - -# Get the bounds of a geometry -_extent(geom; kw...)::XYExtent = _extent(GI.trait(geom), geom; kw...) -function _extent(::Nothing, data::AbstractVector; kw...)::XYExtent - g1 = first(data) - if GI.trait(g1) isa GI.PointTrait - xs = extrema(p -> GI.x(p), data) - ys = extrema(p -> GI.y(p), data) - return _float64_xy_extent(Extents.Extent(X=xs, Y=ys)) - else - ext = reduce(data; init=_extent(first(data))) do ext, geom - Extents.union(ext, _extent(geom)) - end - return _float64_xy_extent(ext) - end -end -_extent(::Nothing, data::RasterStackOrArray; kw...)::XYExtent = _float64_xy_extent(Extents.extent(data)) -function _extent(::Nothing, data::T; geometrycolumn=nothing)::XYExtent where T - if Tables.istable(T) - singlecolumn = isnothing(geometrycolumn) ? first(GI.geometrycolumns(data)) : geometrycolumn - cols = Tables.columns(data) - if singlecolumn isa Symbol && singlecolumn in Tables.columnnames(cols) - # Table of geometries - geoms = Tables.getcolumn(data, singlecolumn) - return _extent(nothing, geoms) - else - multicolumn = isnothing(geometrycolumn) ? DEFAULT_POINT_ORDER : geometrycolumn - # TODO: test this branch - # Table of points with dimension columns - bounds = reduce(multicolumn; init=(;)) do acc, key - if key in Tables.columnnames(cols) - merge(acc, (; key=extrema(cols[key]))) - else - acc - end - end - return _float64_xy_extent(Extensts.Extent(bounds)) - end - else - ext = Extents.extent(data) - ext isa Extents.Extent || throw(ArgumentError("object returns `nothing` from `Extents.extent`.")) - return _float64_xy_extent(ext) - end -end -function _extent(::GI.AbstractPointTrait, point; kw...)::XYExtent - x, y = Float64(GI.x(point)), Float64(GI.y(point)) - Extents.Extent(X=(x, x), Y=(y, y)) -end -function _extent(::GI.AbstractGeometryTrait, geom; kw...)::XYExtent - geomextent = GI.extent(geom; fallback=false) - if isnothing(geomextent) - points = GI.getpoint(geom) - xbounds = extrema(GI.x(p) for p in points) - ybounds = extrema(GI.y(p) for p in points) - return _float64_xy_extent(Extents.Extent(X=xbounds, Y=ybounds)) - else - return _float64_xy_extent(geomextent) - end -end -_extent(::GI.AbstractFeatureTrait, feature; kw...)::XYExtent = _extent(GI.geometry(feature)) -function _extent(::GI.AbstractFeatureCollectionTrait, features; kw...)::XYExtent - features = GI.getfeature(features) - init = _float64_xy_extent(_extent(first(features))) - ext = reduce(features; init) do acc, f - Extents.union(acc, _extent(f)) - end - return _float64_xy_extent(ext) -end - -function _float64_xy_extent(ext::Extents.Extent) - xbounds = map(Float64, ext.X) - ybounds = map(Float64, ext.Y) - return Extents.Extent(X=xbounds, Y=ybounds) -end diff --git a/test/rasterize.jl b/test/rasterize.jl index 9d22e0861..8c166524e 100644 --- a/test/rasterize.jl +++ b/test/rasterize.jl @@ -10,11 +10,13 @@ B = [1.0 0.4; 2.0 missing] ga = Raster(A, (X, Y); missingval=missing) st = RasterStack((a=A, b=B), (X, Y); missingval=(a=missing,b=missing)) -pointvec = [(-20.0, 30.0), - (-20.0, 10.0), - (0.0, 10.0), - (0.0, 30.0), - (-20.0, 30.0)] +pointvec = [ + (-20.0, 30.0), + (-20.0, 10.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]]) @@ -23,7 +25,7 @@ multi_point = ArchGDAL.createmultipoint(pointvec) linestring = ArchGDAL.createlinestring(pointvec) multi_linestring = ArchGDAL.createmultilinestring([pointvec]) linearring = ArchGDAL.createlinearring(pointvec) -line_collection = GeoInterface.GeometryCollection([linestring, multi_linestring]) +line_collection = GeoInterface.GeometryCollection([multi_linestring]) poly_collection = GeoInterface.GeometryCollection([polygon]) pointfc = map(GeoInterface.getpoint(polygon), vals) do geom, v (geometry=geom, val1=v, val2=2.0f0v) @@ -44,8 +46,10 @@ st = RasterStack((A1, copy(A1))) geom = polygon geom = linearring geom = pointvec + 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) + for A in (A1, A2), geom in (pointvec, pointfc, multi_point, linestring, multi_linestring, linearring, polygon, multi_polygon, table, line_collection, poly_collection) fill!(A, 0) rasterize!(sum, A, geom; shape=:point, fill=1); @test sum(A) == 5.0 @@ -61,30 +65,36 @@ st = RasterStack((A1, copy(A1))) 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) + 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)) st.layer1 .= st.layer2 .= 0 rasterize!(sum, st, geom; shape=:point, fill=(layer1=2, layer2=3)) - @test sum(st[:layer1]) == 10 - @test sum(st[:layer2]) == 15 - @test parent(st[:layer1]) isa Array{Float64,2} + @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)) + @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)) + @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) st[:layer1] .= 0; st[:layer2] .= 0 rasterize!(sum, st, geom; shape=:point, fill=(layer1=1:5, layer2=6:10)) @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=2, layer2=3)) - @test sum(st[:layer1]) == 8 - @test sum(st[:layer2]) == 12 rasterize!(last, st, geom; shape=:point, fill=(layer1=1:5, layer2=6:10)) @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=2, layer2=3)) - @test sum(st[:layer1]) == 8 - @test sum(st[:layer2]) == 12 - st[:layer1] .= 0; st[:layer2] .= 0 rasterize!(first, st, geom; shape=:point, fill=(layer1=1:5, layer2=6:10)) @test sum(st[:layer1]) == sum(1:4) @test sum(st[:layer2]) == sum(6:9) @@ -96,18 +106,19 @@ end @testset "all line and polygon geoms work as :line" begin A = A1 geom = linestring - for A in (A1, A2), geom in (linestring, multi_linestring, linearring, polygon, multi_polygon, line_collection) + geom = line_collection + for A in (A1, A2), geom in (linestring, multi_linestring, linearring, polygon, multi_polygon, line_collection, poly_collection) A .= 0 rasterize!(sum, A, geom; shape=:line, fill=1) @test sum(A) == 20 + 20 + 20 + 20 @test sum(rasterize(sum, geom; to=A, shape=:line, fill=1, missingval=0)) == 80 end @testset ":line is detected for line geometries" begin - for A in (A1, A2), geom in (linestring, multi_linestring, line_collection) + for A in (A1, A2), geom in (linestring, multi_linestring) A .= 0 rasterize!(A, geom; fill=1) @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)) == 80 end end end @@ -115,29 +126,30 @@ end @testset "polygon geoms work as :polygon" begin A = A1 poly = polygon + poly = poly_collection for A in (A1, A2), poly in (polygon, multi_polygon, poly_collection) 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) @test parent(ra) isa Matrix{Int} @test sum(ra) == sum(ra_res) === 20 * 20 - ra = rasterize(poly; to=A, shape=:polygon, fill=1, boundary=:touches) + ra = rasterize(last, poly; to=A, shape=:polygon, fill=1, boundary=:touches) @test parent(ra) isa Array{Union{Missing,Int},2} @test sum(skipmissing(ra)) === 21 * 21 - rasterize!(A, poly; shape=:polygon, fill=1, boundary=:inside) + rasterize!(last, A, poly; shape=:polygon, fill=1, boundary=:inside) @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(poly; to=A, fill=1) + R = rasterize(last, poly; to=A, fill=1) @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) end - st1 = rasterize(poly; fill=(layer1=1, layer2=2), to=st) + st1 = rasterize(last, poly; fill=(layer1=1, layer2=2), to=st) @test sum(skipmissing(st1[:layer1])) == 400 # The last value overwrites the first @test sum(skipmissing(st1[:layer2])) == 800 # Missing size / res