Support lon/lat bounds in QuadMesh

[zxdawn]

pcolormesh will work correctly for 2D data with bounds according to the definition.

I tested the Dataset below:

<xarray.Dataset>
Dimensions:                              (y_bounds: 4173, x_bounds: 451, y: 4172, x: 450)
Coordinates:
    lon_bounds                           (y_bounds, x_bounds) float32 dask.array<chunksize=(4096, 451), meta=np.ndarray>
    lat_bounds                           (y_bounds, x_bounds) float32 dask.array<chunksize=(4096, 451), meta=np.ndarray>
    time                                 datetime64[ns] 2019-08-10
    latitude                             (y, x) float32 dask.array<chunksize=(4096, 450), meta=np.ndarray>
    longitude                            (y, x) float32 dask.array<chunksize=(4096, 450), meta=np.ndarray>
    crs                                  object +proj=latlong +datum=WGS84 +e...
Dimensions without coordinates: y_bounds, x_bounds, y, x
Data variables:
    nitrogendioxide_tropospheric_column  (y, x) float32 dask.array<chunksize=(4096, 450), meta=np.ndarray>

import geoviews as gv

gvds = gv.Dataset(ds.unify_chunks()).persist()

mesh = gvds.to(gv.QuadMesh,
        kdims=["lon_bounds", "lat_bounds"],
        vdims=['nitrogendioxide_tropospheric_column'])

Here's the error:

---------------------------------------------------------------------------
KeyError                                  Traceback (most recent call last)
/tmp/ipykernel_163767/1539571452.py in <module>
      2 gvds = gv.Dataset(ds.unify_chunks()).persist()
      3 
----> 4 mesh = gvds.to(gv.QuadMesh,
      5         kdims=["lon_bounds", "lat_bounds"],
      6         vdims=['nitrogendioxide_tropospheric_column'])

~/new/miniconda3/envs/arctic/lib/python3.8/site-packages/geoviews/element/__init__.py in __call__(self, *args, **kwargs)
     37                 args = (Dataset,)
     38                 kwargs['kdims'] = []
---> 39         converted = super(GeoConversion, self).__call__(*args, **kwargs)
     40         if is_gpd:
     41             if kdims is None: kdims = group_type.kdims

~/new/miniconda3/envs/arctic/lib/python3.8/site-packages/holoviews/core/data/__init__.py in __call__(self, new_type, kdims, vdims, groupby, sort, **kwargs)
    142             else:
    143                 selected = self._element.reindex(groupby+kdims, vdims)
--> 144         params = {'kdims': [selected.get_dimension(kd, strict=True) for kd in kdims],
    145                   'vdims': [selected.get_dimension(vd, strict=True) for vd in vdims],
    146                   'label': selected.label}

~/new/miniconda3/envs/arctic/lib/python3.8/site-packages/holoviews/core/data/__init__.py in <listcomp>(.0)
    142             else:
    143                 selected = self._element.reindex(groupby+kdims, vdims)
--> 144         params = {'kdims': [selected.get_dimension(kd, strict=True) for kd in kdims],
    145                   'vdims': [selected.get_dimension(vd, strict=True) for vd in vdims],
    146                   'label': selected.label}

~/new/miniconda3/envs/arctic/lib/python3.8/site-packages/holoviews/core/dimension.py in get_dimension(self, dimension, default, strict)
    987             name_map.update({util.dimension_sanitizer(dim.name): dim for dim in all_dims})
    988             if strict and dimension not in name_map:
--> 989                 raise KeyError("Dimension %r not found." % dimension)
    990             else:
    991                 return name_map.get(dimension, default)

KeyError: "Dimension 'lon_bounds' not found."

So, I suppose we can't use bounds at this stage?

Note that if we don't use bounds for pcolormesh, it will produce strange stripes over the polar region like this:

[jbednar]

I can't quite figure out how this is meant to work; why are lon_bounds and lat_bounds coordinates on their own instead of just being the upper and lower ends of the (apparently missing?) coordinates for the x and y dimensions? I'm very confused!

[zxdawn]

The bounds are actually the upper and lower ends, that's why the length of dimension y_bounds and x_bounds is one more than y and x, respectively.

Here're the values of lon and lon_bounds:

array([[-165.15884 , -164.42227 , -163.68459 , ...,  -45.972137,
         -45.890667,  -45.808647],
       [-164.85976 , -164.12515 , -163.38963 , ...,  -46.105793,
         -46.023884,  -45.941418],
       [-164.56575 , -163.83319 , -163.09987 , ...,  -46.239372,
         -46.157024,  -46.074116],
       ...,
       [  75.04791 ,   74.92377 ,   74.80114 , ...,   33.606133,
          33.449783,   33.291237],
       [  75.002075,   74.87801 ,   74.755455, ...,   33.607353,
          33.451225,   33.292904],
       [  74.956795,   74.83281 ,   74.71033 , ...,   33.6091  ,
          33.4532  ,   33.295113]], dtype=float32)

array([[-165.67815 , -164.94127 , -164.2031  , ...,  -45.86466 ,
         -45.783134,  -45.701057],
       [-165.37622 , -164.6412  , -163.90508 , ...,  -45.9982  ,
         -45.916233,  -45.83371 ],
       [-165.07875 , -164.34566 , -163.61162 , ...,  -46.131596,
         -46.04919 ,  -45.966217],
       ...,
       [  75.087616,   74.962746,   74.8394  , ...,   33.528893,
          33.371563,   33.21201 ],
       [  75.04198 ,   74.91718 ,   74.793915, ...,   33.53045 ,
          33.373344,   33.21402 ],
       [  74.99674 ,   74.872025,   74.748825, ...,   33.53234 ,
          33.37547 ,   33.21638 ]], dtype=float32)

BTW, here're the results using pcolormesh only:

Using center lon/lat

Using lon/lat bounds

[zxdawn]

@jbednar I reproduced this using a simple xarray Dataset and this problem actually exists in applying pcolormesh with projection (see the xarray issue).

It seems we can't get the correct image without the bounds data.

import xarray as xr
import geoviews as gv
import geoviews.feature as gf
from holoviews.operation.datashader import datashade

gv.extension("bokeh")

ds = xr.tutorial.open_dataset('rasm').load()

gvds = gv.Dataset(ds.isel(time=0))

mesh = gvds.to(gv.QuadMesh, kdims=["xc", "yc"], vdims=['Tair'])

datashade(mesh,
          precompute=True,
          cmap='Spectral_r',
          dynamic=True,
         )*\
mesh * gf.coastline.opts(width=500, height=500,
                  tools=['hover'],
                  projection=ccrs.PlateCarree(),
                 )

If we use the polar projection, that's the same problem mentioned before:

[zxdawn]

Note that the newest cartopy (0.20.1) fixes the xarray case, but the result generated by geoviews is still wrong. So, if this problem can be fixed, then supporting bounds doesn't matter.

[zxdawn]

Ha, project before rasterizing like #483 can fix it!

import xarray as xr
import geoviews as gv
import cartopy.crs as ccrs
import geoviews.feature as gf
from holoviews.operation.datashader import datashade

gv.extension("bokeh")

ds = xr.tutorial.open_dataset('rasm').load()

gvds = gv.Dataset(ds.isel(time=0))

mesh = gvds.to(gv.QuadMesh, kdims=["xc", "yc"], vdims=['Tair'])

datashade(gv.project(mesh, projection=ccrs.NorthPolarStereo()),
          precompute=True,
          cmap='Spectral_r',
          dynamic=True,
         )*\
gf.coastline.opts(width=500, height=500,
                  tools=['hover'],
                  projection=ccrs.NorthPolarStereo(),
                 )

[maximlt]

Great thanks @zxdawn for reporting how to do that correctly :)

[zxdawn]

@maximlt You're welcome. I found another problem of setting the extent of the project and posted it on the forum. If that's actually an issue, I can submit it on GitHub. Otherwise, if that's a usage topic, I will keep it there.