pythonmatplotlibrasterholoviewsgeoviews

From Matplotlib Raster to Geoviews/ Holoviews / hvplot: How to transform x, y and z

I understand that Geoviews and Holoviews share common attributes, and Hvplot is meant to be a high level plotting API to all three.

Now, coming from Matplotlib, I have still difficulties adapting to the parameters required to display raster images in Geoviews or Holoviews.

Here's an example, where I am doing a Kernel Density Estimation for spatial data:

# add coordinates of observations
xy_train = np.vstack([y, x]).T
print(xy_train)
# [[5654810.66920637  413645.79802685]
# [5654712.51814666  412629.87266155]
# [5656120.03682466  411642.74943511]
# ...
# [5656316.96943554  411795.80163676]
# [5656299.73356505  411795.50717494]
# [5655756.85624901  411734.34680852]]

# create mesh
xbins=100j
ybins=100j
xx, yy = np.mgrid[left_bound:right_bound:xbins, 
                bottom_bound:top_bound:ybins]
xy_sample = np.vstack([yy.ravel(), xx.ravel()]).T
# compute Kernel Density here
# ..
kde = KernelDensity(kernel='gaussian', bandwidth=100, algorithm='ball_tree')
kde.fit(xy_train)
# get results
z = np.exp(kde.score_samples(xy_sample))
# reshape results to mesh
zz = z.reshape(xx.shape)
# plot in matplotlib
fig, ax_lst = plt.subplots(111)
levels = np.linspace(zz.min(), zz.max(), 25)
axis.contourf(xx, yy, zz, levels=levels, cmap='viridis')
axis.plot()
plt.show()

Shows my Image: enter image description here

Now I want to use the pyviz environment for interactive display and map overlay, e.g. using Geoviews.

This somehow works, but gives me an error:

xr_dataset = gv.Dataset(hv.Image((xx, yy, zz.T), datatype=['grid']), crs=ccrs.UTM(zone='33N'))

Image02195: Image dimension(s) x and y are not evenly sampled to relative tolerance of 0.001. Please use the QuadMesh element for irregularly sampled data or set a higher tolerance on hv.config.image_rtol or the rtol parameter in the Image constructor.

I can still display the image (somehow low resolution).

gv.tile_sources.CartoLight.opts(width=800, height=480) * xr_dataset.to.image(['x', 'y']).opts(cmap='viridis', alpha=0.5)

enter image description here

.. but when I try to create FilledContours in Geoviews, it doesn't seem to work like in matplotlib:

gv.FilledContours(xx, yy, zz, levels=levels, cmap='viridis')

ValueError: kdims argument expects a Dimension or list of dimensions, specified as tuples, strings, dictionaries or Dimension instances, not a ndarray type. Ensure you passed the data as the first argument.

The documentation doesn't provide much info on how I should format dimensions (hv.help(gv.FilledContours)). I think I get lost somewhere when I need to create a Raster from the numpy xx/yy coordinate mesh (hv.Image((xx, yy, zz.T), datatype=['grid'])).

Can someone explain the difference in syntax that is required for matplotlib Contourf and Holoviews/Geoviews/Hvplot FilledContours?

[edit]

I found a way to create contours, but the dimensions problem persists:

# get xarray dataset, suited for handling raster data in pyviz
xr_dataset = gv.Dataset(hv.Image((xx.T, yy.T, zz.T), bounds=(left_bound,bottom_bound,right_bound,top_bound), 
        kdims=[hv.Dimension('x'),  hv.Dimension('y')], datatype=['grid']), crs=ccrs.UTM(zone='33N'))
# Error: Image06593: Image dimension(s) x and y are not evenly sampled to relative tolerance of 0.001

# create contours from image
gv.FilledContours(xr_dataset)

# plot
gv.tile_sources.EsriImagery.opts(width=800, height=480) * gv.FilledContours(xr_dataset).opts(cmap='viridis', alpha=0.5)

enter image description here

Solution

  • The main thing to know about HoloViews/GeoViews elements is that the data is almost always specified as the first argument, which is unlike matplotlib where the data is often specified using multiple arguments. In your case you already had the correct syntax for an Image but didn't carry that over to other elements. So to make this concrete, to construct an Image you would do:

    img = gv.Image((xx, yy, zz.T), crs=ccrs.UTM(zone='33N'))

    However since you have 2D coordinate arrays rather than the 1D coordinates that an Image expects (in the next release this will error), you actually have a QuadMesh, which gets constructed in the same way:

    qmesh = gv.QuadMesh((xx, yy, zz.T), crs=ccrs.UTM(zone='33N'))

    And the same is also true for the geoviews FilledContours:

    contours = gv.FilledContours((xx, yy, zz.T), crs=ccrs.UTM(zone='33N'))

    So to summarize, the difference between HoloViews elements and matplotlib calls is that HoloViews is a lightweight wrapper around your data, which lets you give each coordinate and value array semantic meaning by assigning them to a key or value dimension, while matplotlib makes that mapping more explicit.

    HoloViews understands a number of formats for defining gridded data like yours, the simplest is a tuple of x/y-coordinate arrays and the value array, it also understands xarray objects and dictionaries of the different arrays which would look like this:

    contours = gv.FilledContours({'x': xx, 'y': yy, 'z': zz.T}, ['x', 'y'], 'z', crs=ccrs.UTM(zone='33N'))

    In this format we can explicitly see how the 'x' and 'y' coordinate arrays are mapped to the key dimensions and the 'z' value array to the value dimensions.