ENH: Tile warping (#48)

[darribas]

Initial commit to add tile warping functionality

[coveralls]

Coverage increased (+1.4%) to 92.652% when pulling 8877dd4 on tile_warping#48 into 09db67e on master.

[darribas]

Anybody recognises the errors popping up? Never seen them and they seem to do with pytest, not contextily? @jorisvandenbossche @ljwolf

[ljwolf]

Yeah, that's really odd, and it's only affecting the 3.6 runs? never seen this before, but it's probably not anything you've changed.

[jorisvandenbossche]

The error you're seeing was a temporary problem with pytest-cov (the latest pytest release yesterday broke that package, but they already released an update), so it you restart the travis build (or push again), it should work again.

[jorisvandenbossche]

Quick feedback: I tried this branch and I think there is still some error in the tile url determination, as it is downloading the wrong tiles:

df = geopandas.read_file(geopandas.datasets.get_path('nybb'))
ax = df.plot()
contextily.add_basemap(ax, crs=df.crs)

gives me

(instead of New York boroughs)

[jorisvandenbossche]

Something like

--- a/contextily/plotting.py
+++ b/contextily/plotting.py
@@ -84,6 +84,12 @@ def add_basemap(ax, zoom=ZOOM, url=sources.ST_TERRAIN,
         # Extent
         left, right = ax.get_xlim()
         bottom, top = ax.get_ylim()
+        if crs is not None:
+            import pyproj
+            proj_in = pyproj.Proj(crs, preserve_units=True)
+            proj_out = pyproj.Proj({'init': 'epsg:3857'}, preserve_units=True)
+            left, bottom = pyproj.transform(proj_in, proj_out, left, bottom)
+            right, top = pyproj.transform(proj_in, proj_out, right, top)
         # Zoom
         if isinstance(zoom, str) and (zoom.lower() == 'auto'):
             min_ll = _sm2ll(left, bottom)

fixes it. But the question is maybe a bit which library we want to use here exactly.

[darribas]

OK I've worked on this a bit more and I think it actually now fully works as one would expect it. Things I've included:

• User methods in tile.py to warp tiles (warp_tiles) and rasterio-like rasters (warp_img_transform)
• Internal method for the warping (_warper) in tile.py
• Added functionality to add_basemap to warp both tiles pulled from the web and local files
• A demonstration notebook (warping_guide.ipynb) to show you how this works

One thing I wanted to check with you. As @jorisvandenbossche mentions, for this to work, we need reprojection capabilities for points. Currently, I'm doing this in _reproj_bb and I use geopandas. This is probably overkill to reproject two points but I wanted a quick prototype. Ideally, we'd use whatever rasterio uses and we don't need to add any extra dependencies. Do you know what this would be? It should be an easy swap as both the import and the functionality are encapsulated inside _reproj_bb.

If you OK this (@jorisvandenbossche and/or @ljwolf ), I'll add tests and make sure they pass, it should be pretty straightforward.

[ljwolf]

Personally, I think it's fine to use geopandas's crs stuff for this. It's not a particularly onerous dependency given we will already have rasterio as one for the resampling.

[jorisvandenbossche]

Some user feedback:

import geopandas
import contextily

countries = geopandas.read_file(geopandas.datasets.get_path('naturalearth_lowres'))
countries = countries[countries['name'] != 'Antarctica'].copy()

Using WGS84 lon/lat, I need to limit the axis extent manually, otherwise pyproj cannot project the image bounds to WebMercator.

ax = countries.plot(facecolor='none', edgecolor='k')
ax.axis((-180, 180, -60, 84))
contextily.add_basemap(ax, attribution=None)

Should we try to automatically limit the bounds? (we know the bounds of web mercator)
Or, we should actually check if the providers also provide tiles outside of those bounds. I actually suppose they do, which means we should do the conversion of lon/lat to web mercator ourselves (or specify to pyproj to allow out of bounds conversion)

Converting the countries first to WebMercator gives the following bug:

countries2 = countries.to_crs(epsg=3857)

ax = countries2.plot(facecolor='none', edgecolor='k')
contextily.add_basemap(ax, attribution=None, crs=countries2.crs)

(which doesn't happen on master)

In general, I am also not sure we should by default warp the web mercator images to WGS84 / plate carree coordinates. This will give quite distorted images.
The better options would be do automatically convert the axis data to web mercator coordinates under the hood, and keeping the images in their original state (I think this is what all web tile + lon/lat (geojson) based plotting apps do). But, that might not be that easy to do .. (as you probably need to set up such transform when creating the ax ..)

[darribas]

New commits pass locally and implement warping with a default to no warping.

[darribas]

@jorisvandenbossche / @ljwolf, any idea of the error that's returning now? I've not seen it before and it's not the rasterio install from the weekend. Maybe something related to #44?

[jorisvandenbossche]

@darribas seems something went wrong with your rebase / merge of master, as the diff now shows many things of other PRs. Maybe try to squash it into a single commit to clean-up?

[jorisvandenbossche]

any idea of the error that's returning now?

The first one not directly, the second one if from a missing geopandas

[darribas]

How can I squash? I think it was my merge with master, which was effective but not elegant? :-)

[darribas]

hey @jorisvandenbossche this is passing after the last commit to add geopandas as dependency. How should I do the squashing? Do you think it's compulsory? 😄

[darribas]

Re-up of this one for @jorisvandenbossche and @ljwolf to have a look. if It works we can merge in and we're closer to the release by the end of July :)

[jorisvandenbossche]

@darribas I rebased / squashed everything to fix the conflicts, and forced pushed. The diff looks better now (and I hope I fixed everything correctly, if you can quickly check the diff?). Best be careful how to pull it again locally for you.

[jorisvandenbossche]

Some feedback:

• the notebook is not working for me (might be that it needs to be updated for the latest changes, eg that 4326 is no longer the default crs?).
  In addition, I would personally also not show it with EPSG:4326, but rather with a projected CRS (I think that is a better use case to show)
• I still see the bug with the black background I showed above (https://github.com/darribas/contextily/pull/49#issuecomment-490089359). Now it less urgent, as the default is again Mercator, so you don't need to specify the Mercator CRS. But if you do, you get thus the above bug, and it might be pointing to a problem in the reprojection, I don't know.

[jorisvandenbossche]

I think we can rather simply implement this with just pyproj (or even simpler, also no need for creating points of it), and so not need to add geopandas as a dependency (only pyproj). I think the following should be a drop in replacement:

def _reproj_bb(left, right, bottom, top, s_crs, t_crs):
    proj_source = pyproj.Proj(s_crs, preserve_units=True)
    proj_target = pyproj.Proj(t_crs, preserve_units=True)
    x, y = pyproj.transform(proj_source, proj_target, [left, right], [bottom, top])
    left, right = x
    bottom, top = y
    return left, right, bottom, top

[jorisvandenbossche]

Actually, rasterio also has a rasterio.warp.transform_bounds(src_crs, dst_crs, left, bottom, right, top) which I think does exactly what we need here.

[darribas]

Good catch! I don't think that was there last time I checked on rasterio! Addressed in b0354e6

[jorisvandenbossche]

ext -> extent ?

(I find that more readable, and it's also called like that in the return section of the bounds2img docstring)

[darribas]

Good catch! Fixed in 8877dd4

[jorisvandenbossche]

Merged!