non standard projections for gv.WMTS tiles

[scottyhq]

Is your feature request related to a problem? Please describe.

I'd like to be able to utilize WMTS tiles in a projection other Web Mercator (EPSG:3857), which currently does not seem possible. This is especially useful for Polar regions, where Polar Stereographic projections are commonly used. This functionality is available via ipyleaflet (https://github.com/jupyter-widgets/ipyleaflet/blob/master/examples/CustomProjections.ipynb), so I'm assuming the same should be possible with geoviews...

Describe the solution you'd like

I'd like to either be able to point to a pre-tiled non-EPSG:3857 WMTS such as https://wiki.earthdata.nasa.gov/display/GIBS/GIBS+API+for+Developers#GIBSAPIforDevelopers-Projections&Resolution , or have MapTiles reprojected on the fly.

import geoviews as gv
import geopandas as gpd
from shapely import geometry
import cartopy.crs as ccrs
import hvplot.pandas

S,N,W,E = [74.5, 80, -102, -98]
bbox = geometry.box(W,S,E,N)
gf = gpd.GeoDataFrame(geometry=[bbox], crs=4326)

crs = ccrs.NorthPolarStereo(central_longitude=-45, true_scale_latitude=70)
print(crs.proj4_init)
print(gf.crs)

# Tiles in EPSG3413 ?
tiles = gv.WMTS('https://gibs.earthdata.nasa.gov/wmts/epsg3413/best/BlueMarble_NextGeneration/default/EPSG3413_500m/{Z}/{Y}/{X}.jpeg',
                crs=crs)

box = gf.hvplot.polygons(line_color='m', fill_color='m', fill_alpha=0.1, 
                         ylim=(70,90), xlim=(-180, 180),
                         coastline=True, projection=crs)
tiles * box

Describe alternatives you've considered

Using tiles=True keyword also produces an incorrect map without any tracebacks or warnings, so perhaps this should also be considered a bug and an error or warning is needed. Or maybe reprojecting on the fly is possible and I'm just unable to figure out the right incantation of settings?

Additional context

related issues: #427, #357

[philippjfr]

I would love to have this, and I'll have to look at how ipyleaflet does it. I suspect leaflet.js itself has support for reprojecting tile sources. The main difficulty is finding a way to get this into Bokeh, however in theory we could simply bundle it into GeoViews.js since GeoViews now ships its own JS extension. So we'd have to subclass the tile source implementation from bokeh and add reprojection in there somewhere.

Using tiles=True keyword also produces an incorrect map without any tracebacks or warnings, so perhaps this should also be considered a bug and an error or warning is needed.

A warning should definitely be generated here for the time being.

[jlstevens]

This would really great to have in the next release (personally, I would really prefer it if we could get this functionality into Bokeh and @tonyfast is wondering if this could be done via numba e.g xarray-spatial) . I'll assign this feature request to 1.9.0 but we can bump it if it is too ambitious.

[philippjfr]

One trick now that I've merged the gv.get_tile_rgb utility would be to have a dynamicmap which fetches the tiles server-side and then projects them:

stream = hv.streams.RangeXY()
size = hv.streams.PlotSize(width=400, height=400)
def tile_xy(x_range, y_range, width, height, scale):
    x_range = x_range or ccrs.GOOGLE_MERCATOR.x_limits
    y_range = y_range or ccrs.GOOGLE_MERCATOR.y_limits
    x0, x1 = x_range
    y0, y1 = y_range
    bounds = (x0, y0, x1, y1)
    lat_lon = gv.util.project_extents(bounds, ccrs.GOOGLE_MERCATOR, ccrs.PlateCarree())
    zoom = gv.util.zoom_level(lat_lon, width, height)
    return gv.util.get_tile_rgb(gv.tile_sources.ESRI(), bounds, zoom, ccrs.GOOGLE_MERCATOR)
hv.DynamicMap(tile_xy, streams=[stream, size]).opts(projection=ccrs.Orthographic(), width=500, height=500)

[cwerner]

Can I ask a question about the current state of on-the-fly tile reprojection? Is this currently possible (I'm not super familiar with holoviz/ geoviews so maybe I am missing something)?

I have a large stack of projected xarray data sources (satellite data) that I want to analyse/ select on (and in the end datashade and plot with geoviews together with StamenLabels tiles). Currently, I reproject my xarray data with rioxarray early in the pipeline - however ideally I'd keep my data as is for as long as possible and only reproject it when the plot is generated or even better have geoviews reproject the tile data automatically when it detects it needs to.

Cheers,
C

[jbednar]

I'm a bit confused by that description. This issue is about reprojecting tiles (not currently supported by geoviews), not reprojecting xarray data sources (always supported by geoviews). Xarray data sources will be reprojected as needed if they are gv.Image objects; gv.Image is the same as hv.Image from holoviews but essentially adds that one capability of dynamic reprojection. So you can already keep your data in your preferred projection, if it's not a tile source; projecting early saves time for repeated viewing, but is not required. What you can't yet do without the code in this issue is to leave your data in your preferred projection and overlay it in that same projection on a tile source, because without this code we have no way to reproject tile sources from Web Mercator into your preferred projection. Hope that makes sense?

[cwerner]

Thanks for the explanation and sorry for any confusion. Yes, I’d like to keep my datas’ projection and have a set of tiles projected to the projection of my data.

[jbednar]

Ok, then that specific bit is not supported; instead if you want to delay reprojection, you'll currently need the final plotting projection to be Web Mercator. Hopefully we can get to this issue sometime soon and remove that limitation!

[zxdawn]

I have tried the similar method mentioned by @philippjfr

stream = hv.streams.RangeXY()
size = hv.streams.PlotSize(width=400, height=400)
def tile_xy(x_range, y_range, width, height, scale):
    x_range = x_range or ccrs.GOOGLE_MERCATOR.x_limits
    y_range = y_range or ccrs.GOOGLE_MERCATOR.y_limits
    x0, x1 = x_range
    y0, y1 = y_range
    bounds = (x0, y0, x1, y1)
    lat_lon = gv.util.project_extents(bounds, ccrs.GOOGLE_MERCATOR, ccrs.PlateCarree())
    zoom = gv.util.zoom_level(lat_lon, width, height)
    return gv.util.get_tile_rgb(gv.tile_sources.OSM(), bounds, zoom, ccrs.GOOGLE_MERCATOR)
hv.DynamicMap(tile_xy, streams=[stream, size]).opts(projection=ccrs.NorthPolarStereo(), width=500, height=500)

The resolution looks quite low when I enlarge:

Here's the result of gv.tile_sources.OSM():

[prusswan]

I'd like to either be able to point to a pre-tiled non-EPSG:3857 WMTS such as https://wiki.earthdata.nasa.gov/display/GIBS/GIBS+API+for+Developers#GIBSAPIforDevelopers-Projections&Resolution

Without this, the WMTS implementation is only partial. Aside from the supported projection, in a typical WMTSCapabilities.xml, it is actually possible to specify different tile origins and tile grids/resolutions for each of the zoom levels (usually they share the same tile origin and the resolutions differ by powers of 2, but this is not always the case). To my knowledge, out of the common web mapping libraries, only OpenLayers fully support this. So the current implementation in Geoviews may be better understood as the common XYZ (with EPSG:3857 web mercator projection assumed).

Edit: related #568