search

holoviz / geoviews

Simple, concise geographical visualization in Python
http://geoviews.org
BSD 3-Clause "New" or "Revised" License
587 stars 75 forks source link

gv.Polygons (default)PlateCarree IndexError #577

Open DmitriyLeybel opened 2 years ago

DmitriyLeybel commented 2 years ago

ALL software version info

python == 3.9.13

geoviews == 1.9.5 holoviews == 1.14.8 shapely == 1.8.2 cartopy == 0.20.3 bokeh == 2.4.2 matplotlib == 3.5.1

Description of expected behavior and the observed behavior

In attempting to plot a geometry from a shape file from the census.gov repo, I'm experiencing errors when using the default crs argument(Plate Carree) of gv.Polygons. The same error persists when assigning Plate Carree as the crs to use with hvplot, when changing the plotting backend from bokeh to matplotlib, and when changing the crs with geopandas(to_crs). Other crs projections do appear to work when set in the plotting function, such as Mercator, GOOGLE_MERCATOR, and EckertII, and default hvplot method.

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

From Notebook:

import geopandas as gpd
import geoviews as gv
from cartopy import crs

gv.extension('bokeh')

t = gpd.read_file('test_geo.json')

gv.Polygons(t) #Index Error 
# t.hvplot(crs=crs.PlateCarree()) #Another Index Error

# gv.Polygons(t, crs=crs.Mercator()) # Works
# t.hvplot() # Works
IndexError ```python IndexError Traceback (most recent call last) File ~\.conda\envs\geo\lib\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 ~\.conda\envs\geo\lib\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 ~\.conda\envs\geo\lib\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 ~\.conda\envs\geo\lib\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 ~\.conda\envs\geo\lib\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 ~\.conda\envs\geo\lib\site-packages\holoviews\ipython\display_hooks.py:146, in display_hook..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 ~\.conda\envs\geo\lib\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 ~\.conda\envs\geo\lib\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 ~\.conda\envs\geo\lib\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 ~\.conda\envs\geo\lib\site-packages\panel\viewable.py:460, in Renderable._render_model(self, doc, comm) 458 if comm is None: 459 comm = state._comm_manager.get_server_comm() --> 460 model = self.get_root(doc, comm) 462 if config.embed: 463 embed_state(self, model, doc, 464 json=config.embed_json, 465 json_prefix=config.embed_json_prefix, 466 save_path=config.embed_save_path, 467 load_path=config.embed_load_path, 468 progress=False) File ~\.conda\envs\geo\lib\site-packages\panel\viewable.py:518, in Renderable.get_root(self, doc, comm, preprocess) 501 """ 502 Returns the root model and applies pre-processing hooks 503 (...) 515 Returns the bokeh model corresponding to this panel object 516 """ 517 doc = init_doc(doc) --> 518 root = self._get_model(doc, comm=comm) 519 if preprocess: 520 self._preprocess(root) File ~\.conda\envs\geo\lib\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 ~\.conda\envs\geo\lib\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 ~\.conda\envs\geo\lib\site-packages\panel\pane\holoviews.py:248, in HoloViews._get_model(self, doc, root, parent, comm) 246 plot = self.object 247 else: --> 248 plot = self._render(doc, comm, root) 250 plot.pane = self 251 backend = plot.renderer.backend File ~\.conda\envs\geo\lib\site-packages\panel\pane\holoviews.py:321, in HoloViews._render(self, doc, comm, root) 318 if comm: 319 kwargs['comm'] = comm --> 321 return renderer.get_plot(self.object, **kwargs) File ~\.conda\envs\geo\lib\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 ~\.conda\envs\geo\lib\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 ~\.conda\envs\geo\lib\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 ~\.conda\envs\geo\lib\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 ~\.conda\envs\geo\lib\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 ~\.conda\envs\geo\lib\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 ~\.conda\envs\geo\lib\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 ~\.conda\envs\geo\lib\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 ~\.conda\envs\geo\lib\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 ~\.conda\envs\geo\lib\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 ~\.conda\envs\geo\lib\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 ~\.conda\envs\geo\lib\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 ~\.conda\envs\geo\lib\site-packages\shapely\algorithms\cga.py:14, in is_ccw_impl..is_ccw_op(ring) 13 def is_ccw_op(ring): ---> 14 return signed_area(ring) >= 0.0 File ~\.conda\envs\geo\lib\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 :Polygons [Longitude,Latitude] (id,name) ```

test_geo.zip

gv.Polygons with Mercator crs: image

hoxbro commented 2 years ago

gv.Polygons defaults to PlateCarree, which is why both of those examples give an IndexError. hvplot without any inputs works because it only outputs an Holoviews plot. If you set geo=True it will also give an IndexError.

I think you get an IndexError because your data is projected coordinates with a geographic coordinates system.

Your solution of setting crs looks fine to me. Alternative, you can change the CRS of the GeoDataFrame with t = t.set_crs("EPSG:3857", allow_override=True).to_crs("EPSG:4326") (not sure if those are the correct CRS).

When #581 is merged, the error will be better.