Better error message and exception for out-of-bound data

[hoxbro]

Fixes holoviz/hvplot#738
Fixes #577

This change will give a better error message and exception for out-of-bound data.

Added from None to avoid During handling of the above exception, another exception occurred message.

[jbednar]

Given that this is a viz tool designed to support data exploration, isn't there an argument for simply warning about data that can't be projected into these coordinates, while projecting the remainder of the data? E.g. if people have data extending to the poles and use a projection not defined for the full extent, seems like they should be able to see the rest of the data anyway. Otherwise they'll have to figure out the extent of this projection and figure out how to clip out that data, right?

[hoxbro]

This is actually what it already does. The IndexError is happening because there is no data after the conversion. I wanted to show a small example of it happening, but I could not get it to work without changing the code. See comment below.

[hoxbro]

This error can occur with the current implementation for coordinates close to the bounds but still inside the bounds.

Code used for the examples

import xarray as xr
import cartopy.crs as ccrs
import panel as pn
import geoviews as gv
import numpy as np

pn.extension()
gv.extension("bokeh", logo=False)

tiles = gv.tile_sources.Wikipedia().opts(width=1000, height=500, padding=0.1)
crs = ccrs.epsg(32635)

def plot(y1):
    da = xr.DataArray(np.arange(10000).reshape(100,100),
                      coords=dict(x=np.linspace(166021, 166021 + 3000, 100),
                                  y=np.linspace(y1, 3000, 100))
                     )
    return gv.Image(da, crs=crs) * tiles


w = pn.widgets.IntSlider(value=0, start=-2000, end=2000, step=100)
pn.Column(w, gv.DynamicMap(pn.bind(plot, w)))



da = xr.DataArray(np.arange(10000).reshape(100,100),
                  coords=dict(x=np.linspace(166021, 166021 + 3000, 100),
                              y=np.linspace(-3000, 0, 100)))
crs = ccrs.epsg(32635)
gv.Image(da, crs=crs) * tiles

An example is if I use CRS-32635. Which haves bounds from (x1, y1, x2, y2) = (166021.44, 0.00, 534994.66, 9329005.18), which means that y1 should always be able to plot if it is above 0. It will, with the current implementation, raise an error even though it shouldn't:

current.mp4

This is because the srj_proj.threshold is 6680 and therefore ignores the first part of the data when doing the conversion.

geoviews/geoviews/util.py

Lines 88 to 90 in a315de7

	# Erode boundary by threshold to avoid transform issues.
	# This is a workaround for numerical issues at the boundary.
	eroded_boundary = boundary_poly.buffer(-src_proj.threshold)

By setting buffer(0), we get the desired result, but when we have data just outside the bounds, we get a thin line. I think this is somewhat related to the padding of the plot, but I'm not entirely sure.

first_attempt.mp4

If I add a keyword argument to project_extents of threshold=0 and then change the code above to eroded_boundary = boundary_poly.buffer(-threshold) and calculate the threshold before passing it to the function, here:

geoviews/geoviews/plotting/plot.py

Lines 72 to 73 in a315de7

	else:
	extents = project_extents(extents, element.crs, proj)

The way I calculate the threshold is very unelegant but is done with:

threshold = min(
    element._SheetCoordinateSystem__xstep,
    element._SheetCoordinateSystem__ystep,
    element.crs.threshold,
)

I get the desired error for the last example:

second_attempt.mp4

[hoxbro]

This will now return this error message:

da = xr.DataArray(np.arange(10000).reshape(100,100),
                  coords=dict(x=t + np.linspace(166021, 166021 + 3000, 100),
                              y=t + np.linspace(-3000, -15, 100)))
crs = ccrs.epsg(32635)
gv.Image(da, crs=crs)# * tiles

ValueError: Could not project data from 'WGS 84 / UTM zone 35N' projection to 'Mercator' projection. 
Ensure some data is inside the bounds +/- the threshold of the CRS. 
For 'WGS 84 / UTM zone 35N' this is (xmin, xmax, ymin, ymax) = (172701, 827299, 6680, 9322326).

---

The examples shown in the previous comment can be archived with the current implementation with:

import xarray as xr
import cartopy.crs as ccrs
import panel as pn
import geoviews as gv
import numpy as np

pn.extension()
gv.extension("bokeh", logo=False)

tiles = gv.tile_sources.Wikipedia().opts(width=1000, height=500, padding=0.1)
crs = ccrs.epsg(32635)
t = int(crs.threshold)


def plot(y1):
    da = xr.DataArray(np.arange(10000).reshape(100,100),
                      coords=dict(x=t + np.linspace(166021, 166021 + 3000, 100),
                                  y=t + np.linspace(y1, 3000, 100))
                     )
    return gv.Image(da, crs=crs) * tiles


w = pn.widgets.IntSlider(value=0, start=-2000, end=2000 , step=100)
pn.Column(w, gv.DynamicMap(pn.bind(plot, w)))

screenrecord-2022-08-08_15.16.25.mp4