Rasters.jl
Rasters.jl copied to clipboard
Can't correctly recognize netCDF file
For the following sample netCDF file described by CDL:
netcdf sample {
dimensions:
lat = 2 ;
lon = 2 ;
variables:
int time ;
time:long_name = "time" ;
time:units = "hours since 1900-01-01" ;
time:calendar = "gregorian" ;
float lat(lat) ;
lat:_FillValue = NaNf ;
lat:units = "degrees_north" ;
float lon(lon) ;
lon:_FillValue = NaNf ;
lon:units = "degress_east" ;
double snowc(lat, lon) ;
snowc:_FillValue = NaN ;
snowc:units = "%" ;
snowc:long_name = "Snow cover" ;
snowc:regrid_method = "conservative" ;
snowc:coordinates = "time" ;
data:
time = 438288 ;
lat = 40.125, 40.375 ;
lon = 90.125, 90.375 ;
snowc =
0.383202466116373, 0.382934974725521,
0.417284701516663, 0.483263764216981 ;
}
It can not be opened:
using Rasters, NCDatasets
Raster("sample.nc","snowc")
It prompts:
ERROR: `Rasters.jl` requires backends to be loaded externally as of version 0.8. Run `import ArchGDAL` to fix this error.
Stacktrace:
[1] _open(f::Function, s::Rasters.GDALsource, filename::String; kw::@Kwargs{name::Rasters.NoKW, group::Rasters.NoKW})
@ Rasters ~/.julia/packages/Rasters/TJQZ7/src/sources/sources.jl:87
[2] Raster(ds::String, filename::String; dims::Rasters.NoKW, refdims::Tuple{}, name::Rasters.NoKW, group::Rasters.NoKW, metadata::Rasters.NoKW, missingval::Rasters.NoKW, crs::Rasters.NoKW, mappedcrs::Rasters.NoKW, source::Rasters.NoKW, replace_missing::Bool, write::Bool, lazy::Bool, dropband::Bool, checkmem::Bool)
@ Rasters ~/.julia/packages/Rasters/TJQZ7/src/array.jl:326
[3] Raster(ds::String, filename::String)
@ Rasters ~/.julia/packages/Rasters/TJQZ7/src/array.jl:308
[4] top-level scope
If specifying data source:
Raster("sample.nc","snowc",source=:netcdf)
It prompts:
ERROR: MethodError: no method matching _fix_missingval(::NCDataset{Nothing, Missing}, ::Rasters.NoKW)
Closest candidates are:
_fix_missingval(::Type, ::Union{Nothing, Rasters.NoKW})
@ Rasters ~/.julia/packages/Rasters/TJQZ7/src/array.jl:347
_fix_missingval(::Type{T}, ::M) where {T, M}
@ Rasters ~/.julia/packages/Rasters/TJQZ7/src/array.jl:351
_fix_missingval(::AbstractArray, ::Rasters.NoKW)
@ Rasters ~/.julia/packages/Rasters/TJQZ7/src/array.jl:349
If using ArchGDAL, it does open the file, but it gives wrong coordinates:
Raster("snowc_xesmf_sample.nc","snowc")
╭───────────────────────╮
│ 2×2 Raster{Float64,2} │
├───────────────────────┴──────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────── dims ┐
↓ X Projected{Float64} LinRange{Float64}(90.0, 90.25, 2) ForwardOrdered Regular Intervals{Start},
→ Y Projected{Float64} LinRange{Float64}(40.25, 40.0, 2) ReverseOrdered Regular Intervals{Start}
Where did you see that syntax? It's not documented like that, it's just an accident that it doesn't error ;)
The syntax is:
Raster(filename; name=:snowc)
I will fix so your syntax throws an error
Thank you for pointing out my error!
The syntax does work if using ArchGDAL:
using Rasters, ArchGDAL
Raster(“sample.nc”, ”snowc”)
But it outputs wrong coordinates:)
I mean did you get the right coordinates with the right syntax?
The problem with your syntax is it only "works" by mistake, it's not doing what you think it is so the values are not what you want.
I need to make it throw an error
(I think it's just using GDAL to load the file, and it's not good at loading netcdf)
Yes, it provides the correct coordinates with the proper syntax. However, with incorrect syntax using ArchGDAL, the coordinates only shift by half a grid width, making it difficult to detect the error. It would be beneficial if throwing an error in such cases.
Yes I will stop that syntax from working at all, its just the fallback function somehow ignoring the second argument.
Thanks. I have updated the title of this issue to avoid misleading other users.
Additionally, there is a minor issue where the variable in the sample NetCDF file is not recognized by Rasters.jl without explicitly specifying the variable name. Here is the code I attempted:
Raster("sample.nc")
It prompts:
╭─────────────────────────────────────────────╮
│ 0-dimensional Raster{Dates.DateTime,0} time │
├─────────────────────────────────────────────┴──────────── metadata ┐
Metadata{Rasters.NCDsource} of Dict{String, Any} with 3 entries:
"units" => "hours since 1900-01-01"
"calendar" => "gregorian"
"long_name" => "time"
├──────────────────────────────────────────────────────────── raster ┤
extent: nothing
└────────────────────────────────────────────────────────────────────┘
1950-01-01T00:00:00
Right, it's ignoring the dimension variables lat/lon but not the coordinates variable. So it takes the first of time and snowc. Will fix.