Plotting the unstuctured grid data of the WaveWatch III wave model output

Liu-Jincan commented 1 year ago

First of all, thank you so much for developing this tool.

Summary

I have recently been introduced to psyplot and would like to visualize WAVEWATCH III (WW3) wave model unstructured grid data.

Reason

I have tried to plot the WW3 unstructured grid data (hs or fp) using the code from this website, and my code is as follows:

import psyplot.project as psy
import matplotlib as mpl

psy.rcParams['plotter.maps.xgrid'] = False
psy.rcParams['plotter.maps.ygrid'] = False
mpl.rcParams['figure.figsize'] = [10., 8.]

tsunami = psy.plot.mapplot(
    'ww3.20220815T00Z.nc',
    name='hs',
    load=True,
    cmap='Blues')

#
tsunami.show()

The ww3.20220815T00Z.nc could be downloaded via icloud link.

And the running error is as follows:

>>> tsunami = psy.plot.mapplot(
...     '/media/jincanliu/Elements/MaterialLibrary-20220815-20220829-BeijingHaidianCMA524/NC/2022/20220815/ww3.20220815T00Z.nc', name='hs', load=True,
...     maskleq=0, cmap='Blues',
...     clabel='{desc}', stock_img=True, lsm='50m')
Traceback (most recent call last):
  File "<stdin>", line 4, in <module>
  File "/1t/Data-Assimilation-for-Ocean-Current-Forecasts/anaconda3/envs/py-ev-psyplot-2/lib/python3.7/site-packages/psyplot/project.py", line 1987, in __call__
    six.iteritems(self._default_dims), six.iteritems(kwargs))))
  File "/1t/Data-Assimilation-for-Ocean-Current-Forecasts/anaconda3/envs/py-ev-psyplot-2/lib/python3.7/site-packages/psyplot/project.py", line 1772, in _add_data
    return self.project._add_data(*args, **kwargs)
  File "/1t/Data-Assimilation-for-Ocean-Current-Forecasts/anaconda3/envs/py-ev-psyplot-2/lib/python3.7/site-packages/psyplot/project.py", line 180, in wrapper
    return func(self, *args, **kwargs)
  File "/1t/Data-Assimilation-for-Ocean-Current-Forecasts/anaconda3/envs/py-ev-psyplot-2/lib/python3.7/site-packages/psyplot/project.py", line 803, in _add_data
    project=self, enable_post=enable_post, **fmt)
  File "/1t/Data-Assimilation-for-Ocean-Current-Forecasts/anaconda3/envs/py-ev-psyplot-2/lib/python3.7/site-packages/psyplot/plotter.py", line 1261, in __init__
    make_plot=make_plot)
  File "/1t/Data-Assimilation-for-Ocean-Current-Forecasts/anaconda3/envs/py-ev-psyplot-2/lib/python3.7/site-packages/psyplot/plotter.py", line 1414, in initialize_plot
    self.plot_data = self.data
  File "/1t/Data-Assimilation-for-Ocean-Current-Forecasts/anaconda3/envs/py-ev-psyplot-2/lib/python3.7/site-packages/psyplot/plotter.py", line 1162, in plot_data
    self._set_data(value)
  File "/1t/Data-Assimilation-for-Ocean-Current-Forecasts/anaconda3/envs/py-ev-psyplot-2/lib/python3.7/site-packages/psy_simple/plotters.py", line 5891, in _set_data
    raise ValueError(f"Can only plot {self.allowed_dims}-dimensional "
ValueError: Can only plot 2-dimensional data!

I know it's probably because the NetCDF data doesn't meet the UGRID Conventions, but I really don't know how to change it to work.

Output of `psyplot -aV`

```bash # Paste the output of the command ``psyplot -aV`` (ran from the command line) $ psyplot -aV INFO:numexpr.utils:Note: NumExpr detected 16 cores but "NUMEXPR_MAX_THREADS" not set, so enforcing safe limit of 8. INFO:numexpr.utils:NumExpr defaulting to 8 threads. psy_maps.plugin: requirements: cartopy: 0.18.0 version: 1.4.2 psy_reg.plugin: requirements: scipy: 1.7.3 statsmodels: 0.13.2 version: 1.4.0 psy_simple.plugin: version: 1.4.1 psyplot: requirements: matplotlib: 3.5.3 numpy: 1.21.6 pandas: 1.3.5 python: 3.7.12 | packaged by conda-forge | (default, Oct 26 2021, 06:08:53) [GCC 9.4.0] xarray: 0.20.2 version: 1.4.3 psyplot_gui: requirements: pyqt: 5.15.7 qt: 5.15.4 qtconsole: 5.3.2 version: 1.4.0 ```

Chilipp commented 1 year ago

thanks for your issue @Liu-Jincan! There are two challenges here:

missing implementations in psyplot

Your netCDF file does contain variables that are defined on the nodes, rather than on the faces as in the example that you used above. Visualization on the nodes is not yet supported by the current release, but there are ways to circumvent this (https://github.com/psyplot/psyplot/pull/31). Here is how you can create a new conda environment with the necessary packages:

Create a file named environment.yml with the following content:

name: psyplot-ugrid
channels:
 - conda-forge
dependencies:
 - psy-maps
 - pyqt
 - netcdf4
 - dask
 - pip
 - pip:
   - git+https://github.com/Chilipp/psyplot.git@ugrid
   - git+https://github.com/Chilipp/gridded.git@dual_node_mesh

create the conda environment via conda env create -f environment.yml
activate the conda environment via conda activate psyplot-ugrid

Missing mesh information

Your netcdf file is missing the mesh variable that is necessary for the UGRID dimensions. This mesh variable (in this example, the variable Mesh2, defines how to find the face-node-connectivty, what are the parameters, and so on. without this variable, psyplot is not able to visualize the data. However, I was able to figure out the structure of this variable for your file. So here is how you can repair it

import xarray as xr
ds = xr.open_dataset("ww3.20220815T00Z.nc")
ds["mesh"] = xr.Variable((), 1, {"face_node_connectivity": "tri", "node_coordinates": "longitude latitude"})
ds["tri"].attrs["start_index"] = 1
ds["hs"].attrs["mesh"] = "mesh"
ds["fp"].attrs["mesh"] = "mesh"
ds.to_netcdf("ww3.20220815T00Z-fixed.nc")

Visualization

Once you fixed your file and installed the necessary libraries, you should be able to visualize the data, e.g. via

import psyplot.project as psy
sp = psy.plot.mapplot("ww3.20220815T00Z-fixed.nc", name="fp")

or via

psyplot ww3.20220815T00Z-fixed.nc -pm mapplot -n fp -o test.png

which gives you the following output

Liu-Jincan commented 1 year ago

Thank you very much for your reply @Chilipp, and the suggestion is very helpful to me. :smiley:

Here's my plot via psyplot, cool~:

Chilipp commented 1 year ago

Awesome! Glad that I could help :smiley:

psyplot / psyplot