On-the-fly reprojection of tile maps

[darribas]

One that I've discussed with @jorisvandenbossche offline and now have made some impromptu progress. Currently, contextily only supports Web Mercator so, whatever CRS the original data is expressed in, if the user wants to combine them with tile maps, they need to move their data to WM. In some cases, that's OK but, in other contexts, that might be expensive and it might be a lot better to reproject (warp) the tile map.

To get started, here's a tentative function that does this:

from rasterio.io import MemoryFile
from rasterio.vrt import WarpedVRT
from rasterio.enums import Resampling

def warp_tiles(img, ext, 
               t_crs='EPSG:4326',
               resampling=Resampling.bilinear):
    '''
    Reproject (warp) a Web Mercator basemap into any CRS on-the-fly
    ...
    
    Arguments
    ---------
    img         : ndarray
                  Image as a 3D array of RGB values (e.g. as returned from
                  `contextily.bounds2img`)
    ext         : tuple
                  Bounding box [minX, maxX, minY, maxY] of the returned image
    t_crs       : str/CRS
                  [Optional. Default='EPSG:4326'] Target CRS, expressed in any
                  format permitted by rasterio. Defaults to WGS84 (lon/lat)
    resampling  : <enum 'Resampling'>
                  [Optional. Default=Resampling.bilinear] Resampling method for
                  executing warping, expressed as a `rasterio.enums.Resampling
                  method
    
    Returns
    -------
    img         : ndarray
                  Warped Image as a 3D array of RGB values
    ext         : tuple
                  Bounding box [minX, maxX, minY, maxY] of the returned (warped)
                  image
    '''
    h, w, b = img.shape
    # --- https://mapbox.github.io/rasterio/quickstart.html#opening-a-dataset-in-writing-mode
    minX, maxX, minY, maxY = ext
    x = np.linspace(minX, maxX, w)
    y = np.linspace(minY, maxY, h)
    resX = (x[-1] - x[0]) / w
    resY = (y[-1] - y[0]) / h
    transform = from_origin(x[0] - resX / 2,
                            y[-1] + resY / 2, resX, resY)
    # --- Write basemap into memory file
    with MemoryFile() as memfile:
        with memfile.open(driver='GTiff', height=h, width=w, \
                          count=b, dtype=str(img.dtype.name), \
                          crs='epsg:3857', transform=transform) as mraster:
            for band in range(b):
                mraster.write(img[:, :, band], band+1)
            # --- Virtual Warp
            with WarpedVRT(mraster, crs=t_crs,
                           resampling=resampling) as vrt:
                src_wm = vrt.read().transpose(1, 2, 0)
    bb = vrt.bounds
    extent = bb.left, bb.right, bb.bottom, bb.top
    return src_wm, extent

I'm putting it in this precarious format, rather than as a PR so:

1. Folks quickly say if this API works

2. We can discuss where it could go (my first thought is on tile.py but haven't thought too much about it

3. Folks can also chip in on how this could be exposed to end-users. For now, this is all I can think of:

   • As a standalone function warp_tiles that can be piped with bounds2img
   • Through a t_crs argument in add_basemap that defaults to EPSG:3857 and adds a warped basemap to a set of provided matplotlib.Axes
   • Through the Place API, where the user can also specify a t_crs argument when initialising and the response is returned in that CRS (maybe @choldgraf has some input on this?).

Ideas/suggestions?

[choldgraf]

This is a cool idea -I don't have a ton of perspective on this one (I haven't used the Place API for some time now...maybe @lawasser does?). However why not try making the lower-level functions first and once that stabilizes, think about providing the higher-level API via place.

[ljwolf]

Thanks for doing this @darribas, looks great.

I've not seen the memfile approach before, but that's a clever way to get the transform applied!

I know I'd use the first two suggested functions.

Do you know if this already accepts a geodataframe's .crs attribute? That'd be useful for a call like:

ax = gdf.plot('feature')
contextily.add_basemap(ax, crs=gdf.crs)

So, if a None crs is passed, assume webmercator (since that's a no-op) and otherwise turn geopandas crs dicts into rasterio format.

I can take a look at the diff to make this work, but that'd be easier in a PR 😁

[darribas]

OK, cool, I've just issued #49 with an initial commit along the lines suggested by @ljwolf. It's tagged as WIP but, once we're set this is what we want and how we wanted, I can add tests for the new functionality. So far I've added tile.py/warp_tiles() and included a crs option in plotting.py/add_basemap() that enables dynamic warping in that method.

@ljwolf, I've not played extensively with CRS formats, but I did toy around with what you suggest (I had the same idea) and it seems to work. So I'm not sure every single CRS format is accepted, but {'init': 'EPSG:XXX'} seems to work just fine.

[jorisvandenbossche]

Cool!

When we talked previously about it, I was also thinking about the following related things:

• Next to the default support for WebMercator, it might also be nice to out of the box support long/lat values (and autodetect this case):

  • Reason: it's a very typical use case (eg any geojson file), making it nice to support it out of the box without specific user action
  • Of course, the biggest disadvantage is a small possibility of ambiguity, but, I would say that this is very small (given that for lat/long, the absolute value < 180, while such values for Mercator are rather uncommonly small and located somewhere in the ocean)

• Support for cartopy: detect that an ax is a GeoAxes, and then also automatically do the appropriate thing (passing a certain transform? need to check further what the best action might be here)

[darribas]

Assuming the user had a GeoDataFrame named gdf in EPSG:4326, all they'd have to do is:

contextily.add_basemap(ax, crs=gdf.crs)

I think that's easy enough to not justify a specific case. What do people think?

On cartopy, I think that's worth exploring, but potentially in a separate PR? I don't use cartopy enough to know how to best approach it, but it strikes me as a different feature.

[jorisvandenbossche]

I think that's easy enough to not justify a specific case. What do people think?

It's certainly quite easy to do. But, still making this special case also avoids the to_crs(epsg=3857) in case you have lat/lon data to start with. Being able to avoid is more worth it IMO.

---

Longer reasoning: the main reason that I might still do it is the following. A lot of people have data in lat/lon, and that is the defacto default for geo data on the web, (interactive) web plotting, geojson, etc (and the WebMercator projection handling almost always happens under the hood in similar applications).
So if we would start over with contextily, I would personally make this the default expected input data instead of WebMercator. Because now, you should almost never convert your data to WebMercator (it's not a useful projection for anything else than plotting, eg local area calculation), only just to be able to plot it.

So in an 'ideal' workflow, I would say: if people have lat/lon data, and they don't need to do area/distance based calculations, just use the data as is in lat/lon for visualisation (no need to take care of converting it to WebMercator just for plotting). If people project their lat/lon data to a local projected CRS for other reasons than plotting (or they already have such data), then they also need to remember to pass that CRS to the add_basemap function in order to have correct plotting (which is also nicer than the current need to convert to WebMercator, for sure).

[darribas]

Discussion on development for this feature currently moved over to the PR for this #49

[darribas]

Closing as implemented in #49