Skip to content

Commit

Permalink
Further fixes for path projection
Browse files Browse the repository at this point in the history
  • Loading branch information
philippjfr committed May 19, 2018
1 parent 8f4e644 commit 2ad395c
Show file tree
Hide file tree
Showing 2 changed files with 39 additions and 7 deletions.
40 changes: 35 additions & 5 deletions geoviews/operation/projection.py
Original file line number Diff line number Diff line change
Expand Up @@ -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,
Expand Down Expand Up @@ -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():
Expand Down Expand Up @@ -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)


Expand Down Expand Up @@ -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)


Expand Down
6 changes: 4 additions & 2 deletions geoviews/util.py
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down Expand Up @@ -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':
Expand Down Expand Up @@ -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.
"""
Expand Down

0 comments on commit 2ad395c

Please sign in to comment.