diff --git a/README.md b/README.md index 0aeac0eb..1e23fa72 100644 --- a/README.md +++ b/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` diff --git a/contextily/tile.py b/contextily/tile.py index e7265531..ac47c83e 100644 --- a/contextily/tile.py +++ b/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) diff --git a/requirements.txt b/requirements.txt index dd414e6f..478a6781 100644 --- a/requirements.txt +++ b/requirements.txt @@ -1,5 +1,5 @@ -cartopy geopy +matplotlib mercantile pandas pillow diff --git a/tests/test_ctx.py b/tests/test_ctx.py index 1d4e5630..639d30a7 100644 --- a/tests/test_ctx.py +++ b/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)