diff --git a/examples/user_guide/Working_with_Bokeh.ipynb b/examples/user_guide/Working_with_Bokeh.ipynb index 100bd325..b7e08509 100644 --- a/examples/user_guide/Working_with_Bokeh.ipynb +++ b/examples/user_guide/Working_with_Bokeh.ipynb @@ -75,6 +75,26 @@ "gts.EsriImagery.opts(width=600, height=570, global_extent=True) * gts.StamenLabels.options(level='annotation')" ] }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Note that since the frontend tile renderer only support Mercator coordinates we have to download the tiles and let Cartopy handle reprojecting them in Python. This allows us view them in any coordinate system:" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "import cartopy.crs as ccrs\n", + "\n", + "tiles = gv.tile_sources.ESRI()\n", + "\n", + "gv.util.get_tile_rgb(tiles, bbox=(-180, -70, 180, 70), zoom_level=1).opts(width=500, height=500, projection=ccrs.Orthographic())" + ] + }, { "cell_type": "markdown", "metadata": {}, diff --git a/geoviews/element/geo.py b/geoviews/element/geo.py index 5d928a53..9cad6d8a 100644 --- a/geoviews/element/geo.py +++ b/geoviews/element/geo.py @@ -214,6 +214,11 @@ class WMTS(_GeoFeature): https://maps.wikimedia.org/osm-intl/{Z}/{X}/{Y}@2x.png """ + crs = param.ClassSelector(default=ccrs.GOOGLE_MERCATOR, class_=ccrs.CRS, doc=""" + Cartopy coordinate-reference-system specifying the + coordinate system of the data. Inferred automatically + when _Element wraps cartopy Feature object.""") + group = param.String(default='WMTS') layer = param.String(doc="The layer on the tile service") diff --git a/geoviews/util.py b/geoviews/util.py index 4bff553e..bfd7f67a 100644 --- a/geoviews/util.py +++ b/geoviews/util.py @@ -8,11 +8,13 @@ import shapely.geometry as sgeom from cartopy import crs as ccrs +from cartopy.io.img_tiles import GoogleTiles, QuadtreeTiles +from holoviews.element import Tiles from holoviews.core.util import basestring from shapely.geometry.base import BaseMultipartGeometry from shapely.geometry import ( MultiLineString, LineString, MultiPolygon, Polygon, LinearRing, - Point, MultiPoint + Point, MultiPoint, box ) geom_types = (MultiLineString, LineString, MultiPolygon, Polygon, @@ -106,6 +108,34 @@ def project_extents(extents, src_proj, dest_proj, tol=1e-6): return geom_in_crs.bounds +def zoom_level(bounds, width, height): + """ + Compute zoom level given bounds and the plot size. + """ + w, s, e, n = bounds + max_width, max_height = 256, 256 + ZOOM_MAX = 21 + ln2 = np.log(2) + + def latRad(lat): + sin = np.sin(lat * np.pi / 180) + radX2 = np.log((1 + sin) / (1 - sin)) / 2 + return np.max([np.min([radX2, np.pi]), -np.pi]) / 2 + + def zoom(mapPx, worldPx, fraction): + return np.floor(np.log(mapPx / worldPx / fraction) / ln2) + + latFraction = (latRad(n) - latRad(s)) / np.pi + + lngDiff = e - w + lngFraction = ((lngDiff + 360) if lngDiff < 0 else lngDiff) / 360 + + latZoom = zoom(height, max_height, latFraction) + lngZoom = zoom(width, max_width, lngFraction) + zoom = np.min([latZoom, lngZoom, ZOOM_MAX]) + return int(zoom) if np.isfinite(zoom) else 0 + + def geom_dict_to_array_dict(geom_dict, coord_names=['Longitude', 'Latitude']): """ Converts a dictionary containing an geometry key to a dictionary @@ -660,3 +690,47 @@ def from_xarray(da, crs=None, apply_transform=False, nan_nodata=False, **kwargs) if hasattr(el.data, 'attrs'): el.data.attrs = da.attrs return el + + +def get_tile_rgb(tile_source, bbox, zoom_level, bbox_crs=ccrs.PlateCarree()): + """ + Returns an RGB element given a tile_source, bounding box and zoom level. + + Parameters + ---------- + tile_source: WMTS element or string URL + The tile source to download the tiles from. + bbox: tuple + A four tuple specifying the (left, bottom, right, top) corners of the + domain to download the tiles for. + zoom_level: int + The zoom level at which to download the tiles + bbox_crs: ccrs.CRs + cartopy CRS defining the coordinate system of the supplied bbox + + Returns + ------- + RGB element containing the tile data in the specified bbox + """ + + from .element import RGB, WMTS + if isinstance(tile_source, (WMTS, Tiles)): + tile_source = tile_source.data + + if bbox_crs is not ccrs.GOOGLE_MERCATOR: + bbox = project_extents(bbox, bbox_crs, ccrs.GOOGLE_MERCATOR) + + if '{Q}' in tile_source: + tile_source = QuadtreeTiles(url=tile_source.replace('{Q}', '{tile}')) + else: + tile_source = GoogleTiles(url=tile_source) + + bounds = box(*bbox) + rgb, extent, orient = tile_source.image_for_domain(bounds, zoom_level) + if orient == 'lower': + rgb = rgb[::-1] + x0, x1, y0, y1 = extent + l, b, r, t = bbox + return RGB( + rgb, bounds=(x0, y0, x1, y1), crs=ccrs.GOOGLE_MERCATOR, vdims=['R', 'G', 'B'], + ).clone(datatype=['grid', 'xarray', 'iris'])[l:r, b:t]