Rasters.jl icon indicating copy to clipboard operation
Rasters.jl copied to clipboard

Can't correctly recognize netCDF file

Open melodyjulia opened this issue 1 year ago • 8 comments
trafficstars

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}

melodyjulia avatar Sep 12 '24 07:09 melodyjulia

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

rafaqz avatar Sep 12 '24 07:09 rafaqz

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:)

melodyjulia avatar Sep 12 '24 07:09 melodyjulia

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)

rafaqz avatar Sep 12 '24 07:09 rafaqz

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.

melodyjulia avatar Sep 12 '24 08:09 melodyjulia

Yes I will stop that syntax from working at all, its just the fallback function somehow ignoring the second argument.

rafaqz avatar Sep 12 '24 08:09 rafaqz

Thanks. I have updated the title of this issue to avoid misleading other users.

melodyjulia avatar Sep 12 '24 08:09 melodyjulia

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

melodyjulia avatar Sep 13 '24 01:09 melodyjulia

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.

rafaqz avatar Sep 13 '24 06:09 rafaqz