Merge pull request #50 from jorisvandenbossche/merge-tiles

Re-implementation of tile merging based on mercantile

Author: darribas

README.md

@@ -21,10 +21,10 @@ the internet into geospatial raster files. Bounding boxes can be passed in both
 
 ## Dependencies
 
-* `cartopy`
 * `mercantile`
 * `numpy`
 * `pandas`
+- `matplotlib`
 * `pillow`
 * `rasterio`
 * `requests`

contextily/tile.py

@@ -9,7 +9,6 @@
 import pandas as pd
 import rasterio as rio
 from PIL import Image
-from cartopy.io.img_tiles import _merge_tiles as merge_tiles
 from rasterio.transform import from_origin
 from . import tile_providers as sources
 
@@ -146,6 +145,7 @@ def bounds2img(w, s, e, n, zoom='auto',
     if zoom == 'auto':
         zoom = _calculate_zoom(w, e, s, n)
     tiles = []
+    arrays = []
     for t in mt.tiles(w, s, e, n, [zoom]):
         x, y, z = t.x, t.y, t.z
         tile_url = url.replace('tileX', str(x)).replace('tileY', str(y)).replace('tileZ', str(z))
@@ -155,17 +155,15 @@ def bounds2img(w, s, e, n, zoom='auto',
             image = Image.open(image_stream).convert('RGB')
             image = np.asarray(image)
         # ---
-        wt, st, et, nt = mt.bounds(t)
-        xr = np.linspace(wt, et, image.shape[0])
-        yr = np.linspace(st, nt, image.shape[1])
-        tiles.append([image, xr, yr, 'lower'])
-    merged, extent = merge_tiles(tiles)[:2]
+        tiles.append(t)
+        arrays.append(image)
+    merged, extent = _merge_tiles(tiles, arrays)
     # lon/lat extent --> Spheric Mercator
-    minX, maxX, minY, maxY = extent
-    w, s = mt.xy(minX, minY)
-    e, n = mt.xy(maxX, maxY)
-    extent = w, e, s, n
-    return merged[::-1], extent
+    west, south, east, north = extent
+    left, bottom = mt.xy(west, south)
+    right, top = mt.xy(east, north)
+    extent = left, right, bottom, top
+    return merged, extent
 
 
 def _retryer(tile_url, wait, max_retries):
@@ -332,3 +330,50 @@ def _calculate_zoom(w, s, e, n):
     zoom_lat = np.ceil(np.log2(360 * 2. / lat_length))
     zoom = np.max([zoom_lon, zoom_lat])
     return int(zoom)
+
+
+def _merge_tiles(tiles, arrays):
+    """
+    Merge a set of tiles into a single array.
+
+    Parameters
+    ---------
+    tiles  : list of mercantile.Tile objects
+             The tiles to merge.
+    arrays : list of numpy arrays
+             The corresponding arrays (image pixels) of the tiles. This list
+             has the same length and order as the `tiles` argument.
+
+    Returns
+    -------
+    img : np.ndarray
+        Merged arrays.
+    extent : tuple
+         Bounding box [west, south, east, north] of the returned image
+         in long/lat.
+    """
+    # create (n_tiles x 2) array with column for x and y coordinates
+    tile_xys = np.array([(t.x, t.y) for t in tiles])
+
+    # get indices starting at zero
+    indices = tile_xys - tile_xys.min(axis=0)
+
+    # the shape of individual tile images
+    h, w, d = arrays[0].shape
+
+    # number of rows and columns in the merged tile
+    n_x, n_y = (indices+1).max(axis=0)
+
+    # empty merged tiles array to be filled in
+    img = np.zeros((h * n_y, w * n_x, d), dtype=np.uint8)
+
+    for ind, arr in zip(indices, arrays):
+        x, y = ind
+        img[y*h:(y+1)*h, x*w:(x+1)*w, :] = arr
+
+    bounds = np.array([mt.bounds(t) for t in tiles])
+    west, south, east, north = (
+        min(bounds[:, 0]), min(bounds[:, 1]),
+        max(bounds[:, 2]), max(bounds[:, 3]))
+
+    return img, (west, south, east, north)

requirements.txt

@@ -1,5 +1,5 @@
-cartopy
 geopy
+matplotlib
 mercantile
 pandas
 pillow

tests/test_ctx.py

@@ -43,10 +43,8 @@ def test_bounds2raster():
     assert_array_almost_equal(rimg.mean(), img.mean())
     assert_array_almost_equal(ext, (0.0, 939258.2035682457, 
                                     6261721.35712164, 6887893.492833804))
-    assert_array_almost_equal(ext, (0.0, 939258.2035682457, 
-                                    6261721.35712164, 6887893.492833804))
-    rtr_bounds = [-613.0928221724841, 6262334.050013727,
-                  938645.1107460733, 6888506.185725891]
+    rtr_bounds = [-611.49622628141, 6262332.853347922,
+                  938646.7073419644, 6888504.989060086]
     assert_array_almost_equal(list(rtr.bounds), rtr_bounds)
 
 def test_bounds2img():
@@ -159,10 +157,10 @@ def test_add_basemap():
     ax_extent = (-11740727.544603072, -11662456.027639052,
                   4852834.0517692715, 4891969.810251278)
     assert_array_almost_equal(ax_extent, ax.images[0].get_extent())
-    assert ax.images[0].get_array().sum() == 75687792
-    assert ax.images[0].get_array().shape == (256, 511, 3)
+    assert ax.images[0].get_array().sum() == 75853866
+    assert ax.images[0].get_array().shape == (256, 512, 3)
     assert_array_almost_equal(ax.images[0].get_array().mean(),
-                              192.86068982387476)
+                              192.90635681152344)
 
     # Test local source
     f, ax = matplotlib.pyplot.subplots(1)
@@ -188,10 +186,10 @@ def test_add_basemap():
     ax_extent = (-11740727.544603072, -11691807.846500559,
                   4852834.0517692715, 4891969.810251278)
     assert_array_almost_equal(ax_extent, ax.images[0].get_extent())
-    assert ax.images[0].get_array().sum() == 719543527
-    assert ax.images[0].get_array().shape == (1021, 1276, 3)
+    assert ax.images[0].get_array().sum() == 723918764
+    assert ax.images[0].get_array().shape == (1024, 1280, 3)
     assert_array_almost_equal(ax.images[0].get_array().mean(),
-                              184.10237852536648)
+                              184.10206197102863)
 
 def test_attribution():
     f, ax = matplotlib.pyplot.subplots(1)