missingval keyword for RasterStack

[tiemvanderdeure]

It doesn't seem possible to pass a missingval keyword to RasterStack at the moment. Would this be hard to implement?

I ran into this with the CHELSA dataset, which has UInt16 type but a missingval of -99999, which throws an error unless missingval is passed explicitly (not sure if this is a bug or not in itself?).

E.g.

using RasterDataSources, Rasters
RasterStack(CHELSA{BioClim}, (1,2); missingval = -99999.0)

This workaround does work, of course:
rs = RasterStack(Raster.(CHELSA{BioClim}, (1,2); missingval = -99999.0)...)

[rafaqz]

How can there be any missing values if its below zero? It shouldnt have a missing value set at all. (so set it to nothing)

But the file should still load. (So we will fix that part - but the misssingval will still be wrong)

[rafaqz]

Ok looking at it now its actually Int16, and its has the correct missingval of -32768...

But you are right, its not possible to change it to -9999. Will fix that

julia> st = RasterStack(CHELSA{BioClim}, (1, 2); missingval = -9999)
╭─────────────────────────╮
│ 43200×20880 RasterStack │
├─────────────────────────┴─────────────────────────────────────────────────────────────────────────────────────────────────── dims ┐
  ↓ X Projected{Float64} LinRange{Float64}(-180.00013888885002, 179.99152633785002, 43200) ForwardOrdered Regular Intervals{Start},
  → Y Projected{Float64} LinRange{Float64}(83.99152708185001, -90.00013888884999, 20880) ReverseOrdered Regular Intervals{Start}
├─────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────── layers ┤
  :bio1 eltype: Int16 dims: X, Y size: 43200×20880
  :bio2 eltype: Int16 dims: X, Y size: 43200×20880
├─────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────── raster ┤
  extent: Extent(X = (-180.00013888885002, 179.99985967115), Y = (-90.00013888884999, 83.99986041515001))
  missingval: (bio1 = -32768, bio2 = -32768)
  crs: GEOGCS["WGS 84",DATUM["WGS_1984",SPHEROID["WGS 84",6378137,298.257223563,AUTHORITY["EPSG","7030"]],AUTHORITY["EPSG","6326"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.01745329251994
33,AUTHORITY["EPSG","9122"]],AXIS["Latitude",NORTH],AXIS["Longitude",EAST],AUTHORITY["EPSG","4326"]]
└───────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────┘

julia> st.bio1
╭──────────────────────────────────╮
│ 43200×20880 Raster{Int16,2} bio1 │
├──────────────────────────────────┴────────────────────────────────────────────────────────────────────────────────────────── dims ┐
  ↓ X Projected{Float64} LinRange{Float64}(-180.00013888885002, 179.99152633785002, 43200) ForwardOrdered Regular Intervals{Start},
  → Y Projected{Float64} LinRange{Float64}(83.99152708185001, -90.00013888884999, 20880) ReverseOrdered Regular Intervals{Start}
├───────────────────────────────────────────────────────────────────────────────────────────────────────────────────────── metadata ┤
  Metadata{Rasters.GDALsource} of Dict{String, Any} with 4 entries:
  "units"    => ""
  "offset"   => 0.0
  "filepath" => "/home/raf/Data/Raster/CHELSA/BioClim/CHELSA_bio10_01.tif"
  "scale"    => 1.0
├─────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────── raster ┤
  extent: Extent(X = (-180.00013888885002, 179.99985967115), Y = (-90.00013888884999, 83.99986041515001))
  missingval: -32768
  crs: GEOGCS["WGS 84",DATUM["WGS_1984",SPHEROID["WGS 84",6378137,298.257223563,AUTHORITY["EPSG","7030"]],AUTHORITY["EPSG","6326"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.01745329251994
33,AUTHORITY["EPSG","9122"]],AXIS["Latitude",NORTH],AXIS["Longitude",EAST],AUTHORITY["EPSG","4326"]]
└───────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────┘
    ↓ →        83.9915      83.9832      83.9749      83.9665      83.9582      83.9499      83.9415  …     -89.9501     -89.9585     -89.9668     -89.9751     -89.9835     -89.9918     -90.0001
 -180.0    -32768       -32768       -32768       -32768       -32768       -32768       -32768          -32768       -32768       -32768       -32768       -32768       -32768       -32768
 -179.992  -32768       -32768       -32768       -32768       -32768       -32768       -32768          -32768       -32768       -32768       -32768       -32768       -32768       -32768
    ⋮                                                               ⋮                                 ⋱                                              ⋮                                 
  179.975  -32768       -32768       -32768       -32768       -32768       -32768       -32768          -32768       -32768       -32768       -32768       -32768       -32768       -32768
  179.983  -32768       -32768       -32768       -32768       -32768       -32768       -32768       …  -32768       -32768       -32768       -32768       -32768       -32768       -32768
  179.992  -32768       -32768       -32768       -32768       -32768       -32768       -32768          -32768       -32768       -32768       -32768       -32768       -32768       -32768

[rafaqz]

Ahh it turns out this is a deeper issue of using nothing to mean both "no keyword was passed" and "the obect property has no value".

So we need to pivot to a specific NoKW type as default accross Raster/RasterStack constructors and let users force setting things to nothing. I'll try to do it comprehensively for all the keywords. Its kind of a mess in there...

[tiemvanderdeure]

About the missing values below zero: it seems somethign changed between CHELSA v1 and CHELSA v2. On CHELSA v2, somehow the missingvalue really is -99999, while the type is UInt16, which throws an error.

There aren't actually any missing values though, so setting missingval to whatever value fixes the error.

To reproduce, on this branch of RasterDataSources: EcoJulia/RasterDataSources.jl#62 run:

using RasterDataSources, Rasters, ArchGDAL
Raster(CHELSA{BioClim}, 1; version = 2) # errors
Raster(CHELSA{BioClim}, 1; version = 2, missingval = missing) # works

I don't really understand why missingval is -99999.0. But then the error is because _missingval_from_gdal tries to convert this to UInt16 and there are no checks if this is possible. I'll open a (very tiny) PR to fix this.

It seems a bit inconsistent that a missingval is converted to the datatype of the raster when it's derived, but not when it's set manually, but maybe there are good reasons for that?

Full stacktrace:

ERROR: InexactError: UInt16(-99999.0)
Stacktrace:
  [1] UInt16
    @ .\float.jl:888 [inlined]
  [2] convert
    @ .\number.jl:7 [inlined]
  [3] _missingval_from_gdal(T::Type{UInt16}, x::Float64)
    @ RastersArchGDALExt C:\Users\tsh371\.julia\packages\Rasters\juYeI\ext\RastersArchGDALExt\gdal_source.jl:305
  [4] missingval(::ArchGDAL.RasterDataset{UInt16, ArchGDAL.Dataset})
    @ RastersArchGDALExt C:\Users\tsh371\.julia\packages\Rasters\juYeI\ext\RastersArchGDALExt\gdal_source.jl:269
  [5] (::Rasters.var"#47#48"{@Kwargs{crs::Nothing}, String})(ds::ArchGDAL.RasterDataset{UInt16, ArchGDAL.Dataset})
    @ Rasters C:\Users\tsh371\.julia\packages\Rasters\juYeI\src\array.jl:254
  [6] call_composed
    @ .\operators.jl:1045 [inlined]
  [7] call_composed
    @ .\operators.jl:1044 [inlined]
  [8] (::ComposedFunction{…})(x::ArchGDAL.RasterDataset{…}; kw::@Kwargs{})
    @ Base .\operators.jl:1041
  [9] readraster(f::ComposedFunction{typeof(Rasters.cleanreturn), Rasters.var"#47#48"{…}}, args::String; kwargs::@Kwargs{})
    @ ArchGDAL C:\Users\tsh371\.julia\packages\ArchGDAL\lgE4A\src\context.jl:268
 [10] readraster(f::Function, args::String)
    @ ArchGDAL C:\Users\tsh371\.julia\packages\ArchGDAL\lgE4A\src\context.jl:265
 [11] _open(f::Function, ::Type{Rasters.GDALsource}, filename::String; write::Bool, kw::@Kwargs{})
    @ RastersArchGDALExt C:\Users\tsh371\.julia\packages\Rasters\juYeI\ext\RastersArchGDALExt\gdal_source.jl:113
 [12] _open(f::Function, ::Type{Rasters.GDALsource}, filename::String)
    @ RastersArchGDALExt C:\Users\tsh371\.julia\packages\Rasters\juYeI\ext\RastersArchGDALExt\gdal_source.jl:92
 [13] #_open#136
    @ C:\Users\tsh371\.julia\packages\Rasters\juYeI\src\sources\sources.jl:78 [inlined]
 [14] _open(f::Function, T::Type, filename::String)
    @ Rasters C:\Users\tsh371\.julia\packages\Rasters\juYeI\src\sources\sources.jl:77 [inlined]    
 [15] Raster(filename::String; name::Symbol, key::Symbol, source::Nothing, kw::@Kwargs{crs::Nothing})
    @ Rasters C:\Users\tsh371\.julia\packages\Rasters\juYeI\src\array.jl:252
 [16] Raster(T::Type{CHELSA{BioClim}}, layer::Int64; crs::Nothing, kw::@Kwargs{version::Int64})    
    @ RastersRasterDataSourcesExt C:\Users\tsh371\.julia\packages\Rasters\juYeI\ext\RastersRasterDataSourcesExt\constructors.jl:27
 [17] top-level scope

[rafaqz]

No good reasons, I'll fix that too in the kw PR.