Visualising NetCDF using Raster.jl

Michaela_Vollmer · August 15, 2022, 6:57am

Hello! I am trying to plot spatial data of a NetCDF file from a WFlow example. I have followed the steps outlined in the following question Raster.jl discourse post where a similar example is shown where the data can be plotted. However, if I apply the same principle to the WFLow example I get the following error:

ERROR: MethodError: no method matching reverse(::DimensionalData.Dimensions.LookupArrays.AutoOrder)
Closest candidates are:

The following is a minimal working example:

using NCDatasets
using Rasters, GBIF, Plots

fn = download("https://github.com/visr/wflow-artifacts/releases/download/v0.2.6/forcing-moselle.nc",joinpath("./", "forcing-moselle.nc"))
fn="forcing-moselle.nc"
fn=Raster(fn)
wkt = WellKnownText(Dataset(fn)["spatial_ref"].attrib["crs_wkt"])

# Choose the prec layer from the netcdf file
prec = Raster(fn; name=:precip, crs=wkt, mappedcrs=wkt)

# Use typemin(Int32) as the missing value.
prec = replace_missing(prec, typemin(Int32))

#resampled = resample(prec, 1; crs=EPSG(4326))

plot(prec)

How can this be fixed? Thank you very much in advance!

digital_carver · August 15, 2022, 8:13am

Please post the full stacktrace of the error, since this doesn’t tell us which line the error is coming from.

visr · August 15, 2022, 8:42am

Here is the full stacktrace. Tagging @Raf since this may be an issue in DimensionalData.jl.

ERROR: MethodError: no method matching reverse(::DimensionalData.Dimensions.LookupArrays.AutoOrder)
Closest candidates are:
  reverse(::Union{SubString{String}, String}) at strings/substring.jl:171
  reverse(::DimensionalData.Dimensions.LookupArrays.AbstractCategorical) at d:\visser_mn\.julia\packages\DimensionalData\1iOoH\src\LookupArrays\methods.jl:5
  reverse(::DimensionalData.Dimensions.LookupArrays.AbstractSampled) at d:\visser_mn\.julia\packages\DimensionalData\1iOoH\src\LookupArrays\methods.jl:10
  ...
Stacktrace:
  [1] reverse(lookup::DimensionalData.Dimensions.LookupArrays.Categorical{Union{Missing, Float64}, Vector{Union{Missing, Float64}}, DimensionalData.Dimensions.LookupArrays.AutoOrder, DimensionalData.Dimensions.LookupArrays.Metadata{Rasters.NCDfile, Dict{Symbol, Any}}})
    @ DimensionalData.Dimensions.LookupArrays d:\visser_mn\.julia\packages\DimensionalData\1iOoH\src\LookupArrays\methods.jl:7
  [2] reverse(dim::X{DimensionalData.Dimensions.LookupArrays.Categorical{Union{Missing, Float64}, Vector{Union{Missing, Float64}}, DimensionalData.Dimensions.LookupArrays.AutoOrder, DimensionalData.Dimensions.LookupArrays.Metadata{Rasters.NCDfile, Dict{Symbol, Any}}}})
    @ DimensionalData d:\visser_mn\.julia\packages\DimensionalData\1iOoH\src\array\methods.jl:307
  [3] _reverse(dim::X{DimensionalData.Dimensions.LookupArrays.Categorical{Union{Missing, Float64}, Vector{Union{Missing, Float64}}, DimensionalData.Dimensions.LookupArrays.AutoOrder, DimensionalData.Dimensions.LookupArrays.Metadata{Rasters.NCDfile, Dict{Symbol, Any}}}})
    @ DimensionalData d:\visser_mn\.julia\packages\DimensionalData\1iOoH\src\array\methods.jl:304
  [4] reverse(A::Raster{Float32, 3, Tuple{X{DimensionalData.Dimensions.LookupArrays.Categorical{Union{Missing, Float64}, Vector{Union{Missing, Float64}}, DimensionalData.Dimensions.LookupArrays.AutoOrder, DimensionalData.Dimensions.LookupArrays.Metadata{Rasters.NCDfile, Dict{Symbol, Any}}}}, Y{DimensionalData.Dimensions.LookupArrays.Categorical{Union{Missing, Float64}, Vector{Union{Missing, Float64}}, DimensionalData.Dimensions.LookupArrays.AutoOrder, DimensionalData.Dimensions.LookupArrays.Metadata{Rasters.NCDfile, Dict{Symbol, Any}}}}, Ti{DimensionalData.Dimensions.LookupArrays.Sampled{Dates.DateTime, Vector{Dates.DateTime}, DimensionalData.Dimensions.LookupArrays.ForwardOrdered, DimensionalData.Dimensions.LookupArrays.Irregular{Tuple{Nothing, Nothing}}, DimensionalData.Dimensions.LookupArrays.Points, DimensionalData.Dimensions.LookupArrays.Metadata{Rasters.NCDfile, Dict{Symbol, Any}}}}}, Tuple{}, Array{Float32, 3}, Symbol, DimensionalData.Dimensions.LookupArrays.Metadata{Rasters.NCDfile, Dict{Symbol, Any}}, Float32}; dims::X{DimensionalData.Dimensions.LookupArrays.Categorical{Union{Missing, Float64}, Vector{Union{Missing, Float64}}, DimensionalData.Dimensions.LookupArrays.AutoOrder, DimensionalData.Dimensions.LookupArrays.Metadata{Rasters.NCDfile, Dict{Symbol, Any}}}})
    @ DimensionalData d:\visser_mn\.julia\packages\DimensionalData\1iOoH\src\array\methods.jl:297
  [5] _reorder(#unused#::Type{DimensionalData.Dimensions.LookupArrays.ForwardOrdered}, x::Raster{Float32, 3, Tuple{X{DimensionalData.Dimensions.LookupArrays.Categorical{Union{Missing, Float64}, Vector{Union{Missing, Float64}}, DimensionalData.Dimensions.LookupArrays.AutoOrder, DimensionalData.Dimensions.LookupArrays.Metadata{Rasters.NCDfile, Dict{Symbol, Any}}}}, Y{DimensionalData.Dimensions.LookupArrays.Categorical{Union{Missing, Float64}, Vector{Union{Missing, Float64}}, DimensionalData.Dimensions.LookupArrays.AutoOrder, DimensionalData.Dimensions.LookupArrays.Metadata{Rasters.NCDfile, Dict{Symbol, Any}}}}, Ti{DimensionalData.Dimensions.LookupArrays.Sampled{Dates.DateTime, Vector{Dates.DateTime}, DimensionalData.Dimensions.LookupArrays.ForwardOrdered, DimensionalData.Dimensions.LookupArrays.Irregular{Tuple{Nothing, Nothing}}, DimensionalData.Dimensions.LookupArrays.Points, DimensionalData.Dimensions.LookupArrays.Metadata{Rasters.NCDfile, Dict{Symbol, Any}}}}}, Tuple{}, Array{Float32, 3}, Symbol, DimensionalData.Dimensions.LookupArrays.Metadata{Rasters.NCDfile, Dict{Symbol, Any}}, Float32}, dim::X{DimensionalData.Dimensions.LookupArrays.Categorical{Union{Missing, Float64}, Vector{Union{Missing, Float64}}, DimensionalData.Dimensions.LookupArrays.AutoOrder, DimensionalData.Dimensions.LookupArrays.Metadata{Rasters.NCDfile, Dict{Symbol, Any}}}})
    @ DimensionalData d:\visser_mn\.julia\packages\DimensionalData\1iOoH\src\utils.jl:33
  [6] reorder(x::Raster{Float32, 3, Tuple{X{DimensionalData.Dimensions.LookupArrays.Categorical{Union{Missing, Float64}, Vector{Union{Missing, Float64}}, DimensionalData.Dimensions.LookupArrays.AutoOrder, DimensionalData.Dimensions.LookupArrays.Metadata{Rasters.NCDfile, Dict{Symbol, Any}}}}, Y{DimensionalData.Dimensions.LookupArrays.Categorical{Union{Missing, Float64}, Vector{Union{Missing, Float64}}, DimensionalData.Dimensions.LookupArrays.AutoOrder, DimensionalData.Dimensions.LookupArrays.Metadata{Rasters.NCDfile, Dict{Symbol, Any}}}}, Ti{DimensionalData.Dimensions.LookupArrays.Sampled{Dates.DateTime, Vector{Dates.DateTime}, DimensionalData.Dimensions.LookupArrays.ForwardOrdered, DimensionalData.Dimensions.LookupArrays.Irregular{Tuple{Nothing, Nothing}}, DimensionalData.Dimensions.LookupArrays.Points, DimensionalData.Dimensions.LookupArrays.Metadata{Rasters.NCDfile, Dict{Symbol, Any}}}}}, Tuple{}, Array{Float32, 3}, Symbol, DimensionalData.Dimensions.LookupArrays.Metadata{Rasters.NCDfile, Dict{Symbol, Any}}, Float32}, orderdim::X{DataType})
    @ DimensionalData d:\visser_mn\.julia\packages\DimensionalData\1iOoH\src\utils.jl:29
  [7] _reorder
    @ d:\visser_mn\.julia\packages\DimensionalData\1iOoH\src\utils.jl:26 [inlined]
  [8] reorder
    @ d:\visser_mn\.julia\packages\DimensionalData\1iOoH\src\utils.jl:21 [inlined]
  [9] _prepare(A::Raster{Float32, 3, Tuple{X{DimensionalData.Dimensions.LookupArrays.Categorical{Union{Missing, Float64}, Vector{Union{Missing, Float64}}, DimensionalData.Dimensions.LookupArrays.AutoOrder, DimensionalData.Dimensions.LookupArrays.Metadata{Rasters.NCDfile, Dict{Symbol, Any}}}}, Y{DimensionalData.Dimensions.LookupArrays.Categorical{Union{Missing, Float64}, Vector{Union{Missing, Float64}}, DimensionalData.Dimensions.LookupArrays.AutoOrder, DimensionalData.Dimensions.LookupArrays.Metadata{Rasters.NCDfile, Dict{Symbol, Any}}}}, Ti{DimensionalData.Dimensions.LookupArrays.Sampled{Dates.DateTime, Vector{Dates.DateTime}, DimensionalData.Dimensions.LookupArrays.ForwardOrdered, DimensionalData.Dimensions.LookupArrays.Irregular{Tuple{Nothing, Nothing}}, DimensionalData.Dimensions.LookupArrays.Points, DimensionalData.Dimensions.LookupArrays.Metadata{Rasters.NCDfile, Dict{Symbol, Any}}}}}, Tuple{}, Array{Float32, 3}, Symbol, DimensionalData.Dimensions.LookupArrays.Metadata{Rasters.NCDfile, Dict{Symbol, Any}}, Float32})
    @ Rasters d:\visser_mn\.julia\packages\Rasters\UC1tU\src\plotrecipes.jl:195
 [10] macro expansion
    @ d:\visser_mn\.julia\packages\Rasters\UC1tU\src\plotrecipes.jl:19 [inlined]
 [11] apply_recipe(plotattributes::AbstractDict{Symbol, Any}, A::AbstractRaster)
    @ Rasters d:\visser_mn\.julia\packages\RecipesBase\qpxEX\src\RecipesBase.jl:289
 [12] _process_userrecipes!(plt::Any, plotattributes::Any, args::Any)
    @ RecipesPipeline d:\visser_mn\.julia\packages\RecipesPipeline\OXGmH\src\user_recipe.jl:36
 [13] recipe_pipeline!(plt::Any, plotattributes::Any, args::Any)
    @ RecipesPipeline d:\visser_mn\.julia\packages\RecipesPipeline\OXGmH\src\RecipesPipeline.jl:70
 [14] _plot!(plt::Plots.Plot, plotattributes::Any, args::Any)
    @ Plots d:\visser_mn\.julia\packages\Plots\lW9ll\src\plot.jl:209
 [15] plot(args::Any; kw::Base.Pairs{Symbol, V, Tuple{Vararg{Symbol, N}}, NamedTuple{names, T}} where {V, N, names, T<:Tuple{Vararg{Any, N}}})
    @ Plots d:\visser_mn\.julia\packages\Plots\lW9ll\src\plot.jl:91
 [16] plot(args::Any)
    @ Plots d:\visser_mn\.julia\packages\Plots\lW9ll\src\plot.jl:82

joa-quim · August 15, 2022, 10:08am

You can show it with GMT. But the data has outliers so it need a manually tuned color map

using GMT

C = makecpt(0,5);
imshow("forcing-moselle.nc?precip", proj=:guess, cmap=C, colorbar=true)

But note that GMT automatically installation on Linux is broken due to dependencies conflicts external to the package itself.

Raf · August 15, 2022, 10:28am

Sorry, looks like I’ve introduced a bug with some recent changes.

Should be fixed in the next 24 hours.

Michaela_Vollmer · August 15, 2022, 11:02am

Thanks for this. I will give this a try

Michaela_Vollmer · August 15, 2022, 11:03am

Thanks! That would be great!

Raf · August 15, 2022, 5:06pm

Ok this is partly a bug in NCDatasets.jl - the dimension lookups X and Y are being given the type Vector{Union{Mising,Float64}} - which means they miss dispatch in Rasters.jl and are getting assigned to a fallback and partly broken.

That will be fixed, because they don’t actually contain any missing values. But they shouldn’t have Missing in the type - actual missing values in the dimension variable would break most of the use cases of the dimensions.

~~(I’ll make an issue on NCDatasets.jl github for this)~~

Raf · August 15, 2022, 5:15pm

It seems this is actually a problem with your data - your X and Y coordinate variables have _FillValue = NaN.

This is not allowed in CF standards, (see NetCDF Climate and Forecast (CF) Metadata Conventions and search for _FillValue if you are interested)

We can probably handle it more gracefully, but you should also probably speak to whoever is making this dataset and remove the _FillValue attribute from coordinate variables. (Haha @visr I guess thats you, I just realised where it comes from)

Edit:

One last thing - now this is fixed and plots locally (will push a new version in the next few hours), there are 366 time slices! Unlike imshow in GMT, plot(prec) with Rasters.jl will try to tile them all. Which will take forever and wont be what you want.

Instead, do plot(prec[Ti=1]) or plot(prec[Ti=At(DateTime(2000, 1, 2))]) to plot one of the slices, or use a Selector to choose slices you want to plot, like Where to plot the first day in some specific months:

using Dates
plot(prec[Ti=Where(d -> dayofmonth(d) == 1 && month(d) in (Jan, Apr, Jul, Oct))])

visr · August 16, 2022, 12:18pm

Haha nice job @Raf you managed to trace the issue back to me. I can imagine actual missing values on coordinate are not allowed, but didn’t know having a _FillValue is explicitly disallowed under the CF conventions. I added them because xarray adds them. So I removed them now in remove _FillValue from coordinate variables by visr · Pull Request #215 · Deltares/Wflow.jl · GitHub. The wflow-artifacts datasets linked above have now been fixed by ncatted -O -a _FillValue,lat,d,, inmaps.nc for x/y/lat/lon. There is an upstream xarray issue as well: to_netcdf -> _fill_value without NaN · Issue #2037 · pydata/xarray · GitHub. I’m afraid having fill values on coordinate variables is not uncommon.

Raf · August 16, 2022, 9:29pm

Sweet. I also made an issue at NCDatasets.jl to track the fact that the FillValue field propagates to the dimension variable eltype: https://github.com/Alexander-Barth/NCDatasets.jl/issues/185.

But I also fixed Rasters.jl to handle that when it happens, which should be registered shortly.

Of course some things will break if there are actually any missing values… findfirst and searchsorted will not like that at all. It’s possible Missing in the eltype will break some base functions as well, I’m not sure though.