Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Can't correctly recognize netCDF file #737

Open
melodyjulia opened this issue Sep 12, 2024 · 8 comments
Open

Can't correctly recognize netCDF file #737

melodyjulia opened this issue Sep 12, 2024 · 8 comments

Comments

@melodyjulia
Copy link

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}
@rafaqz
Copy link
Owner

rafaqz commented Sep 12, 2024

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

@melodyjulia
Copy link
Author

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

@rafaqz
Copy link
Owner

rafaqz commented Sep 12, 2024

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)

@melodyjulia
Copy link
Author

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.

@rafaqz
Copy link
Owner

rafaqz commented Sep 12, 2024

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

@melodyjulia melodyjulia changed the title Can't correctly recognize netCDF file outputed from python package xESMF Can't correctly recognize netCDF file Sep 12, 2024
@melodyjulia
Copy link
Author

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

@melodyjulia
Copy link
Author

melodyjulia commented Sep 13, 2024

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

@rafaqz
Copy link
Owner

rafaqz commented Sep 13, 2024

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.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
Development

No branches or pull requests

2 participants