Added projection aware QuadMesh element

doc/Gridded_Datasets_II.ipynb

@@ -56,9 +56,7 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "kdims = ['realization', 'longitude', 'latitude', 'time']\n",
-    "vdims = ['surface_temperature']\n",
-    "dataset = gv.Dataset(xr_ensembles, kdims=kdims, vdims=vdims, crs=ccrs.PlateCarree())\n",
+    "dataset = gv.Dataset(xr_ensembles, vdims='surface_temperature', crs=ccrs.PlateCarree())\n",
     "dataset"
    ]
   },
@@ -143,7 +141,27 @@
    "source": [
     "Using ``dynamic`` mode means that the data for each frame is only extracted when you're actually viewing that part of the data, which can have huge benefits in terms of speed and memory consumption.  However, it relies on having a running Python process to render and serve each image, and so it cannot be used when generating static HTML output such as for the GeoViews web site.\n",
     "\n",
+    "## Irregularly sampled data\n",
     "\n",
+    "Often gridded datasets are not regularly sampled, instead providing irregularly sampled multi-dimensional coordinates. Such datasets can be easily visualized using the QuadMesh element."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "%%opts QuadMesh [colorbar=True fig_size=300 projection=ccrs.Robinson()] (cmap='RdBu_r')\n",
+    "xrds = xr.tutorial.load_dataset('RASM_example_data')\n",
+    "qmesh = gv.Dataset(xrds.Tair).to(gv.QuadMesh, ['xc', 'yc'], dynamic=True)\n",
+    "qmesh.redim.range(Tair=(-30, 30)) * gf.coastline"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
     "## Non-geographic views\n",
     "\n",
     "So far we have focused entirely on geographic views of the data, plotting the data on a projection. However the conversion interface is completely general, allowing us to slice and dice the data in any way we like. For these views we will switch to the bokeh plotting extension:"

geoviews/__init__.py

@@ -9,7 +9,8 @@
 from .element import (_Element, Feature, Tiles,     # noqa (API import)
                       WMTS, LineContours, FilledContours, Text, Image,
                       Points, Path, Polygons, Shape, Dataset, RGB,
-                      Contours, Graph, TriMesh, Nodes, EdgePaths )
+                      Contours, Graph, TriMesh, Nodes, EdgePaths,
+                      QuadMesh)
 from . import data                                  # noqa (API import)
 from . import operation                             # noqa (API import)
 from . import plotting                              # noqa (API import)

geoviews/element/__init__.py

@@ -3,7 +3,7 @@
 from .geo import (_Element, Feature, Tiles, is_geographic,     # noqa (API import)
                   WMTS, Points, Image, Text, LineContours, RGB,
                   FilledContours, Path, Polygons, Shape, Dataset,
-                  Contours, TriMesh, Graph, Nodes, EdgePaths)
+                  Contours, TriMesh, Graph, Nodes, EdgePaths, QuadMesh)
 
 
 class GeoConversion(ElementConversion):

geoviews/element/geo.py

@@ -5,11 +5,13 @@
 from cartopy.io.img_tiles import GoogleTiles as cGoogleTiles
 from cartopy.io.shapereader import Reader
 from holoviews.core import Element2D, Dimension, Dataset as HvDataset, NdOverlay
-from holoviews.core.util import basestring, pd, max_extents, dimension_range
+from holoviews.core.util import (basestring, pd, max_extents,
+                                 dimension_range, get_param_values)
 from holoviews.element import (
     Contours as HvContours, Graph as HvGraph, Image as HvImage,
     Nodes as HvNodes, Path as HvPath, Polygons as HvPolygons,
-    RGB as HvRGB, Text as HvText, TriMesh as HvTriMesh)
+    RGB as HvRGB, Text as HvText, TriMesh as HvTriMesh,
+    QuadMesh as HvQuadMesh)
 
 from shapely.geometry.base import BaseGeometry
 
@@ -238,6 +240,36 @@ class Image(_Element, HvImage):
     group = param.String(default='Image')
 
 
+class QuadMesh(_Element, HvQuadMesh):
+    """
+    QuadMesh is a Raster type to hold x- and y- bin values
+    with associated values. The x- and y-values of the QuadMesh
+    may be supplied either as the edges of each bin allowing
+    uneven sampling or as the bin centers, which will be converted
+    to evenly sampled edges.
+
+    As a secondary but less supported mode QuadMesh can contain
+    a mesh of quadrilateral coordinates that is not laid out in
+    a grid. The data should then be supplied as three separate
+    2D arrays for the x-/y-coordinates and grid values.
+    """
+
+    datatype = param.List(default=['grid', 'xarray'])
+
+    vdims = param.List(default=[Dimension('z')], bounds=(1, 1))
+
+    group = param.String(default='QuadMesh')
+
+    _binned = True
+
+    def trimesh(self):
+        trimesh = super(QuadMesh, self).trimesh()
+        node_params = get_param_values(trimesh.nodes)
+        nodes = TriMesh.node_type(trimesh.nodes.data, **node_params)
+        return TriMesh((trimesh.data, nodes), crs=self.crs,
+                       **get_param_values(trimesh))
+
+
 class RGB(_Element, HvRGB):
     """
     An RGB element is a Image containing channel data for the the

geoviews/operation.py

@@ -3,11 +3,12 @@
 
 from cartopy import crs as ccrs
 from cartopy.img_transform import warp_array, _determine_bounds
-from holoviews.core.util import cartesian_product
+from holoviews.core.util import cartesian_product, get_param_values
 from holoviews.operation import Operation
 from shapely.geometry import Polygon, LineString
 
-from .element import Image, Shape, Polygons, Path, Points, Contours, RGB, Graph, Nodes, EdgePaths
+from .element import (Image, Shape, Polygons, Path, Points, Contours,
+                      RGB, Graph, Nodes, EdgePaths, QuadMesh)
 from .util import project_extents, geom_to_array
 
 
@@ -101,6 +102,58 @@ def _process_element(self, element):
         return element.clone(data, crs=self.projection)
 
 
+class project_quadmesh(_project_operation):
+
+    supported_types = [QuadMesh]
+
+    def _process_element(self, element):
+        proj = self.p.projection
+        irregular = any(element.interface.irregular(element, kd)
+                        for kd in element.kdims)
+        zs = element.dimension_values(2, flat=False)
+        if irregular:
+            X, Y = [np.asarray(element.interface.coords(element, kd, expanded=True))
+                    for kd in element.kdims]
+        else:
+            X = element.dimension_values(0, expanded=True)
+            Y = element.dimension_values(1, expanded=True)
+            zs = zs.T
+
+        coords = proj.transform_points(element.crs, X, Y)
+        PX, PY = coords[..., 0], coords[..., 1]
+
+        # Mask quads which are wrapping around the x-axis
+        wrap_proj_types = (ccrs._RectangularProjection,
+                           ccrs._WarpedRectangularProjection,
+                           ccrs.InterruptedGoodeHomolosine,
+                           ccrs.Mercator)
+        if isinstance(proj, wrap_proj_types):
+            with np.errstate(invalid='ignore'):
+                edge_lengths = np.hypot(
+                    np.diff(PX , axis=1),
+                    np.diff(PY, axis=1)
+                )
+                to_mask = (
+                    (edge_lengths >= abs(proj.x_limits[1] -
+                                         proj.x_limits[0]) / 2) |
+                    np.isnan(edge_lengths)
+                )
+            if np.any(to_mask):
+                mask = np.zeros(zs.shape, dtype=np.bool)
+                mask[:, 1:][to_mask] = True
+                mask[:, 2:][to_mask[:, :-1]] = True
+                mask[:, :-1][to_mask] = True
+                mask[:, :-2][to_mask[:, 1:]] = True
+                mask[1:, 1:][to_mask[:-1]] = True
+                mask[1:, :-1][to_mask[:-1]] = True
+                mask[:-1, 1:][to_mask[1:]] = True
+                mask[:-1, :-1][to_mask[1:]] = True
+                zs[mask] = np.NaN
+
+        params = get_param_values(element)
+        return QuadMesh((PX, PY, zs), crs=self.projection, **params)
+
+
 class project_image(_project_operation):
     """
     Projects an geoviews Image to the specified projection,
@@ -229,4 +282,5 @@ def _process(self, element, key=None):
         element = element.map(project_image, project_image.supported_types)
         element = element.map(project_shape, project_shape.supported_types)
         element = element.map(project_graph, project_graph.supported_types)
+        element = element.map(project_quadmesh, project_quadmesh.supported_types)
         return element.map(project_points, project_points.supported_types)

geoviews/plotting/bokeh/__init__.py

@@ -11,13 +11,13 @@
 from holoviews.plotting.bokeh.chart import PointPlot
 from holoviews.plotting.bokeh.graphs import TriMeshPlot, GraphPlot
 from holoviews.plotting.bokeh.path import PolygonPlot, PathPlot, ContourPlot
-from holoviews.plotting.bokeh.raster import RasterPlot, RGBPlot
+from holoviews.plotting.bokeh.raster import RasterPlot, RGBPlot, QuadMeshPlot
 
 from ...element import (WMTS, Points, Polygons, Path, Contours, Shape,
                         Image, Feature, Text, RGB, Nodes, EdgePaths,
-                        Graph, TriMesh)
+                        Graph, TriMesh, QuadMesh)
 from ...operation import (project_image, project_shape, project_points,
-                          project_path, project_graph)
+                          project_path, project_graph, project_quadmesh)
 from ...util import geom_to_array
 from .plot import GeoPlot, OverlayPlot, DEFAULT_PROJ
 from . import callbacks # noqa
@@ -76,6 +76,11 @@ class GeoPointPlot(GeoPlot, PointPlot):
     _project_operation = project_points
 
 
+class GeoQuadMeshPlot(GeoPlot, QuadMeshPlot):
+
+    _project_operation = project_quadmesh
+
+
 class GeoRasterPlot(GeoPlot, RasterPlot):
 
     _project_operation = project_image.instance(fast=False)
@@ -200,7 +205,8 @@ def get_extents(self, element, ranges=None):
                 Graph: GeoGraphPlot,
                 TriMesh: GeoTriMeshPlot,
                 Nodes: GeoPointPlot,
-                EdgePaths: GeoPathPlot}, 'bokeh')
+                EdgePaths: GeoPathPlot,
+                QuadMesh: GeoQuadMeshPlot}, 'bokeh')
 
 options = Store.options(backend='bokeh')
 

geoviews/plotting/mpl/__init__.py

@@ -12,24 +12,25 @@
 from holoviews.core import (Store, HoloMap, Layout, Overlay,
                             CompositeOverlay, Element, NdLayout)
 from holoviews.core import util
-
+from holoviews.core.data import GridInterface
 from holoviews.core.options import SkipRendering, Options
 from holoviews.plotting.mpl import (ElementPlot, ColorbarPlot, PointPlot,
                                     AnnotationPlot, TextPlot,
                                     LayoutPlot as HvLayoutPlot,
                                     OverlayPlot as HvOverlayPlot,
                                     PathPlot, PolygonPlot, ImagePlot,
-                                    ContourPlot, GraphPlot, TriMeshPlot)
+                                    ContourPlot, GraphPlot, TriMeshPlot,
+                                    QuadMeshPlot)
 from holoviews.plotting.mpl.util import get_raster_array
 
 
 from ...element import (Image, Points, Feature, WMTS, Tiles, Text,
                         LineContours, FilledContours, is_geographic,
                         Path, Polygons, Shape, RGB, Contours, Nodes,
-                        EdgePaths, Graph, TriMesh)
+                        EdgePaths, Graph, TriMesh, QuadMesh)
 from ...util import project_extents, geo_mesh
 
-from ...operation import project_points, project_path, project_graph
+from ...operation import project_points, project_path, project_graph, project_quadmesh
 
 
 def _get_projection(el):
@@ -237,6 +238,8 @@ def get_data(self, element, ranges, style):
         self._norm_kwargs(element, ranges, style, element.vdims[0])
         style.pop('interpolation', None)
         xs, ys, zs = geo_mesh(element)
+        xs = GridInterface._infer_interval_breaks(xs)
+        ys = GridInterface._infer_interval_breaks(ys)
         if self.geographic:
             style['transform'] = element.crs
         return (xs, ys, zs), style, {}
@@ -254,6 +257,18 @@ def update_handles(self, *args):
         return GeoPlot.update_handles(self, *args)
 
 
+
+class GeoQuadMeshPlot(GeoPlot, QuadMeshPlot):
+
+    _project_operation = project_quadmesh
+
+    def get_data(self, element, ranges, style):
+        if self._project_operation and self.geographic:
+            element = self._project_operation(element, projection=self.projection)
+        return super(GeoPlot, self).get_data(element, ranges, style)
+
+
+
 class GeoRGBPlot(GeoImagePlot):
     """
     Draws a imshow plot from the data in a RGB Element.
@@ -509,7 +524,8 @@ def draw_annotation(self, axis, data, crs, opts):
                 Graph: GeoGraphPlot,
                 TriMesh: GeoTriMeshPlot,
                 Nodes: GeoPointPlot,
-                EdgePaths: GeoPathPlot}, 'matplotlib')
+                EdgePaths: GeoPathPlot,
+                QuadMesh: GeoQuadMeshPlot}, 'matplotlib')
 
 
 # Define plot and style options