From 242d0c54c10ab6e051f480e5fceac9025135f63f Mon Sep 17 00:00:00 2001
From: Philipp Rudiger <prudiger@anaconda.com>
Date: Sat, 13 Mar 2021 13:40:41 +0100
Subject: [PATCH] Add utility to download tile RGB (#458)

---
 examples/user_guide/Working_with_Bokeh.ipynb | 20 ++++++
 geoviews/element/geo.py                      |  5 ++
 geoviews/util.py                             | 76 +++++++++++++++++++-
 3 files changed, 100 insertions(+), 1 deletion(-)

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]