From 8f4e644f2813381f51c0074159e428e61e51d9f2 Mon Sep 17 00:00:00 2001 From: Philipp Rudiger Date: Sat, 19 May 2018 15:00:08 +0100 Subject: [PATCH] Fixes for projecting path data --- geoviews/element/geo.py | 1 + geoviews/operation/projection.py | 28 ++++++++++++++------- geoviews/plotting/bokeh/plot.py | 4 +-- geoviews/util.py | 42 ++++++++++++++++++++++++++++++++ 4 files changed, 64 insertions(+), 11 deletions(-) diff --git a/geoviews/element/geo.py b/geoviews/element/geo.py index 6ecfca8b..36442acc 100644 --- a/geoviews/element/geo.py +++ b/geoviews/element/geo.py @@ -285,6 +285,7 @@ class QuadMesh(_Element, HvQuadMesh): def trimesh(self): trimesh = super(QuadMesh, self).trimesh() node_params = get_param_values(trimesh.nodes) + node_params['crs'] = self.crs nodes = TriMesh.node_type(trimesh.nodes.data, **node_params) return TriMesh((trimesh.data, nodes), crs=self.crs, **get_param_values(trimesh)) diff --git a/geoviews/operation/projection.py b/geoviews/operation/projection.py index d2ba63ca..75975361 100644 --- a/geoviews/operation/projection.py +++ b/geoviews/operation/projection.py @@ -10,7 +10,7 @@ from ..element import (Image, Shape, Polygons, Path, Points, Contours, RGB, Graph, Nodes, EdgePaths, QuadMesh, VectorField, HexTiles, Labels) -from ..util import project_extents, geom_to_array +from ..util import project_extents, geom_to_array, wrap_path_data class _project_operation(Operation): @@ -65,19 +65,30 @@ def _process_element(self, element): xs, ys = geom[xdim.name], geom[ydim.name] if not len(xs): continue - path = geom_type(np.column_stack([xs, ys])) - path = path.intersection(boundary_poly) - proj = self.p.projection.project_geometry(path, element.crs) - proj_arr = geom_to_array(proj) + + proj_arr = self.p.projection.quick_vertices_transform( + np.column_stack([xs, ys]), element.crs) + + if proj_arr is None: + vertices = np.column_stack([xs, ys]) + if isinstance(element.crs, ccrs.PlateCarree): + vertices = wrap_path_data(vertices, element.crs, element.crs) + path = geom_type(vertices) + path = path.intersection(boundary_poly) + proj = self.p.projection.project_geometry(path, element.crs) + proj_arr = geom_to_array(proj) geom[xdim.name] = proj_arr[:, 0] geom[ydim.name] = proj_arr[:, 1] projected.append(geom) else: # Handle iso-contour case data = {k: vals[0] for k, vals in data.items()} - geom = path.geom() - if boundary_poly: - geom = geom.intersection(boundary_poly) + 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: continue proj = self.p.projection.project_geometry(geom, element.crs) @@ -87,7 +98,6 @@ def _process_element(self, element): return element.clone(projected, crs=self.p.projection) - class project_shape(_project_operation): """ Projects Shape Element from the source coordinate reference system diff --git a/geoviews/plotting/bokeh/plot.py b/geoviews/plotting/bokeh/plot.py index 46feaf90..d50c313c 100644 --- a/geoviews/plotting/bokeh/plot.py +++ b/geoviews/plotting/bokeh/plot.py @@ -98,7 +98,7 @@ def _postprocess_hover(self, renderer, source): except: CustomJSHover = None if (not self.geographic or None in (hover, CustomJSHover) or - isinstance(hover.tooltips, basestring)): + isinstance(hover.tooltips, basestring) or self.projection is not GOOGLE_MERCATOR): return element = self.current_frame xdim, ydim = [dimension_sanitizer(kd.name) for kd in element.kdims] @@ -147,7 +147,7 @@ def get_extents(self, element, ranges): def get_data(self, element, ranges, style): proj = self.projection - if self._project_operation and self.geographic and element.crs != proj: + if self._project_operation and self.geographic: element = self._project_operation(element, projection=proj) return super(GeoPlot, self).get_data(element, ranges, style) diff --git a/geoviews/util.py b/geoviews/util.py index 59fa524a..b2924331 100644 --- a/geoviews/util.py +++ b/geoviews/util.py @@ -1,4 +1,6 @@ import numpy as np +import shapely.geometry as sgeom + from cartopy import crs as ccrs from shapely.geometry import (MultiLineString, LineString, MultiPolygon, Polygon) @@ -114,9 +116,19 @@ def polygon_to_geom(poly, multi=True): else: path = path lines.append(path) + lines = [to_ccw(line) for line in lines] return MultiPolygon(lines) if multi else lines +def to_ccw(geom): + """ + Reorients polygon to be wound counter-clockwise. + """ + if isinstance(geom, sgeom.Polygon) and not geom.exterior.is_ccw: + geom = sgeom.polygon.orient(geom) + return geom + + def geom_to_arr(geom): arr = geom.array_interface_base['data'] if (len(arr) % 2) != 0: @@ -164,3 +176,33 @@ def geo_mesh(element): zs = np.ma.concatenate([zs, zs[:, 0:1]], axis=1) return xs, ys, zs + +def wrap_data(vertices, src_crs, tgt_crs): + """ + Wraps path coordinates along the longitudinal axis. + """ + self_params = tgt_crs.proj4_params.copy() + src_params = src_crs.proj4_params.copy() + self_params.pop('lon_0'), src_params.pop('lon_0') + + xs, ys = vertices[:, 0], vertices[:, 1] + potential = (self_params == src_params and + tgt_crs.y_limits[0] <= ys.min() and + tgt_crs.y_limits[1] >= ys.max()) + if not potential: + return vertices + + bboxes, proj_offset = tgt_crs._bbox_and_offset(src_crs) + mod = np.diff(src_crs.x_limits)[0] + x_lim = xs.min(), xs.max() + for poly in bboxes: + # Arbitrarily choose the number of moduli to look + # above and below the -180->180 range. If data is beyond + # this range, we're not going to transform it quickly. + for i in [-1, 0, 1, 2]: + offset = mod * i - proj_offset + if ((poly[0] + offset) <= x_lim[0] and + (poly[1] + offset) >= x_lim[1]): + vertices = vertices + [[-offset, 0]] + break + return vertices