mask fails when geometry encompasses raster

[asinghvi17]

MWE:

using Rasters, GeoInterface
geom = GeoInterface.Polygon([GeoInterface.LinearRing([(0.0, 0.0), (0.0, 10.0), (10.0, 10.0), (10.0, 0.0), (0.0, 0.0)])])
raster = Raster(ones(11, 11), (X(1:0.1:2), Y(1:0.1:2)))

mask(raster, with=geom)

╭─────────────────────────────────────────╮
│ 11×11 Raster{Union{Missing, Float64},2} │
├─────────────────────────────────────────┴───────────────────────────────────────────────────────────────────────────────────────────────────── dims ┐
  ↓ X Sampled{Float64} 1.0:0.1:2.0 ForwardOrdered Regular Points,
  → Y Sampled{Float64} 1.0:0.1:2.0 ForwardOrdered Regular Points
├───────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────── raster ┤
  extent: Extent(X = (1.0, 2.0), Y = (1.0, 2.0))
  missingval: missing
└─────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────┘
 ↓ →  1.0       1.1       1.2       1.3       1.4       1.5       1.6       1.7       1.8       1.9       2.0
 1.0   missing   missing   missing   missing   missing   missing   missing   missing   missing   missing   missing
 1.1   missing   missing   missing   missing   missing   missing   missing   missing   missing   missing   missing
 1.2   missing   missing   missing   missing   missing   missing   missing   missing   missing   missing   missing
 1.3   missing   missing   missing   missing   missing   missing   missing   missing   missing   missing   missing
 1.4   missing   missing   missing   missing   missing   missing   missing   missing   missing   missing   missing
 1.5   missing   missing   missing   missing   missing   missing   missing   missing   missing   missing   missing
 1.6   missing   missing   missing   missing   missing   missing   missing   missing   missing   missing   missing
 1.7   missing   missing   missing   missing   missing   missing   missing   missing   missing   missing   missing
 1.8   missing   missing   missing   missing   missing   missing   missing   missing   missing   missing   missing
 1.9   missing   missing   missing   missing   missing   missing   missing   missing   missing   missing   missing
 2.0   missing   missing   missing   missing   missing   missing   missing   missing   missing   missing   missing

that seems a bit strange?

Even rasterize has a similarly weird result:

 rasterize(geom, to=raster, fill=1)

╭───────────────────────────────────────────────╮
│ 11×11 Raster{Union{Missing, Int64},2} nothing │
├───────────────────────────────────────────────┴────────────────────────────────────────────────────────────────────────────────────────────────────────── dims ┐
  ↓ X Sampled{Float64} 1.0:0.1:2.0 ForwardOrdered Regular Intervals{Center},
  → Y Sampled{Float64} 1.0:0.1:2.0 ForwardOrdered Regular Intervals{Center}
├────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────── metadata ┤
  Metadata of Dict{Any, Any}()
├──────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────── raster ┤
  extent: Extent(X = (0.95, 2.05), Y = (0.95, 2.05))
  missingval: missing
└────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────┘
 ↓ →  1.0       1.1       1.2       1.3       1.4       1.5       1.6       1.7       1.8       1.9       2.0
 1.0   missing   missing   missing   missing   missing   missing   missing   missing   missing   missing   missing
 1.1   missing   missing   missing   missing   missing   missing   missing   missing   missing   missing   missing
 1.2   missing   missing   missing   missing   missing   missing   missing   missing   missing   missing   missing
 ⋮                                            ⋮                                                 ⋮
 1.7   missing   missing   missing   missing   missing   missing   missing   missing   missing   missing   missing
 1.8   missing   missing   missing   missing   missing   missing   missing   missing   missing   missing   missing
 1.9   missing   missing   missing   missing   missing   missing   missing   missing   missing   missing   missing
 2.0   missing   missing   missing   missing   missing   missing   missing   missing   missing   missing   missing

[asinghvi17]

This also causes zonal to fail, which is where I found this bug. If the geometry intersects but does not cover the raster, it works correctly.

[rafaqz]

Yeah, that's a bug. Something wrong in the polygon burning algorithm

[asinghvi17]

Maybe

Rasters.jl/src/methods/burning/edges.jl

Lines 89 to 91 in 251bbbd

	# 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

? That's the only suspicious code I found so far

[asinghvi17]

Trying to troubleshoot further. I can construct the Edges object manually after some sleuthing:

allocs = Rasters._burning_allocs(raster; threaded = false)
Rasters.Edges(geom, dims(raster); allocs)

which yields

2-element Rasters.Edges:
 Rasters.Edge((-9.0, -9.0), 0.0, -8, 91)
 Rasters.Edge((91.0, -9.0), 0.0, -8, 91)

Not sure what this means though. The fields are:

  start    :: Tuple{Float64, Float64}
  gradient :: Float64
  iystart  :: Int32
  iystop   :: Int32

so iystart and iystop both appear to encompass the raster, and seem broadly correct as well.

[asinghvi17]

to dig in further,

receptacle = fill(false, dims(raster))
Rasters._burn_polygon!(receptacle, geom; geomextent = GeoInterface.extent(geom))
receptacle

yields

julia> receptacle
╭────────────────────────╮
│ 11×11 DimArray{Bool,2} │
├────────────────────────┴─────────────────────────────────────────────────────────────────────────────────────────────── dims ┐
  ↓ X Sampled{Float64} 1.0:0.1:2.0 ForwardOrdered Regular Points,
  → Y Sampled{Float64} 1.0:0.1:2.0 ForwardOrdered Regular Points
└──────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────┘
 ↓ →  1.0  1.1  1.2  1.3  1.4  1.5  1.6  1.7  1.8  1.9  2.0
 1.0  0    0    0    0    0    0    0    0    0    0    0
 1.1  0    0    0    0    0    0    0    0    0    0    0
 1.2  0    0    0    0    0    0    0    0    0    0    0
 1.3  0    0    0    0    0    0    0    0    0    0    0
 ⋮                        ⋮                        ⋮    
 1.6  0    0    0    0    0    0    0    0    0    0    0
 1.7  0    0    0    0    0    0    0    0    0    0    0
 1.8  0    0    0    0    0    0    0    0    0    0    0
 1.9  0    0    0    0    0    0    0    0    0    0    0
 2.0  0    0    0    0    0    0    0    0    0    0    0

when it should presumably be filled with ones.

[rafaqz]

Yeah, I'm guessing there is something wrong with the algorithm for this case, where it never hits a y crossing at all so it's never triggered to fill up the row with true values. Or something.