Reading certain tiffs with geotiff is wrong if the CRS is not WGS84

[gottacatchenall]

While working on this, GEO-BON/bon-in-a-box-pipelines#10, we realized that some types of missing values throw an error when trying to load

What happened?

Reading layers in with geotiff, we run into several issues with missing/NaN values in geotiffs.

As an example

layerpath = "chelsa_clim_bio11981.tif"
geotiff(SimpleSDMPredictor, layerpath)

fails with the below error where the chelsa tiff is available here

Stacktrace

ERROR: MethodError: no method matching iterate(::Nothing)
Closest candidates are:
  iterate(::Union{LinRange, StepRangeLen}) at range.jl:872
  iterate(::Union{LinRange, StepRangeLen}, ::Integer) at range.jl:872
  iterate(::T) where T<:Union{Base.KeySet{<:Any, <:Dict}, Base.ValueIterator{<:Dict}} at dict.jl:712
  ...
Stacktrace:
 [1] indexed_iterate(I::Nothing, i::Int64)
   @ Base ./tuple.jl:91
 [2] (::SimpleSDMLayers.var"#53#54"{SimpleSDMPredictor, Int64})(dataset::ArchGDAL.Dataset)
   @ SimpleSDMLayers ~/.julia/packages/SimpleSDMLayers/lYDIT/src/datasets/geotiff.jl:81
 [3] read(f::SimpleSDMLayers.var"#53#54"{SimpleSDMPredictor, Int64}, args::String; kwargs::Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}})
   @ ArchGDAL ~/.julia/packages/ArchGDAL/zkx2f/src/context.jl:267
 [4] read
   @ ~/.julia/packages/ArchGDAL/zkx2f/src/context.jl:264 [inlined]
 [5] #geotiff#52
   @ ~/.julia/packages/SimpleSDMLayers/lYDIT/src/datasets/geotiff.jl:46 [inlined]
 [6] geotiff (repeats 2 times)
   @ ~/.julia/packages/SimpleSDMLayers/lYDIT/src/datasets/geotiff.jl:32 [inlined]
 [7] top-level scope
   @ ~/Code/Tmp/foo.jl:9

[gottacatchenall]

i've tracked this down to

https://github.com/EcoJulia/SimpleSDMLayers.jl/blob/dcc49145c2c844bfe63fff99b96cf86720076512/src/datasets/geotiff.jl#L5

in _find_span.

pos is less than m, the function returns nothing, and https://github.com/EcoJulia/SimpleSDMLayers.jl/blob/dcc49145c2c844bfe63fff99b96cf86720076512/src/datasets/geotiff.jl#L81

tries to unpack nothing which gives the resulting stacktrace.

This is because top < minlat in the input args.

[gottacatchenall]

The issue is fixed by changing

https://github.com/EcoJulia/SimpleSDMLayers.jl/blob/07ecfb6495d7bcd186210a77c0781e6fc1397c3d/src/datasets/geotiff.jl#L68
to top = isnothing(top) ? maxlat : max(top, maxlat), but this doesn't really make sense. perhaps something in

minlon = transform[1]
maxlat = transform[4]
maxlon = minlon + size(band,1)*transform[2]
minlat = maxlat - abs(size(band,2)*transform[6])

is wrong? in the example tiff, maxlat is from the ArchGDAL.getgeotransform call is 2.45e6, which doesn't make any sense @tpoisot

[tpoisot]

Yeah, unproperly formatted geotiffs happen (you should see the TerraClimate data). It's weird because the CHELSA data are usually very clean.

I'll take a look - very bad code smell to have a return nothing here, so I'll probably end up refactoring this part.

[gottacatchenall]

It looks like the issue with the specific layers we are using (this and this) that they are in NAD83 but SimpleSDMLayers assumes WGS84 when reading geotiffs in,

[tpoisot]

@gottacatchenall / @karinorman - I am actively trying to fix this today, FYI. There's a reprojection function in ArchGDAL which should do most of the work for us.

[tpoisot]

Here we go