geoviews icon indicating copy to clipboard operation
geoviews copied to clipboard

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

Open hoxbro opened this issue 1 year ago • 2 comments

Fixes https://github.com/holoviz/hvplot/issues/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.

hoxbro avatar Jul 15 '22 13:07 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:

https://user-images.githubusercontent.com/19758978/181490725-fdc7bfc4-d0e4-4636-be19-91db787b42ac.mp4

This is because the srj_proj.threshold is 6680 and therefore ignores the first part of the data when doing the conversion. https://github.com/holoviz/geoviews/blob/a315de762e1a3300877e210aac91798edc3ced29/geoviews/util.py#L88-L90

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.

https://user-images.githubusercontent.com/19758978/181492466-b02d2856-8565-4b01-b35a-0205253a8ebd.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:

https://github.com/holoviz/geoviews/blob/a315de762e1a3300877e210aac91798edc3ced29/geoviews/plotting/plot.py#L72-L73

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:

https://user-images.githubusercontent.com/19758978/181493623-02e1d4e5-005d-4bde-9f8a-7c69b885241c.mp4

hoxbro avatar Jul 28 '22 11:07 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)))

https://user-images.githubusercontent.com/19758978/183426898-17fe7824-074e-40e7-b01e-912add62ddb4.mp4

hoxbro avatar Aug 08 '22 13:08 hoxbro