Python mercantile.Tile() Examples

def __init__(self, url, zoom=18, bounds=None):
        self.zoom_level = zoom
        self._name = "image-{}".format(str(uuid.uuid4()))
        self._url = url

        _first_tile = mercantile.Tile(z=self.zoom_level, x=0, y=0)
        _last_tile = mercantile.Tile(z=self.zoom_level, x=180, y=-85.05)
        g = box(*mercantile.xy_bounds(_first_tile)).union(box(*mercantile.xy_bounds(_last_tile)))

        self._full_bounds = g.bounds

        # TODO: populate rest of fields automatically
        self._tile_size = 256
        self._nbands = 3
        self._dtype = "uint8"
        self.bounds = self._expand_bounds(bounds) 

def get_assets_from_json(
    tiles: Dict, quadkey_zoom: int, x: int, y: int, z: int
) -> List[str]:
    """Find assets."""
    mercator_tile = mercantile.Tile(x=x, y=y, z=z)
    quadkeys = find_quadkeys(mercator_tile, quadkey_zoom)

    assets = list(itertools.chain.from_iterable([tiles.get(qk, []) for qk in quadkeys]))

    # check if we have a mosaic in the url (.json/.gz)
    return list(
        itertools.chain.from_iterable(
            [
                get_assets_from_json(tiles, quadkey_zoom, x, y, z)
                if os.path.splitext(asset)[1] in [".json", ".gz"]
                else [asset]
                for asset in assets
            ]
        )
    ) 

def test_tile_read_extmask():
    """Read masked area."""
    # non-boundless tile covering the masked part
    mercator_tile = mercantile.Tile(x=876431, y=1603669, z=22)
    bounds = mercantile.xy_bounds(mercator_tile)
    with rasterio.Env(GDAL_DISABLE_READDIR_ON_OPEN="TRUE"):
        with rasterio.open(S3_EXTMASK_PATH) as src_dst:
            arr, mask = reader.part(src_dst, bounds, 256, 256)
        assert arr.shape == (3, 256, 256)
        assert mask.shape == (256, 256)
        assert not mask.all()

    # boundless tile covering the masked part
    mercator_tile = mercantile.Tile(x=876431, y=1603668, z=22)
    bounds = mercantile.xy_bounds(mercator_tile)
    with rasterio.Env(GDAL_DISABLE_READDIR_ON_OPEN="EMPTY_DIR"):
        with rasterio.open(S3_MASK_PATH) as src_dst:
            arr, mask = reader.part(src_dst, bounds, 256, 256)
        assert arr.shape == (3, 256, 256)
        assert not mask.all() 

def tile(z, x, y):

    # Todo: predictor should take care of zoom levels
    if z != 18:
        abort(404)

    tile = mercantile.Tile(x, y, z)

    url = tiles.format(x=tile.x, y=tile.y, z=tile.z)
    res = fetch_image(session, url)

    if not res:
        abort(500)

    image = Image.open(res)

    mask = predictor.segment(image)

    return send_png(mask) 

def tiles_from_csv(path):
    """Read tiles from a line-delimited csv file.

    Args:
      file: the path to read the csv file from.

    Yields:
      The mercantile tiles from the csv file.
    """

    with open(path) as fp:
        reader = csv.reader(fp)

        for row in reader:
            if not row:
                continue

            yield mercantile.Tile(*map(int, row)) 

def adjacent_tile(tile, dx, dy, tiles):
    """Retrieves an adjacent tile from a tile store.

    Args:
      tile: the original tile to get an adjacent tile for.
      dx: the offset in tile x direction.
      dy: the offset in tile y direction.
      tiles: the tile store to get tiles from; must support `__getitem__` with tiles.

    Returns:
      The adjacent tile's image or `None` if it does not exist.
    """

    x, y, z = map(int, [tile.x, tile.y, tile.z])
    other = mercantile.Tile(x=x + dx, y=y + dy, z=z)

    try:
        path = tiles[other]
        return Image.open(path).convert("RGB")
    except KeyError:
        return None 

def get_tile_data(driver: Driver,
                  keys: Union[Sequence[str], Mapping[str, str]],
                  tile_xyz: Tuple[int, int, int] = None,
                  *, tile_size: Tuple[int, int] = (256, 256),
                  preserve_values: bool = False,
                  asynchronous: bool = False) -> Any:
    """Retrieve raster image from driver for given XYZ tile and keys"""

    if tile_xyz is None:
        # read whole dataset
        return driver.get_raster_tile(
            keys, tile_size=tile_size, preserve_values=preserve_values,
            asynchronous=asynchronous
        )

    # determine bounds for given tile
    metadata = driver.get_metadata(keys)
    wgs_bounds = metadata['bounds']

    tile_x, tile_y, tile_z = tile_xyz

    if not tile_exists(wgs_bounds, tile_x, tile_y, tile_z):
        raise exceptions.TileOutOfBoundsError(
            f'Tile {tile_z}/{tile_x}/{tile_y} is outside image bounds'
        )

    mercator_tile = mercantile.Tile(x=tile_x, y=tile_y, z=tile_z)
    target_bounds = mercantile.xy_bounds(mercator_tile)

    return driver.get_raster_tile(
        keys, tile_bounds=target_bounds, tile_size=tile_size,
        preserve_values=preserve_values, asynchronous=asynchronous
    ) 

def get_assets(self, x: int, y: int, z: int) -> List[str]:
        """Find assets."""
        mercator_tile = mercantile.Tile(x=x, y=y, z=z)
        quadkeys = find_quadkeys(mercator_tile, self.quadkey_zoom)

        assets = list(
            itertools.chain.from_iterable(
                [self._fetch_dynamodb(qk).get("assets", []) for qk in quadkeys]
            )
        )

        # Find mosaics recursively?
        return assets 

def tiles_to_bounds(tiles: List[mercantile.Tile]) -> List[float]:
    """Get bounds from a set of mercator tiles."""
    zoom = tiles[0].z
    xyz = numpy.array([[t.x, t.y, t.z] for t in tiles])
    extrema = {
        "x": {"min": xyz[:, 0].min(), "max": xyz[:, 0].max() + 1},
        "y": {"min": xyz[:, 1].min(), "max": xyz[:, 1].max() + 1},
    }
    ulx, uly = mercantile.ul(extrema["x"]["min"], extrema["y"]["min"], zoom)
    lrx, lry = mercantile.ul(extrema["x"]["max"], extrema["y"]["max"], zoom)
    return [ulx, lry, lrx, uly] 

def default_filter(
    tile: mercantile.Tile,
    dataset: Sequence[Dict],
    geoms: Sequence[polygons],
    minimum_tile_cover=None,
    tile_cover_sort=False,
    maximum_items_per_tile: Optional[int] = None,
) -> List:
    """Filter and/or sort dataset per intersection coverage."""
    indices = list(range(len(dataset)))

    if minimum_tile_cover or tile_cover_sort:
        tile_geom = polygons(mercantile.feature(tile)["geometry"]["coordinates"][0])
        int_pcts = _intersect_percent(tile_geom, geoms)

        if minimum_tile_cover:
            indices = [ind for ind in indices if int_pcts[ind] > minimum_tile_cover]

        if tile_cover_sort:
            # https://stackoverflow.com/a/9764364
            indices, _ = zip(*sorted(zip(indices, int_pcts), reverse=True))

    if maximum_items_per_tile:
        indices = indices[:maximum_items_per_tile]

    return [dataset[ind] for ind in indices] 

def test_tiles_to_bounds():
    """Get tiles bounds for zoom level."""
    tiles = [mercantile.Tile(x=150, y=182, z=9), mercantile.Tile(x=151, y=182, z=9)]
    assert len(utils.tiles_to_bounds(tiles)) == 4 

def test_find_quadkeys():
    """Get correct quadkeys"""
    tile = mercantile.Tile(150, 182, 9)
    assert utils.find_quadkeys(tile, 8) == ["03023033"]
    assert utils.find_quadkeys(tile, 9) == ["030230330"]
    assert utils.find_quadkeys(tile, 10) == [
        "0302303300",
        "0302303301",
        "0302303303",
        "0302303302",
    ] 

def read_tile(self, z, x, y):
        """Read raster tile data and mask."""
        mercator_tile = mercantile.Tile(x=x, y=y, z=z)
        tile_bounds = mercantile.xy_bounds(mercator_tile)
        return tile_read(
            self.path,
            tile_bounds,
            self.tiles_size,
            indexes=self.indexes,
            nodata=self.nodata,
        ) 

def read_tile(self, z, x, y):
        """Read raster tile data and mask."""
        mercator_tile = mercantile.Tile(x=x, y=y, z=z)
        tile_bounds = mercantile.xy_bounds(mercator_tile)

        data, mask = tile_read(
            self.path,
            tile_bounds,
            self.tiles_size,
            indexes=self.indexes,
            nodata=self.nodata,
        )
        data = (data[0] + data[1]) / 2
        return data.astype(numpy.uint8), mask 

def _expand_bounds(self, bounds):
        if bounds is None:
            return bounds
        min_tile_x, min_tile_y, max_tile_x, max_tile_y = self._tile_coords(bounds)

        ul = box(*mercantile.xy_bounds(mercantile.Tile(z=self.zoom_level, x=min_tile_x, y=max_tile_y)))
        lr = box(*mercantile.xy_bounds(mercantile.Tile(z=self.zoom_level, x=max_tile_x, y=min_tile_y)))

        return ul.union(lr).bounds 

def tiles_from_slippy_map(root):
    """Loads files from an on-disk slippy map directory structure.

    Args:
      root: the base directory with layout `z/x/y.*`.

    Yields:
      The mercantile tiles and file paths from the slippy map directory.
    """

    # The Python string functions (.isdigit, .isdecimal, etc.) handle
    # unicode codepoints; we only care about digits convertible to int
    def isdigit(v):
        try:
            _ = int(v)  # noqa: F841
            return True
        except ValueError:
            return False

    for z in os.listdir(root):
        if not isdigit(z):
            continue

        for x in os.listdir(os.path.join(root, z)):
            if not isdigit(x):
                continue

            for name in os.listdir(os.path.join(root, z, x)):
                y = os.path.splitext(name)[0]

                if not isdigit(y):
                    continue

                tile = mercantile.Tile(x=int(x), y=int(y), z=int(z))
                path = os.path.join(root, z, x, name)
                yield tile, path 

def test_burn_without_feature(self):
        parking_fc = get_parking()

        # This tile does not have a parking lot in our fixtures.
        tile = mercantile.Tile(69623, 104946, 18)

        rasterized = burn(tile, parking_fc["features"], 512)
        rasterized = Image.fromarray(rasterized, mode="P")

        self.assertEqual(rasterized.size, (512, 512))

        # Tile does not have a parking feature in our fixture, the sum of pixels is zero.
        self.assertEqual(np.sum(rasterized), 0) 

def test_burn_with_feature(self):
        parking_fc = get_parking()

        # The tile below has a parking lot in our fixtures.
        tile = mercantile.Tile(70762, 104119, 18)

        rasterized = burn(tile, parking_fc["features"], 512)
        rasterized = Image.fromarray(rasterized, mode="P")

        # rasterized.save('rasterized.png')

        self.assertEqual(rasterized.size, (512, 512))

        # Tile has a parking feature in our fixtures, thus sum should be non-zero.
        self.assertNotEqual(np.sum(rasterized), 0) 

def test_read_tiles(self):
        filename = "tests/fixtures/tiles.csv"
        tiles = [tile for tile in tiles_from_csv(filename)]

        self.assertEqual(len(tiles), 3)
        self.assertEqual(tiles[0], mercantile.Tile(69623, 104945, 18)) 

def test_slippy_map_directory(self):
        root = "tests/fixtures/images"
        tiles = [tile for tile in tiles_from_slippy_map(root)]
        self.assertEqual(len(tiles), 3)

        tile, path = tiles[0]
        self.assertEqual(type(tile), mercantile.Tile)
        self.assertEqual(path, "tests/fixtures/images/18/69105/105093.jpg") 

def test_getitem(self):
        inputs = ["tests/fixtures/images/"]
        target = "tests/fixtures/labels/"

        transform = JointCompose([JointTransform(ImageToTensor(), MaskToTensor())])
        dataset = SlippyMapTilesConcatenation(inputs, target, transform)

        images, mask, tiles = dataset[0]
        self.assertEqual(tiles[0], mercantile.Tile(69105, 105093, 18))
        self.assertEqual(type(images), torch.Tensor)
        self.assertEqual(type(mask), torch.Tensor) 

def main(args):
    dataset = load_config(args.dataset)

    classes = dataset["common"]["classes"]
    colors = dataset["common"]["colors"]
    assert len(classes) == len(colors), "classes and colors coincide"

    assert len(colors) == 2, "only binary models supported right now"
    bg = colors[0]
    fg = colors[1]

    os.makedirs(args.out, exist_ok=True)

    # We can only rasterize all tiles at a single zoom.
    assert all(tile.z == args.zoom for tile in tiles_from_csv(args.tiles))

    with open(args.features) as f:
        fc = json.load(f)

    # Find all tiles the features cover and make a map object for quick lookup.
    feature_map = collections.defaultdict(list)
    for i, feature in enumerate(tqdm(fc["features"], ascii=True, unit="feature")):

        if feature["geometry"]["type"] != "Polygon":
            continue

        try:
            for tile in burntiles.burn([feature], zoom=args.zoom):
                feature_map[mercantile.Tile(*tile)].append(feature)
        except ValueError as e:
            print("Warning: invalid feature {}, skipping".format(i), file=sys.stderr)
            continue

    # Burn features to tiles and write to a slippy map directory.
    for tile in tqdm(list(tiles_from_csv(args.tiles)), ascii=True, unit="tile"):
        if tile in feature_map:
            out = burn(tile, feature_map[tile], args.size)
        else:
            out = np.zeros(shape=(args.size, args.size), dtype=np.uint8)

        out_dir = os.path.join(args.out, str(tile.z), str(tile.x))
        os.makedirs(out_dir, exist_ok=True)

        out_path = os.path.join(out_dir, "{}.png".format(tile.y))

        if os.path.exists(out_path):
            prev = np.array(Image.open(out_path))
            out = np.maximum(out, prev)

        out = Image.fromarray(out, mode="P")

        palette = make_palette(bg, fg)
        out.putpalette(palette)

        out.save(out_path, optimize=True) 

def process_tile(tile):
    """Process a single MBTiles tile

    Parameters
    ----------
    tile : mercantile.Tile

    Returns
    -------

    tile : mercantile.Tile
        The input tile.
    bytes : bytearray
        Image bytes corresponding to the tile.

    """
    global base_kwds, resampling, src

    # Get the bounds of the tile.
    ulx, uly = mercantile.xy(
        *mercantile.ul(tile.x, tile.y, tile.z))
    lrx, lry = mercantile.xy(
        *mercantile.ul(tile.x + 1, tile.y + 1, tile.z))

    kwds = base_kwds.copy()
    kwds['transform'] = transform_from_bounds(ulx, lry, lrx, uly,
                                              kwds['width'], kwds['height'])
    src_nodata = kwds.pop('src_nodata', None)
    dst_nodata = kwds.pop('dst_nodata', None)

    warnings.simplefilter('ignore')

    with MemoryFile() as memfile:

        with memfile.open(**kwds) as tmp:

            # determine window of source raster corresponding to the tile
            # image, with small buffer at edges
            try:
                west, south, east, north = transform_bounds(TILES_CRS, src.crs, ulx, lry, lrx, uly)
                tile_window = window_from_bounds(west, south, east, north, transform=src.transform)
                adjusted_tile_window = Window(
                    tile_window.col_off - 1, tile_window.row_off - 1,
                    tile_window.width + 2, tile_window.height + 2)
                tile_window = adjusted_tile_window.round_offsets().round_shape()

                # if no data in window, skip processing the tile
                if not src.read_masks(1, window=tile_window).any():
                    return tile, None

            except ValueError:
                log.info("Tile %r will not be skipped, even if empty. This is harmless.", tile)

            reproject(rasterio.band(src, tmp.indexes),
                      rasterio.band(tmp, tmp.indexes),
                      src_nodata=src_nodata,
                      dst_nodata=dst_nodata,
                      num_threads=1,
                      resampling=resampling)

        return tile, memfile.read() 

def test_metadata():
    """Should return correct metadata."""
    with rasterio.open(COG_CMAP) as src_dst:
        meta = reader.metadata(src_dst)
        assert meta["dtype"] == "int8"
        assert meta["colorinterp"] == ["palette"]
        assert not meta.get("scale")
        assert not meta.get("ofsset")
        assert meta.get("colormap")

    with rasterio.open(COG_SCALE) as src_dst:
        meta = reader.metadata(src_dst)
        assert meta["dtype"] == "int16"
        assert meta["colorinterp"] == ["gray"]
        assert meta["scale"] == 0.0001
        assert meta["offset"] == 1000.0
        assert meta["band_descriptions"] == [(1, "Green")]
        assert not meta.get("colormap")
        assert meta["nodata_type"] == "Nodata"

        meta = reader.metadata(src_dst, indexes=1)
        assert meta["colorinterp"] == ["gray"]

        bounds = mercantile.bounds(mercantile.Tile(x=218, y=99, z=8))
        meta = reader.metadata(src_dst, bounds)
        assert meta["colorinterp"] == ["gray"]
        assert meta["bounds"] == bounds

    with rasterio.open(S3_ALPHA_PATH) as src_dst:
        with pytest.warns(AlphaBandWarning):
            meta = reader.metadata(src_dst)
            assert len(meta["band_descriptions"]) == 3
            assert meta["colorinterp"] == ["red", "green", "blue"]
            assert meta["nodata_type"] == "Alpha"

        meta = reader.metadata(src_dst, indexes=(1, 2, 3, 4))
        assert len(meta["band_descriptions"]) == 4
        assert meta["colorinterp"] == ["red", "green", "blue", "alpha"]
        assert meta["nodata_type"] == "Alpha"

    with rasterio.open(S3_MASK_PATH) as src_dst:
        meta = reader.metadata(src_dst)
        assert meta["nodata_type"] == "Mask"