search

ecmwf / earthkit-plots

Visualisation tools and templates designed for earth science data.
https://earthkit-plots.readthedocs.io/en/latest/
7 stars 1 forks source link

ValueError with GRIB data #34

Closed SF-N closed 6 days ago

SF-N commented 6 days ago

What happened?

Hi, after updating earthkit-plots, I get the following error when trying to plot GRIB data, which is available here:

import xarray as xr
import earthkit.plots
import earthkit.plots.quickmap
import numpy as np

ds = xr.open_dataset("hplp/hplp_ml_params_133_203_level_1-2.grib", backend_kwargs=dict(indexpath=""))

print(ds["q"].sel(hybrid=1).sel(step=np.timedelta64(0)))

earthkit.plots.quickmap.plot(
    ds["q"].sel(hybrid=1).sel(step=np.timedelta64(0)),
    x="longitude", y="latitude",
)
<xarray.DataArray 'q' (values: 654400)>
[654400 values with dtype=float32]
Coordinates:
    time        datetime64[ns] ...
    step        timedelta64[ns] 00:00:00
    hybrid      float64 1.0
    latitude    (values) float64 ...
    longitude   (values) float64 ...
    valid_time  datetime64[ns] ...
Dimensions without coordinates: values
Attributes: (12/20)
    GRIB_paramId:                    133
    GRIB_dataType:                   fc
    GRIB_numberOfPoints:             654400
    GRIB_typeOfLevel:                hybrid
    GRIB_stepUnits:                  1
    GRIB_stepType:                   instant
    ...                              ...
    GRIB_pl:                         [  20   24   28   32   36   40   44   48...
    GRIB_shortName:                  q
    GRIB_units:                      kg kg**-1
    long_name:                       Specific humidity
    units:                           kg kg**-1
    standard_name:                   specific_humidity

---------------------------------------------------------------------------
ValueError                                Traceback (most recent call last)
Cell In[12], line 6
      3 import numpy as np
      5 print(ds["q"].sel(hybrid=1).sel(step=np.timedelta64(0)))
----> 6 earthkit.plots.quickmap.plot(
      7     ds["q"].sel(hybrid=1).sel(step=np.timedelta64(0)),
      8     x="longitude", y="latitude",
      9 )

File [~/Documents/ECMWF/Code/data_compression/venv/lib/python3.10/site-packages/earthkit/plots/quickmap.py:9](http://localhost:8888/lab/tree/compression-lab-notebooks/~/Documents/ECMWF/Code/data_compression/venv/lib/python3.10/site-packages/earthkit/plots/quickmap.py#line=8), in _quickmap.<locals>.wrapper(return_subplot, domain, *args, **kwargs)
      7 figure = Figure()
      8 subplot = figure.add_map(domain=domain)
----> 9 getattr(subplot, function.__name__)(*args, **kwargs)
     10 for method in schema.quickmap_workflow:
     11     getattr(subplot, method)()

File [~/Documents/ECMWF/Code/data_compression/venv/lib/python3.10/site-packages/earthkit/plots/components/subplots.py:597](http://localhost:8888/lab/tree/compression-lab-notebooks/~/Documents/ECMWF/Code/data_compression/venv/lib/python3.10/site-packages/earthkit/plots/components/subplots.py#line=596), in Subplot.plot(self, style, *args, **kwargs)
    595 else:
    596     method = self.block
--> 597 return method(*args, style=style, **kwargs)

File [~/Documents/ECMWF/Code/data_compression/venv/lib/python3.10/site-packages/earthkit/plots/components/subplots.py:227](http://localhost:8888/lab/tree/compression-lab-notebooks/~/Documents/ECMWF/Code/data_compression/venv/lib/python3.10/site-packages/earthkit/plots/components/subplots.py#line=226), in Subplot.plot_3D.<locals>.decorator.<locals>.wrapper(self, data, x, y, z, style, every, **kwargs)
    217 def wrapper(
    218     self,
    219     data=None,
   (...)
    225     **kwargs,
    226 ):
--> 227     return self._extract_plottables(
    228         method_name or method.__name__,
    229         args=tuple(),
    230         data=data,
    231         x=x,
    232         y=y,
    233         z=z,
    234         style=style,
    235         every=every,
    236         extract_domain=extract_domain,
    237         **kwargs,
    238     )

File [~/Documents/ECMWF/Code/data_compression/venv/lib/python3.10/site-packages/earthkit/plots/components/subplots.py:396](http://localhost:8888/lab/tree/compression-lab-notebooks/~/Documents/ECMWF/Code/data_compression/venv/lib/python3.10/site-packages/earthkit/plots/components/subplots.py#line=395), in Subplot._extract_plottables(self, method_name, args, data, x, y, z, style, units, every, source_units, extract_domain, **kwargs)
    392     warnings.warn(
    393         "The 'interpolation_method' argument is only valid for unstructured data."
    394     )
    395 try:
--> 396     mappable = getattr(style, method_name)(
    397         self.ax, x_values, y_values, z_values, **kwargs
    398     )
    399 except TypeError as err:
    400     if not grids.is_structured(x_values, y_values):

File [~/Documents/ECMWF/Code/data_compression/venv/lib/python3.10/site-packages/earthkit/plots/styles/__init__.py:429](http://localhost:8888/lab/tree/compression-lab-notebooks/~/Documents/ECMWF/Code/data_compression/venv/lib/python3.10/site-packages/earthkit/plots/styles/__init__.py#line=428), in Style.pcolormesh(self, ax, x, y, values, *args, **kwargs)
    427 kwargs.pop("transform_first", None)
    428 kwargs = {**self.to_pcolormesh_kwargs(values), **kwargs}
--> 429 result = ax.pcolormesh(x, y, values, *args, **kwargs)
    430 return result

File [~/Documents/ECMWF/Code/data_compression/venv/lib/python3.10/site-packages/cartopy/mpl/geoaxes.py:307](http://localhost:8888/lab/tree/compression-lab-notebooks/~/Documents/ECMWF/Code/data_compression/venv/lib/python3.10/site-packages/cartopy/mpl/geoaxes.py#line=306), in _add_transform.<locals>.wrapper(self, *args, **kwargs)
    302     raise ValueError(f'Invalid transform: Spherical {func.__name__} '
    303                      'is not supported - consider using '
    304                      'PlateCarree[/RotatedPole.](http://localhost:8888/RotatedPole.)')
    306 kwargs['transform'] = transform
--> 307 return func(self, *args, **kwargs)

File [~/Documents/ECMWF/Code/data_compression/venv/lib/python3.10/site-packages/cartopy/mpl/geoaxes.py:1767](http://localhost:8888/lab/tree/compression-lab-notebooks/~/Documents/ECMWF/Code/data_compression/venv/lib/python3.10/site-packages/cartopy/mpl/geoaxes.py#line=1766), in GeoAxes.pcolormesh(self, *args, **kwargs)
   1756 """
   1757 Add the "transform" keyword to :func:`~matplotlib.pyplot.pcolormesh`.
   1758 
   (...)
   1763 
   1764 """
   1765 # Add in an argument checker to handle Matplotlib's potential
   1766 # interpolation when coordinate wraps are involved
-> 1767 args, kwargs = self._wrap_args(*args, **kwargs)
   1768 result = super().pcolormesh(*args, **kwargs)
   1769 # Wrap the quadrilaterals if necessary

File [~/Documents/ECMWF/Code/data_compression/venv/lib/python3.10/site-packages/cartopy/mpl/geoaxes.py:1794](http://localhost:8888/lab/tree/compression-lab-notebooks/~/Documents/ECMWF/Code/data_compression/venv/lib/python3.10/site-packages/cartopy/mpl/geoaxes.py#line=1793), in GeoAxes._wrap_args(self, *args, **kwargs)
   1792 X = np.asanyarray(args[0])
   1793 Y = np.asanyarray(args[1])
-> 1794 nrows, ncols = np.asanyarray(args[2]).shape[:2]
   1795 Nx = X.shape[-1]
   1796 Ny = Y.shape[0]

ValueError: not enough values to unpack (expected 2, got 1)

The code used to work fine with a previous earthkit version (unfortunately I don't recall which one anymore).

Currently installed package versions: Cartopy 0.23.0 earthkit 0.8.4 earthkit-data 0.11.1 earthkit-geo 0.3.0 earthkit-maps 0.0.19 earthkit-meteo 0.3.0 earthkit-plots 0.2.7 earthkit-plots-default-styles 0.0.1 earthkit-regrid 0.3.4 earthkit-transforms 0.3.4.1

Could you please let me know which package versions to use and how to plot this data?

What are the steps to reproduce the bug?

Please see the code above.

Version

0.2.7

Platform (OS and architecture)

macOS 14.4.1

Relevant log output

No response

Accompanying data

No response

Organisation

No response

juntyr commented 6 days ago

This should be fixed by #31

SF-N commented 6 days ago

Thanks, this fixes it.