Wrong coordinates returned when clipping with exact coordinates

[gabrieldansereau]

@tpoisot I ran into the following issue with clip while trying to subset a layer for Canada into smaller subregions. The boundary returned for some exact coordinates (e.g. 56.0 as a top boundary) sometimes leaves out one pixel row/column. I'm running into this issue with SpeciesDistributionToolkit v0.0.2 in a temporary environment, but I also had it with SimpleSDMLayers v0.8.3.

Here is an example:

using SpeciesDistributionToolkit

spatialrange = (left = -145.0, right = -50.0, bottom = 40.0, top = 89.0);
layer = SimpleSDMPredictor(RasterData(WorldClim2, BioClim); spatialrange...); # also happens with 2.5 arcmin

# Wrong coordinates while clipping
clip(layer; top=56.0).top # 55.833333333333336
clip(layer; top=56.0 - 0.0000001).top # 56.0
clip(layer; top=56.0 + 0.0000001).top # 56.16666666666667

I would expect the first two clip calls to return 56.0 exactly, otherwise it leaves out one cell row. I tried with other exact coordinates and the problem seems to happen with no specific pattern:

julia> exact_lats = (spatialrange.bottom+1.0):1.0:(spatialrange.top - 1.0)
41.0:1.0:88.0

julia> show([clip(layer; top=l).top for l in exact_lats]) # patterns seems random...
[40.833333333333336, 41.833333333333336, 43.0, 43.833333333333336, 45.0, 46.0, 47.0, 47.833333333333336, 49.0, 50.0, 51.0, 52.0, 53.0, 54.0, 55.0, 55.833333333333336, 57.0, 58.0, 59.0, 60.0, 61.0, 62.0, 63.0, 64.0, 65.0, 66.0, 66.83333333333333, 68.0, 69.0, 70.0, 71.0, 71.83333333333333, 73.0, 74.0, 75.0, 76.0, 76.83333333333333, 78.0, 79.0, 80.0, 81.0, 82.0, 83.0, 84.0, 85.0, 86.0, 87.0, 88.0]

The problem also happens with the bottom bound and with longitudes. On the other hand, it does not happen with geotiff reading calls.

julia> exact_lats = (spatialrange.bottom+1.0):1.0:(spatialrange.top - 1.0)
41.0:1.0:88.0

julia> exact_lons = (spatialrange.left+1.0):1.0:(spatialrange.right - 1.0)
-144.0:1.0:-51.0

julia> all(isinteger.([clip(layer; top=l).top for l in exact_lats])) # patterns seems random...
false

julia> all(isinteger.([clip(layer; bottom=l).bottom for l in exact_lats])) # only 1 case
false

julia> all(isinteger.([clip(layer; right=l).right for l in exact_lons])) # happening with longitudes too
false

julia> all(isinteger.([clip(layer; left=l).left for l in exact_lons]))
false

julia> all(isinteger.([SimpleSDMPredictor(RasterData(WorldClim2, BioClim); left = -145.0, right = -50.0, bottom = 40.0, top = t).top for t in exact_lats]))
true

I tracked the problem down to the _match_latitude function. It seems to be a rounding issue. I have a fix I'll make in a PR and let you review.

SpeciesDistributionToolkit.jl/SimpleSDMLayers/src/lib/coordinateconversion.jl

Line 16 in 965c23c

function _match_latitude(layer::T, lat::K; lower::Bool=true) where {T <: SimpleSDMLayer, K <: AbstractFloat}