IndexError when using geographic data

[hxawax]

Hey,

I'm reporting a bug occurring when using hvplot with geo=True option or crs= option

ALL software version info

hvplot version : 0.7.3
xarray version 2022.3.0
geoviews version 1.9.5

Description of expected behavior and the observed behavior

without option the dataset is showing fine xr.open_dataset('bug.nc').band_data.hvplot(x='x',y='y')

Complete, minimal, self-contained example code that reproduces the issue

If I add the geo=Trueoption or crs= option or both the following error occur
xr.open_dataset('bug.nc').band_data.hvplot(crs=cartopy.crs.epsg(2154))

---------------------------------------------------------------------------
IndexError                                Traceback (most recent call last)
File ~/miniconda3/envs/baleine/lib/python3.9/site-packages/IPython/core/formatters.py:973, in MimeBundleFormatter.__call__(self, obj, include, exclude)
    970     method = get_real_method(obj, self.print_method)
    972     if method is not None:
--> 973         return method(include=include, exclude=exclude)
    974     return None
    975 else:

File ~/miniconda3/envs/baleine/lib/python3.9/site-packages/holoviews/core/dimension.py:1316, in Dimensioned._repr_mimebundle_(self, include, exclude)
   1309 def _repr_mimebundle_(self, include=None, exclude=None):
   1310     """
   1311     Resolves the class hierarchy for the class rendering the
   1312     object using any display hooks registered on Store.display
   1313     hooks.  The output of all registered display_hooks is then
   1314     combined and returned.
   1315     """
-> 1316     return Store.render(self)

File ~/miniconda3/envs/baleine/lib/python3.9/site-packages/holoviews/core/options.py:1405, in Store.render(cls, obj)
   1403 data, metadata = {}, {}
   1404 for hook in hooks:
-> 1405     ret = hook(obj)
   1406     if ret is None:
   1407         continue

File ~/miniconda3/envs/baleine/lib/python3.9/site-packages/holoviews/ipython/display_hooks.py:282, in pprint_display(obj)
    280 if not ip.display_formatter.formatters['text/plain'].pprint:
    281     return None
--> 282 return display(obj, raw_output=True)

File ~/miniconda3/envs/baleine/lib/python3.9/site-packages/holoviews/ipython/display_hooks.py:252, in display(obj, raw_output, **kwargs)
    250 elif isinstance(obj, (CompositeOverlay, ViewableElement)):
    251     with option_state(obj):
--> 252         output = element_display(obj)
    253 elif isinstance(obj, (Layout, NdLayout, AdjointLayout)):
    254     with option_state(obj):

File ~/miniconda3/envs/baleine/lib/python3.9/site-packages/holoviews/ipython/display_hooks.py:146, in display_hook.<locals>.wrapped(element)
    144 try:
    145     max_frames = OutputSettings.options['max_frames']
--> 146     mimebundle = fn(element, max_frames=max_frames)
    147     if mimebundle is None:
    148         return {}, {}

File ~/miniconda3/envs/baleine/lib/python3.9/site-packages/holoviews/ipython/display_hooks.py:192, in element_display(element, max_frames)
    189 if type(element) not in Store.registry[backend]:
    190     return None
--> 192 return render(element)

File ~/miniconda3/envs/baleine/lib/python3.9/site-packages/holoviews/ipython/display_hooks.py:68, in render(obj, **kwargs)
     65 if renderer.fig == 'pdf':
     66     renderer = renderer.instance(fig='png')
---> 68 return renderer.components(obj, **kwargs)

File ~/miniconda3/envs/baleine/lib/python3.9/site-packages/holoviews/plotting/renderer.py:410, in Renderer.components(self, obj, fmt, comm, **kwargs)
    408 doc = Document()
    409 with config.set(embed=embed):
--> 410     model = plot.layout._render_model(doc, comm)
    411 if embed:
    412     return render_model(model, comm)

File ~/miniconda3/envs/baleine/lib/python3.9/site-packages/panel/viewable.py:455, in Renderable._render_model(self, doc, comm)
    453 if comm is None:
    454     comm = state._comm_manager.get_server_comm()
--> 455 model = self.get_root(doc, comm)
    457 if config.embed:
    458     embed_state(self, model, doc,
    459                 json=config.embed_json,
    460                 json_prefix=config.embed_json_prefix,
    461                 save_path=config.embed_save_path,
    462                 load_path=config.embed_load_path,
    463                 progress=False)

File ~/miniconda3/envs/baleine/lib/python3.9/site-packages/panel/viewable.py:512, in Renderable.get_root(self, doc, comm, preprocess)
    495 """
    496 Returns the root model and applies pre-processing hooks
    497 
   (...)
    509 Returns the bokeh model corresponding to this panel object
    510 """
    511 doc = init_doc(doc)
--> 512 root = self._get_model(doc, comm=comm)
    513 if preprocess:
    514     self._preprocess(root)

File ~/miniconda3/envs/baleine/lib/python3.9/site-packages/panel/layout/base.py:122, in Panel._get_model(self, doc, root, parent, comm)
    120 if root is None:
    121     root = model
--> 122 objects = self._get_objects(model, [], doc, root, comm)
    123 props = dict(self._init_params(), objects=objects)
    124 model.update(**self._process_param_change(props))

File ~/miniconda3/envs/baleine/lib/python3.9/site-packages/panel/layout/base.py:112, in Panel._get_objects(self, model, old_objects, doc, root, comm)
    110 else:
    111     try:
--> 112         child = pane._get_model(doc, root, model, comm)
    113     except RerenderError:
    114         return self._get_objects(model, current_objects[:i], doc, root, comm)

File ~/miniconda3/envs/baleine/lib/python3.9/site-packages/panel/pane/holoviews.py:239, in HoloViews._get_model(self, doc, root, parent, comm)
    237     plot = self.object
    238 else:
--> 239     plot = self._render(doc, comm, root)
    241 plot.pane = self
    242 backend = plot.renderer.backend

File ~/miniconda3/envs/baleine/lib/python3.9/site-packages/panel/pane/holoviews.py:312, in HoloViews._render(self, doc, comm, root)
    309     if comm:
    310         kwargs['comm'] = comm
--> 312 return renderer.get_plot(self.object, **kwargs)

File ~/miniconda3/envs/baleine/lib/python3.9/site-packages/holoviews/plotting/bokeh/renderer.py:73, in BokehRenderer.get_plot(self_or_cls, obj, doc, renderer, **kwargs)
     66 @bothmethod
     67 def get_plot(self_or_cls, obj, doc=None, renderer=None, **kwargs):
     68     """
     69     Given a HoloViews Viewable return a corresponding plot instance.
     70     Allows supplying a document attach the plot to, useful when
     71     combining the bokeh model with another plot.
     72     """
---> 73     plot = super(BokehRenderer, self_or_cls).get_plot(obj, doc, renderer, **kwargs)
     74     if plot.document is None:
     75         plot.document = Document() if self_or_cls.notebook_context else curdoc()

File ~/miniconda3/envs/baleine/lib/python3.9/site-packages/holoviews/plotting/renderer.py:243, in Renderer.get_plot(self_or_cls, obj, doc, renderer, comm, **kwargs)
    240     defaults = [kd.default for kd in plot.dimensions]
    241     init_key = tuple(v if d is None else d for v, d in
    242                      zip(plot.keys[0], defaults))
--> 243     plot.update(init_key)
    244 else:
    245     plot = obj

File ~/miniconda3/envs/baleine/lib/python3.9/site-packages/holoviews/plotting/plot.py:991, in DimensionedPlot.update(self, key)
    989 def update(self, key):
    990     if len(self) == 1 and ((key == 0) or (key == self.keys[0])) and not self.drawn:
--> 991         return self.initialize_plot()
    992     item = self.__getitem__(key)
    993     self.traverse(lambda x: setattr(x, '_updated', True))

File ~/miniconda3/envs/baleine/lib/python3.9/site-packages/geoviews/plotting/bokeh/plot.py:111, in GeoPlot.initialize_plot(self, ranges, plot, plots, source)
    109 def initialize_plot(self, ranges=None, plot=None, plots=None, source=None):
    110     opts = {} if isinstance(self, HvOverlayPlot) else {'source': source}
--> 111     fig = super(GeoPlot, self).initialize_plot(ranges, plot, plots, **opts)
    112     if self.geographic and self.show_bounds and not self.overlaid:
    113         from . import GeoShapePlot

File ~/miniconda3/envs/baleine/lib/python3.9/site-packages/holoviews/plotting/bokeh/element.py:1409, in ElementPlot.initialize_plot(self, ranges, plot, plots, source)
   1407 # Initialize plot, source and glyph
   1408 if plot is None:
-> 1409     plot = self._init_plot(key, style_element, ranges=ranges, plots=plots)
   1410     self._init_axes(plot)
   1411 else:

File ~/miniconda3/envs/baleine/lib/python3.9/site-packages/holoviews/plotting/bokeh/element.py:496, in ElementPlot._init_plot(self, key, element, plots, ranges)
    489 """
    490 Initializes Bokeh figure to draw Element into and sets basic
    491 figure and axis attributes including axes types, labels,
    492 titles and plot height and width.
    493 """
    494 subplots = list(self.subplots.values()) if self.subplots else []
--> 496 axis_types, labels, plot_ranges = self._axes_props(plots, subplots, element, ranges)
    497 xlabel, ylabel, _ = labels
    498 x_axis_type, y_axis_type = axis_types

File ~/miniconda3/envs/baleine/lib/python3.9/site-packages/holoviews/plotting/bokeh/element.py:407, in ElementPlot._axes_props(self, plots, subplots, element, ranges)
    405 # Get the Element that determines the range and get_extents
    406 range_el = el if self.batched and not isinstance(self, OverlayPlot) else element
--> 407 l, b, r, t = self.get_extents(range_el, ranges)
    408 if self.invert_axes:
    409     l, b, r, t = b, l, t, r

File ~/miniconda3/envs/baleine/lib/python3.9/site-packages/geoviews/plotting/plot.py:73, in ProjectionPlot.get_extents(self, element, ranges, range_type)
     71     extents = None
     72 else:
---> 73     extents = project_extents(extents, element.crs, proj)
     74 return (np.NaN,)*4 if not extents else extents

File ~/miniconda3/envs/baleine/lib/python3.9/site-packages/geoviews/util.py:102, in project_extents(extents, src_proj, dest_proj, tol)
    100     geom_in_src_proj = geom_clipped_to_dest_proj
    101 try:
--> 102     geom_in_crs = dest_proj.project_geometry(geom_in_src_proj, src_proj)
    103 except ValueError:
    104     src_name =type(src_proj).__name__

File ~/miniconda3/envs/baleine/lib/python3.9/site-packages/cartopy/crs.py:805, in Projection.project_geometry(self, geometry, src_crs)
    803 if not method_name:
    804     raise ValueError(f'Unsupported geometry type {geom_type!r}')
--> 805 return getattr(self, method_name)(geometry, src_crs)

File ~/miniconda3/envs/baleine/lib/python3.9/site-packages/cartopy/crs.py:941, in Projection._project_polygon(self, polygon, src_crs)
    939     is_ccw = True
    940 else:
--> 941     is_ccw = polygon.exterior.is_ccw
    942 # Project the polygon exterior/interior rings.
    943 # Each source ring will result in either a ring, or one or more
    944 # lines.
    945 rings = []

File ~/miniconda3/envs/baleine/lib/python3.9/site-packages/shapely/geometry/polygon.py:99, in LinearRing.is_ccw(self)
     96 @property
     97 def is_ccw(self):
     98     """True is the ring is oriented counter clock-wise"""
---> 99     return bool(self.impl['is_ccw'](self))

File ~/miniconda3/envs/baleine/lib/python3.9/site-packages/shapely/algorithms/cga.py:14, in is_ccw_impl.<locals>.is_ccw_op(ring)
     13 def is_ccw_op(ring):
---> 14     return signed_area(ring) >= 0.0

File ~/miniconda3/envs/baleine/lib/python3.9/site-packages/shapely/algorithms/cga.py:7, in signed_area(ring)
      3 """Return the signed area enclosed by a ring in linear time using the
      4 algorithm at: https://web.archive.org/web/20080209143651/http://cgafaq.info:80/wiki/Polygon_Area
      5 """
      6 xs, ys = ring.coords.xy
----> 7 xs.append(xs[1])
      8 ys.append(ys[1])
      9 return sum(xs[i]*(ys[i+1]-ys[i-1]) for i in range(1, len(ring.coords)))/2.0

IndexError: array index out of range
:Image   [y,x]   (band_data)

the same error is happening with xr.open_dataset('bug.nc').band_data.hvplot(geo=True)

I can give you the bug.nc netcdf file if you wish.

[maximlt]

Yes it would be nice if you could share this NetCDF file (EPSG 2154, thanks for the memory 🇫🇷).

[hxawax]

Yes it would be nice if you could share this NetCDF file (EPSG 2154, thanks for the memory fr).

Thanks for your answer !
here is the NETCDF file https://lufi.ethibox.fr/r/JaKj_s6MPH#a87n9jsfTLrTGLD/1QVKTjAOeQ2bSEvudwEdcioPpzs=
(30 days link)

[kearney-sp]

I am getting the same error. Using an xarray DataArray in CRS 'epsg:3857'. Works fine with rasterize=True and Project=False. Fails with Project=True.

[maximlt]

@hxawax I can get your plot working by declaring explicitly x and y:

import xarray as xr
import hvplot.xarray

ds = xr.open_dataset('bug.nc')
ds.band_data.hvplot(x='x', y='y', crs=cartopy.crs.epsg(2154))

(I think what happens when you don't declare them explicitly is that hvPlot assumes you want to plot (y, x) since this is the convention for arrays. The axes being now messed up their extents is also messed up which lead to the error message)

[maximlt]

@kearney-sp did you experience exactly the same issue due to not declaring explicitly x and y? If not, could you please provide an example to reproduce your issue?

[kearney-sp]

@maximlt I was declaring x and y explicitly when I got the error. I will see if I can create a reproducible example, or at least just paste the code snippet (it's part of a much larger app). Thanks.

[maximlt]

Thanks, that'd be much appreciated @kearney-sp !

[gcaria]

I sometimes have the same error, below you can find a MRE, and these are the package versions I am using (please let me know if some relevant ones are missing).

I believe I only get this error with some CRSs but not with all of them.

geoviews             1.9.5
holoviews            1.14.9
xarray               2022.3.0
rioxarray            0.10.3
hvplot               0.7.3

import numpy as np
import xarray as xr
import rioxarray as rxr
import cartopy.crs as ccrs
import geoviews as gv
import holoviews as hv
import hvplot.xarray

da = xr.DataArray(np.arange(1_000_000).reshape(1000,1000),
                  coords=dict(x=np.linspace(577200, 867600, 1000),
                              y=np.linspace(-1802100, -2036700, 1000))
                 )

da = da.rio.write_crs(32635)
crs = ccrs.epsg(32635)
da.hvplot(x="x", y="y", crs=crs, rasterize=True) 

[hoxbro]

@gcaria, I'm in no way an expert in maps, but could the reason you are seeing the IndexError be because you are exceeding the map projected bounds as seen here: https://epsg.io/32635?

Because if I change the coords in your example to the following, it works:

import numpy as np
import xarray as xr
import rioxarray as rxr
import cartopy.crs as ccrs
import geoviews as gv
import holoviews as hv
import hvplot.xarray

da = xr.DataArray(np.arange(1_000_000).reshape(1000,1000),
                  coords=dict(x=np.linspace(166021, 534994, 1000),
                              y=np.linspace(0.00, 9329005, 1000))
                 )
crs = ccrs.epsg(32635)
da = da.rio.write_crs(crs)
da.hvplot(x="y", y="x", crs=crs, rasterize=True)

(Note there is a weird artifact when you move around, this will be fixed with holoviz/holoviews#5351)

@hxawax, this could also be the reason why you see the IndexError as 2154 is limited to France.

I think a more descriptive error message is needed when exceeding the bounds and getting an IndexError.

[gcaria]

Good catch @hoxbro ! There's a very good chance that's what's going on.

I agree with you that a better error message would be very helpful.

[jsignell]

This looks like the same error that is showing up in #951 (leading me to believe even more firmly that that issue has to do with bad hvplot misunderstanding the crs of the bounds).

[jsignell]

Having finally circled back to this, I think that in order to raise a more informative error there would need to be a change on the holoviews side.

[hoxbro]

Are you sure the change is needed in Holoviews?

Because I did some initial work to give a more informative error message in holoviz/geoviews#581 but never finished it.

[jsignell]

Oh sorry I mean "not on hvplot" yeah it's more likely that it would happen on geoviews.

[hoxbro]

You are more than welcome to continue the work on the PR or rewrite it if that is easier 😅