diff --git a/geoviews/operation/projection.py b/geoviews/operation/projection.py index 75975361..54f77204 100644 --- a/geoviews/operation/projection.py +++ b/geoviews/operation/projection.py @@ -5,7 +5,7 @@ from cartopy.img_transform import warp_array, _determine_bounds from holoviews.core.util import cartesian_product, get_param_values from holoviews.operation import Operation -from shapely.geometry import Polygon, LineString +from shapely.geometry import Polygon, LineString, MultiPolygon, MultiLineString from ..element import (Image, Shape, Polygons, Path, Points, Contours, RGB, Graph, Nodes, EdgePaths, QuadMesh, VectorField, @@ -54,7 +54,10 @@ def _process_element(self, element): boundary_poly = element.crs.project_geometry(Polygon(self.p.projection.boundary), self.p.projection) - geom_type = Polygon if isinstance(element, Polygons) else LineString + if isinstance(element, Polygons): + multi_type, geom_type = MultiPolygon, Polygon + else: + multi_type, geom_type = MultiLineString, LineString xdim, ydim = element.kdims[:2] projected = [] for path in element.split(): @@ -83,18 +86,39 @@ def _process_element(self, element): else: # Handle iso-contour case data = {k: vals[0] for k, vals in data.items()} + + # Wrap longitudes if isinstance(element.crs, ccrs.PlateCarree): vertices = wrap_path_data(path.array([0, 1]), element.crs, element.crs) geom = type(element)([vertices]).geom() else: geom = path.geom() - if not geom: + # Clip path to projection boundaries + geoms = [] + for g in geom: + if np.isinf(np.array(g.array_interface_base['data'])).sum(): + continue + try: + g = g.intersection(boundary_poly) + if 'Multi' in g.geom_type or 'Collection' in g.geom_type: + geoms += list(g) + else: + geoms.append(g) + except: + pass + + if not geoms: continue + geom = multi_type(geoms) + + # Project geometry proj = self.p.projection.project_geometry(geom, element.crs) for geom in proj: - xs, ys = np.array(geom.array_interface_base['data']).reshape(-1, 2).T - projected.append(dict(data, **{xdim.name: xs, ydim.name: ys})) + vertices = np.array(geom.array_interface_base['data']).reshape(-1, 2) + xs, ys = vertices.T + if len(xs): + projected.append(dict(data, **{xdim.name: xs, ydim.name: ys})) return element.clone(projected, crs=self.p.projection) @@ -199,6 +223,12 @@ def _process_element(self, element): PX = PX.reshape(zs.shape) if PY.ndim < 2: PY = PY.reshape(zs.shape) + + xinf, yinf = (np.argwhere((~np.isinf(CS)).sum(axis=i).astype(bool)) + for i, CS in enumerate([PX, PY])) + x0, x1 = xinf.min(), xinf.max()+1 + y0, y1 = yinf.min(), yinf.max()+1 + PX, PY, zs = PX[y0:y1, x0:x1], PY[y0:y1, x0:x1], zs[y0:y1, x0:x1] return QuadMesh((PX, PY, zs), crs=self.projection, **params) diff --git a/geoviews/util.py b/geoviews/util.py index b2924331..b9ce38cc 100644 --- a/geoviews/util.py +++ b/geoviews/util.py @@ -71,7 +71,7 @@ def path_to_geom(path, multi=True): for i, path in enumerate(paths): if i != (len(paths)-1): path = path[:-1] - if not len(path): + if len(path) < 2: continue lines.append(LineString(path[:, :2])) continue @@ -101,6 +101,8 @@ def polygon_to_geom(poly, multi=True): for i, path in enumerate(paths): if i != (len(paths)-1): path = path[:-1] + if len(path) < 3: + continue lines.append(Polygon(path[:, :2])) continue elif path.geom_type == 'MultiLineString': @@ -177,7 +179,7 @@ def geo_mesh(element): return xs, ys, zs -def wrap_data(vertices, src_crs, tgt_crs): +def wrap_path_data(vertices, src_crs, tgt_crs): """ Wraps path coordinates along the longitudinal axis. """