Skip to content

Commit

Permalink
Merge 1c1e2a3 into 16b7bec
Browse files Browse the repository at this point in the history
  • Loading branch information
philippjfr committed May 22, 2018
2 parents 16b7bec + 1c1e2a3 commit ca2fb77
Show file tree
Hide file tree
Showing 5 changed files with 163 additions and 57 deletions.
1 change: 1 addition & 0 deletions geoviews/data/geopandas.py
Expand Up @@ -103,6 +103,7 @@ def shape(cls, dataset):
@classmethod
def length(cls, dataset):
length = 0
if len(dataset.data) == 0: return 0
arr = geom_to_array(dataset.data.geometry.iloc[0])
ds = dataset.clone(arr, datatype=cls.subtypes, vdims=[])
for d in dataset.data.geometry:
Expand Down
1 change: 1 addition & 0 deletions geoviews/element/geo.py
Expand Up @@ -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))
Expand Down
161 changes: 107 additions & 54 deletions geoviews/operation/projection.py
Expand Up @@ -5,12 +5,12 @@
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,
HexTiles, Labels)
from ..util import project_extents, geom_to_array
from ..util import project_extents, geom_to_array, wrap_path_data, is_multi_geometry


class _project_operation(Operation):
Expand All @@ -31,7 +31,6 @@ def _process(self, element, key=None):
return element.map(self._process_element, self.supported_types)



class project_path(_project_operation):
"""
Projects Polygons and Path Elements from their source coordinate
Expand All @@ -40,54 +39,109 @@ class project_path(_project_operation):

supported_types = [Polygons, Path, Contours, EdgePaths]

def _project_path(self, element, path, data, boundary, geom_type, multi_type):
"""
Handle case of continuously varying path
"""
xdim, ydim = path.kdims[:2]
xs, ys = (path.dimension_values(0) for i in range(2))
if not len(xs):
return []

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 hasattr(element.crs, '_bbox_and_offset'):
vertices = wrap_path_data(vertices, element.crs, element.crs)
path = geom_type(vertices)
path = path.intersection(boundary)
proj = self.p.projection.project_geometry(path, element.crs)
proj_arr = geom_to_array(proj)
data[xdim.name] = proj_arr[:, 0]
data[ydim.name] = proj_arr[:, 1]
return [data]

def _project_contour(self, element, contour, data, boundary, geom_type, multi_type):
"""
Handle case of iso-contour
"""
xdim, ydim = contour.kdims[:2]
data = {k: vals[0] for k, vals in data.items()}

# Wrap longitudes
vertices = contour.array([0, 1])
if hasattr(element.crs, '_bbox_and_offset'):
vertices = wrap_path_data(vertices, element.crs, element.crs)
geom = type(element)([vertices]).geom()

# Clip path to projection boundaries
geoms = []
for g in geom:
if np.isinf(np.array(g.array_interface_base['data'])).sum():
# Skip if infinity in path
continue
try:
# Compute boundary intersections
g = g.intersection(boundary)
except:
continue
if is_multi_geometry(g):
for p in g:
try:
geoms.append(geom_type(p))
except:
continue
else:
geoms.append(g)
if not geoms:
return []

# Project geometry
projected = []
geom = multi_type(geoms) if len(geoms) > 1 else geom
proj = self.p.projection.project_geometry(geom, element.crs)
proj = proj if is_multi_geometry(proj) else [proj]
for geom in proj:
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 projected

def _project_geodataframe(self, element):
geoms = element.split(datatype='geom')
projected = [self.p.projection.project_geometry(geom, element.crs)
for geom in geoms]
new_data = element.data.copy()
new_data['geometry'] = projected
return element.clone(new_data, crs=self.p.projection)

def _process_element(self, element):
if not len(element):
return element.clone(crs=self.p.projection)
elif element.interface.datatype == 'geodataframe':
geoms = element.split(datatype='geom')
projected = [self.p.projection.project_geometry(geom, element.crs)
for geom in geoms]
new_data = element.data.copy()
new_data['geometry'] = projected
return element.clone(new_data, crs=self.p.projection)

boundary_poly = element.crs.project_geometry(Polygon(self.p.projection.boundary),
self.p.projection)

geom_type = Polygon if isinstance(element, Polygons) else LineString
xdim, ydim = element.kdims[:2]
return self._project_geodataframe(element)

boundary = element.crs.project_geometry(Polygon(self.p.projection.boundary),
self.p.projection)

if isinstance(element, Polygons):
multi_type, geom_type = MultiPolygon, Polygon
else:
multi_type, geom_type = MultiLineString, LineString

projected = []
for path in element.split():
data = {vd.name: path.dimension_values(vd, expanded=False) for vd in path.vdims}
if any(len(vals) > 1 for vals in data.values()):
# Handle continuously varying path case
geom = path.columns()
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)
geom[xdim.name] = proj_arr[:, 0]
geom[ydim.name] = proj_arr[:, 1]
projected.append(geom)
projected += self._project_path(element, path, data, boundary, geom_type, multi_type)
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 not geom:
continue
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}))
projected += self._project_contour(element, path, data, boundary, geom_type, multi_type)
return element.clone(projected, crs=self.p.projection)



class project_shape(_project_operation):
"""
Projects Shape Element from the source coordinate reference system
Expand All @@ -99,9 +153,13 @@ class project_shape(_project_operation):
def _process_element(self, element):
if not len(element):
return element.clone(crs=self.p.projection)
geom = self.p.projection.project_geometry(element.geom(), element.crs)
if isinstance(geom, tuple):
geom = [g for g in geom if g != []][0]
geom = element.geom()
vertices = geom_to_array(geom)
if isinstance(geom, (MultiPolygon, Polygon)):
obj = Polygons([vertices])
else:
obj = Path([vertices])
geom = project_path(obj, projection=self.p.projection).geom()
return element.clone(geom, crs=self.p.projection)


Expand Down Expand Up @@ -315,16 +373,11 @@ class project(Operation):
instantiate=False, doc="""
Projection the image type is projected to.""")

_operations = [project_path, project_image, project_shape,
project_graph, project_quadmesh, project_points]

def _process(self, element, key=None):
element = element.map(project_path.instance(projection=self.p.projection),
project_path.supported_types)
element = element.map(project_image.instance(projection=self.p.projection),
project_image.supported_types)
element = element.map(project_shape.instance(projection=self.p.projection),
project_shape.supported_types)
element = element.map(project_graph.instance(projection=self.p.projection),
project_graph.supported_types)
element = element.map(project_quadmesh.instance(projection=self.p.projection),
project_quadmesh.supported_types)
return element.map(project_points.instance(projection=self.p.projection),
project_points.supported_types)
for op in self._operations:
element = element.map(op.instance(projection=self.p.projection),
op.supported_types)
return element
4 changes: 2 additions & 2 deletions geoviews/plotting/bokeh/plot.py
Expand Up @@ -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]
Expand Down Expand Up @@ -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)

Expand Down
53 changes: 52 additions & 1 deletion 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)

Expand Down Expand Up @@ -69,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 @@ -99,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 All @@ -114,9 +118,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:
Expand Down Expand Up @@ -164,3 +178,40 @@ def geo_mesh(element):
zs = np.ma.concatenate([zs, zs[:, 0:1]], axis=1)
return xs, ys, zs


def wrap_path_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


def is_multi_geometry(geom):
"""
Whether the shapely geometry is a Multi or Collection type.
"""
return 'Multi' in geom.geom_type or 'Collection' in geom.geom_type

0 comments on commit ca2fb77

Please sign in to comment.