Python mercantile.xy_bounds() Examples

def burn(tile, features, size):
    """Burn tile with features.

    Args:
      tile: the mercantile tile to burn.
      features: the geojson features to burn.
      size: the size of burned image.

    Returns:
      image: rasterized file of size with features burned.
    """

    # the value you want in the output raster where a shape exists
    burnval = 1
    shapes = ((geometry, burnval) for feature in features for geometry in feature_to_mercator(feature))

    bounds = mercantile.xy_bounds(tile)
    transform = from_bounds(*bounds, size, size)

    return rasterize(shapes, out_shape=(size, size), transform=transform) 

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 __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 geotiff_options(
    x: int,
    y: int,
    z: int,
    tilesize: int = 256,
    dst_crs: CRS = constants.WEB_MERCATOR_CRS,
) -> Dict:
    """
    GeoTIFF options.

    Attributes
    ----------
        x : int
            Mercator tile X index.
        y : int
            Mercator tile Y index.
        z : int
            Mercator tile ZOOM level.
        tilesize : int, optional
            Output tile size. Default is 256.
        dst_crs: CRS, optional
            Target coordinate reference system, default is "epsg:3857".

    Returns
    -------
        dict

    """
    bounds = mercantile.xy_bounds(mercantile.Tile(x=x, y=y, z=z))
    dst_transform = from_bounds(*bounds, tilesize, tilesize)
    return dict(crs=dst_crs, transform=dst_transform) 

def read_tile(src_path, tile):
    """Benchmark rio-tiler.utils._tile_read."""
    tile_bounds = mercantile.xy_bounds(tile)
    # We make sure to not store things in cache.
    with rasterio.Env(GDAL_CACHEMAX=0, NUM_THREADS="all"):
        with rasterio.open(src_path) as src_dst:
            return reader.part(
                src_dst,
                tile_bounds,
                256,
                256,
                resampling_method="nearest",
                dst_crs=constants.WEB_MERCATOR_CRS,
            ) 

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 get_web_optimized_params(
    src_dst,
    tilesize=256,
    latitude_adjustment: bool = True,
    warp_resampling: str = "nearest",
    grid_crs=CRS.from_epsg(3857),
) -> Dict:
    """Return VRT parameters for a WebOptimized COG."""
    bounds = list(
        transform_bounds(
            src_dst.crs, CRS.from_epsg(4326), *src_dst.bounds, densify_pts=21
        )
    )
    center = [(bounds[0] + bounds[2]) / 2, (bounds[1] + bounds[3]) / 2]

    lat = 0 if latitude_adjustment else center[1]
    max_zoom = get_max_zoom(src_dst, lat=lat, tilesize=tilesize)

    extrema = tile_extrema(bounds, max_zoom)

    left, _, _, top = mercantile.xy_bounds(
        extrema["x"]["min"], extrema["y"]["min"], max_zoom
    )
    vrt_res = _meters_per_pixel(max_zoom, 0, tilesize=tilesize)
    vrt_transform = Affine(vrt_res, 0, left, 0, -vrt_res, top)

    vrt_width = (extrema["x"]["max"] - extrema["x"]["min"]) * tilesize
    vrt_height = (extrema["y"]["max"] - extrema["y"]["min"]) * tilesize

    return dict(
        crs=grid_crs,
        transform=vrt_transform,
        width=vrt_width,
        height=vrt_height,
        resampling=ResamplingEnums[warp_resampling],
    ) 

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 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 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 get_tile_tm2(request, x, y, z):
    """
    """
    scale_denominator = zoom_to_scale_denom(z)

    # compute mercator bounds
    bounds = mercantile.xy_bounds(x, y, z)
    bbox = f"st_makebox2d(st_point({bounds.left}, {bounds.bottom}), st_point({bounds.right},{bounds.top}))"

    sql = Config.tm2query.format(
        bbox=bbox,
        scale_denominator=scale_denominator,
        pixel_width=256,
        pixel_height=256,
    )
    logger.debug(sql)

    async with Config.db_pg.acquire() as conn:
        # join tiles into one bytes string except null tiles
        rows = await conn.fetch(sql)
        pbf = b''.join([row[0] for row in rows if row[0]])

    return response.raw(
        pbf,
        headers={"Content-Type": "application/x-protobuf"}
    ) 

def get_tile_postgis(request, x, y, z, layer):
    """
    Direct access to a postgis layer
    """
    if ' ' in layer:
        return response.text('bad layer name: {}'.format(layer), status=404)

    # get fields given in parameters
    fields = ',' + request.raw_args['fields'] if 'fields' in request.raw_args else ''
    # get geometry column name from query args else geom is used
    geom = request.raw_args.get('geom', 'geom')
    # compute mercator bounds
    bounds = mercantile.xy_bounds(x, y, z)

    # make bbox for filtering
    bbox = f"st_setsrid(st_makebox2d(st_point({bounds.left}, {bounds.bottom}), st_point({bounds.right},{bounds.top})), {OUTPUT_SRID})"

    # compute pixel resolution
    scale = resolution(z)

    sql = single_layer.format(**locals(), OUTPUT_SRID=OUTPUT_SRID)

    logger.debug(sql)

    async with Config.db_pg.acquire() as conn:
        rows = await conn.fetch(sql)
        pbf = b''.join([row[0] for row in rows if row[0]])

    return response.raw(
        pbf,
        headers={"Content-Type": "application/x-protobuf"}
    ) 

def tile(
    src_dst: Union[DatasetReader, DatasetWriter, WarpedVRT],
    x: int,
    y: int,
    z: int,
    tilesize: int = 256,
    **kwargs,
) -> Tuple[numpy.ndarray, numpy.ndarray]:
    """
    Read mercator tile from an image.

    Attributes
    ----------
        src_dst : rasterio.io.DatasetReader
            rasterio.io.DatasetReader object
        x : int
            Mercator tile X index.
        y : int
            Mercator tile Y index.
        z : int
            Mercator tile ZOOM level.
        tilesize : int, optional
            Output tile size. Default is 256.
        kwargs : Any, optional
            Additional options to forward to part()

    Returns
    -------
        data : numpy ndarray
        mask: numpy array

    """
    bounds = transform_bounds(
        src_dst.crs, constants.WGS84_CRS, *src_dst.bounds, densify_pts=21
    )
    if not tile_exists(bounds, z, x, y):
        raise TileOutsideBounds(f"Tile {z}/{x}/{y} is outside image bounds")

    tile_bounds = mercantile.xy_bounds(mercantile.Tile(x=x, y=y, z=z))
    return part(
        src_dst,
        tile_bounds,
        tilesize,
        tilesize,
        dst_crs=constants.WEB_MERCATOR_CRS,
        **kwargs,
    ) 

def test_cog_translate_web():
    """
    Test Web-Optimized COG.

    - Test COG size is a multiple of 256 (mercator tile size)
    - Test COG bounds are aligned with mercator grid at max zoom
    """
    runner = CliRunner()
    with runner.isolated_filesystem():

        web_profile = cog_profiles.get("raw")
        web_profile.update({"blockxsize": 256, "blockysize": 256})
        config = dict(GDAL_TIFF_OVR_BLOCKSIZE="128")

        cog_translate(
            raster_path_web,
            "cogeo.tif",
            web_profile,
            quiet=True,
            web_optimized=True,
            config=config,
        )
        with rasterio.open(raster_path_web) as src_dst:
            with rasterio.open("cogeo.tif") as out_dst:
                blocks = list(set(out_dst.block_shapes))
                assert len(blocks) == 1
                ts = blocks[0][0]
                assert not out_dst.width % ts
                assert not out_dst.height % ts
                max_zoom = get_max_zoom(out_dst)

                bounds = list(
                    transform_bounds(
                        src_dst.crs, "epsg:4326", *src_dst.bounds, densify_pts=21
                    )
                )
                ulTile = mercantile.xy_bounds(
                    mercantile.tile(bounds[0], bounds[3], max_zoom)
                )
                assert out_dst.bounds.left == ulTile.left
                assert out_dst.bounds.top == ulTile.top

                lrTile = mercantile.xy_bounds(
                    mercantile.tile(bounds[2], bounds[1], max_zoom)
                )
                assert out_dst.bounds.right == lrTile.right
                assert round(out_dst.bounds.bottom, 6) == round(lrTile.bottom, 6)