holoviz / geoviews

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

TriMesh plot fails when global triangular mesh is used. #318

Open koldunovn opened 5 years ago

koldunovn commented 5 years ago

I am trying to follow this example http://geo.holoviews.org/gallery/bokeh/trimesh_uk.html#bokeh-gallery-trimesh-uk

but plot triangular mesh that covers the whole globe. Below are the code and the error message (I run the code in classical Jupyter notebook):

import numpy as np
import geoviews as gv
import geoviews.feature as gf
import cartopy.crs as ccrs
import pandas as pd
gv.extension('bokeh')

#download the data (less than Mb)
!wget https://github.com/koldunovn/pi_mesh_python/raw/master/xnodes.txt
!wget https://github.com/koldunovn/pi_mesh_python/raw/master/ynodes.txt
!wget https://github.com/koldunovn/pi_mesh_python/raw/master/elements.txt

xnodes = np.loadtxt('./xnodes.txt')
ynodes = np.loadtxt('./ynodes.txt')
elem   = np.loadtxt('./elements.txt').astype('int')

#this part went well
trimesh = gv.TriMesh((elem, (xnodes, ynodes)), crs=ccrs.PlateCarree())

#here the problem starts
trimesh
---------------------------------------------------------------------------
ValueError                                Traceback (most recent call last)
~/miniconda3/envs/pyfesom2/lib/python3.6/site-packages/IPython/core/formatters.py in __call__(self, obj, include, exclude)
    968 
    969             if method is not None:
--> 970                 return method(include=include, exclude=exclude)
    971             return None
    972         else:

~/miniconda3/envs/pyfesom2/lib/python3.6/site-packages/holoviews/core/dimension.py in _repr_mimebundle_(self, include, exclude)
   1291         combined and returned.
   1292         """
-> 1293         return Store.render(self)
   1294 
   1295 

~/miniconda3/envs/pyfesom2/lib/python3.6/site-packages/holoviews/core/options.py in render(cls, obj)
   1360         data, metadata = {}, {}
   1361         for hook in hooks:
-> 1362             ret = hook(obj)
   1363             if ret is None:
   1364                 continue

~/miniconda3/envs/pyfesom2/lib/python3.6/site-packages/holoviews/ipython/display_hooks.py in pprint_display(obj)
    275     if not ip.display_formatter.formatters['text/plain'].pprint:
    276         return None
--> 277     return display(obj, raw_output=True)
    278 
    279 

~/miniconda3/envs/pyfesom2/lib/python3.6/site-packages/holoviews/ipython/display_hooks.py in display(obj, raw_output, **kwargs)
    245     elif isinstance(obj, (CompositeOverlay, ViewableElement)):
    246         with option_state(obj):
--> 247             output = element_display(obj)
    248     elif isinstance(obj, (Layout, NdLayout, AdjointLayout)):
    249         with option_state(obj):

~/miniconda3/envs/pyfesom2/lib/python3.6/site-packages/holoviews/ipython/display_hooks.py in wrapped(element)
    140         try:
    141             max_frames = OutputSettings.options['max_frames']
--> 142             mimebundle = fn(element, max_frames=max_frames)
    143             if mimebundle is None:
    144                 return {}, {}

~/miniconda3/envs/pyfesom2/lib/python3.6/site-packages/holoviews/ipython/display_hooks.py in element_display(element, max_frames)
    186         return None
    187 
--> 188     return render(element)
    189 
    190 

~/miniconda3/envs/pyfesom2/lib/python3.6/site-packages/holoviews/ipython/display_hooks.py in render(obj, **kwargs)
     63         renderer = renderer.instance(fig='png')
     64 
---> 65     return renderer.components(obj, **kwargs)
     66 
     67 

~/miniconda3/envs/pyfesom2/lib/python3.6/site-packages/holoviews/plotting/bokeh/renderer.py in components(self, obj, fmt, comm, **kwargs)
    247         # Bokeh has to handle comms directly in <0.12.15
    248         comm = False if bokeh_version < '0.12.15' else comm
--> 249         return super(BokehRenderer, self).components(obj,fmt, comm, **kwargs)
    250 
    251 

~/miniconda3/envs/pyfesom2/lib/python3.6/site-packages/holoviews/plotting/renderer.py in components(self, obj, fmt, comm, **kwargs)
    319             plot = obj
    320         else:
--> 321             plot, fmt = self._validate(obj, fmt)
    322 
    323         widget_id = None

~/miniconda3/envs/pyfesom2/lib/python3.6/site-packages/holoviews/plotting/renderer.py in _validate(self, obj, fmt, **kwargs)
    218         if isinstance(obj, tuple(self.widgets.values())):
    219             return obj, 'html'
--> 220         plot = self.get_plot(obj, renderer=self, **kwargs)
    221 
    222         fig_formats = self.mode_formats['fig'][self.mode]

~/miniconda3/envs/pyfesom2/lib/python3.6/site-packages/holoviews/plotting/bokeh/renderer.py in get_plot(self_or_cls, obj, doc, renderer, **kwargs)
    132             curdoc().theme = self_or_cls.theme
    133         doc.theme = self_or_cls.theme
--> 134         plot = super(BokehRenderer, self_or_cls).get_plot(obj, renderer, **kwargs)
    135         plot.document = doc
    136         return plot

~/miniconda3/envs/pyfesom2/lib/python3.6/site-packages/holoviews/plotting/renderer.py in get_plot(self_or_cls, obj, renderer, **kwargs)
    205             init_key = tuple(v if d is None else d for v, d in
    206                              zip(plot.keys[0], defaults))
--> 207             plot.update(init_key)
    208         else:
    209             plot = obj

~/miniconda3/envs/pyfesom2/lib/python3.6/site-packages/holoviews/plotting/plot.py in update(self, key)
    593     def update(self, key):
    594         if len(self) == 1 and ((key == 0) or (key == self.keys[0])) and not self.drawn:
--> 595             return self.initialize_plot()
    596         item = self.__getitem__(key)
    597         self.traverse(lambda x: setattr(x, '_updated', True))

~/miniconda3/envs/pyfesom2/lib/python3.6/site-packages/geoviews/plotting/bokeh/plot.py in initialize_plot(self, ranges, plot, plots, source)
     86     def initialize_plot(self, ranges=None, plot=None, plots=None, source=None):
     87         opts = {} if isinstance(self, HvOverlayPlot) else {'source': source}
---> 88         fig = super(GeoPlot, self).initialize_plot(ranges, plot, plots, **opts)
     89         if self.geographic and self.show_bounds and not self.overlaid:
     90             from . import GeoShapePlot

~/miniconda3/envs/pyfesom2/lib/python3.6/site-packages/holoviews/plotting/bokeh/element.py in initialize_plot(self, ranges, plot, plots, source)
   1240         self.handles['plot'] = plot
   1241 
-> 1242         self._init_glyphs(plot, element, ranges, source)
   1243         if not self.overlaid:
   1244             self._update_plot(key, plot, style_element)

~/miniconda3/envs/pyfesom2/lib/python3.6/site-packages/holoviews/plotting/bokeh/graphs.py in _init_glyphs(self, plot, element, ranges, source)
    523     def _init_glyphs(self, plot, element, ranges, source):
    524         element = self._process_vertices(element)
--> 525         super(TriMeshPlot, self)._init_glyphs(plot, element, ranges, source)
    526 
    527     def _update_glyphs(self, element, ranges, style):

~/miniconda3/envs/pyfesom2/lib/python3.6/site-packages/holoviews/plotting/bokeh/graphs.py in _init_glyphs(self, plot, element, ranges, source)
    337         # Get data and initialize data source
    338         style = self.style[self.cyclic_index]
--> 339         data, mapping, style = self.get_data(element, ranges, style)
    340         self.handles['previous_id'] = element._plot_id
    341 

~/miniconda3/envs/pyfesom2/lib/python3.6/site-packages/geoviews/plotting/bokeh/plot.py in get_data(self, element, ranges, style)
    132         if self._project_operation and self.geographic:
    133             element = self._project_operation(element, projection=proj)
--> 134         return super(GeoPlot, self).get_data(element, ranges, style)
    135 
    136 

~/miniconda3/envs/pyfesom2/lib/python3.6/site-packages/holoviews/plotting/bokeh/graphs.py in get_data(self, element, ranges, style)
    225         self._get_edge_colors(element, ranges, path_data, edge_mapping, style)
    226         if not static:
--> 227             pdata, pmapping = self._get_edge_paths(element, ranges)
    228             path_data.update(pdata)
    229             edge_mapping.update(pmapping)

~/miniconda3/envs/pyfesom2/lib/python3.6/site-packages/holoviews/plotting/bokeh/graphs.py in _get_edge_paths(self, element, ranges)
    160             else:
    161                 raise ValueError("Edge paths do not match the number of supplied edges."
--> 162                                  "Expected %d, found %d paths." % (len(element), len(edges)))
    163         elif self.directed:
    164             xdim, ydim = element.nodes.kdims[:2]

ValueError: Edge paths do not match the number of supplied edges.Expected 5839, found 5793 paths.

:TriMesh   [node1,node2,node3]

If I use a small subset of elements, for example:

trimesh = gv.TriMesh((elem[:1000,:], (xnodes, ynodes)), crs=ccrs.PlateCarree())

The plot is shown covered by points at nodes locations that are connected in some places (as they should be).

I would appreciate any help.

BTW the mesh should look something like this alt text

philippjfr commented 5 years ago

Suspect something goes wrong at the boundaries when projecting the TriMesh. Will have to have a more detailed look later.

koldunovn commented 5 years ago

Thanks a lot for looking at this!

In order to avoid problems at the boundaries we usually just remove elements that cross +-180, with something like this:

d=xnodes[elem].max(axis=1) - xnodes[elem].min(axis=1)
no_cyclic_elem = [i for (i, val) in enumerate(d) if val < 100]

But when I use:

trimesh = gv.TriMesh((elem[no_cyclic_elem], (xnodes, ynodes)), crs=ccrs.PlateCarree())
trimesh

the error is the same, while in Basemap and cartopy it usually solves the problem of horizontal stripes that go across the globe.

koldunovn commented 5 years ago

One more thing that might help is that I have discovered the code below works (and OMG it's fast!):

from holoviews.operation.datashader import datashade

trimesh = gv.TriMesh((elem[no_cyclic_elem], (xnodes, ynodes)), crs=ccrs.PlateCarree())
datashade(trimesh ) 
sameerCoder commented 5 years ago

Hi koldunovn any solution you got for this ?

koldunovn commented 5 years ago

@sameerCoder No, no progress on this.

sameerCoder commented 5 years ago

Hi Koldunovn can you help me in solving getting colorbar and colormap in my trimesh plot . plz see the this question of stackoverflow and let me know how to solve this issue . thank you.