Python rasterio.warp.transform_bounds() Examples

def spatial_info(address: str) -> Dict:
    """
    Return COGEO spatial info.

    Attributes
    ----------
        address : str or PathLike object
            A dataset path or URL. Will be opened in "r" mode.

    Returns
    -------
        out : dict.

    """
    with rasterio.open(address) as src_dst:
        minzoom, maxzoom = get_zooms(src_dst)
        bounds = transform_bounds(
            src_dst.crs, constants.WGS84_CRS, *src_dst.bounds, densify_pts=21
        )
        center = [(bounds[0] + bounds[2]) / 2, (bounds[1] + bounds[3]) / 2, minzoom]

    return dict(
        address=address, bounds=bounds, center=center, minzoom=minzoom, maxzoom=maxzoom
    ) 

def bounds(address: str) -> Dict:
    """
    Retrieve image bounds.

    Attributes
    ----------
        address : str
            file url.

    Returns
    -------
        out : dict
            dictionary with image bounds.

    """
    with rasterio.open(address) as src_dst:
        bounds = transform_bounds(
            src_dst.crs, constants.WGS84_CRS, *src_dst.bounds, densify_pts=21
        )
    return dict(address=address, bounds=bounds) 

def bounds(sceneid: str) -> Dict:
    """
    Retrieve image bounds.

    Attributes
    ----------
        sceneid : str
            Sentinel-2 sceneid.

    Returns
    -------
        out : dict
            dictionary with image bounds.

    """
    scene_params = sentinel2_parser(sceneid)
    preview_file = "{scheme}://{bucket}/{prefix}/{preview_file}".format(**scene_params)
    with rasterio.open(preview_file) as src_dst:
        bounds = transform_bounds(
            src_dst.crs, constants.WGS84_CRS, *src_dst.bounds, densify_pts=21
        )

    return dict(sceneid=sceneid, bounds=bounds) 

def tilejson_handler(url: str, tile_format: str = "png", **kwargs: Dict):
    """Handle /tilejson.json requests."""
    qs = urllib.parse.urlencode(list(kwargs.items()))
    tile_url = f"{APP.host}/tiles/{{z}}/{{x}}/{{y}}.{tile_format}?url={url}"
    if qs:
        tile_url += f"&{qs}"

    with rasterio.open(url) as src_dst:
        bounds = warp.transform_bounds(
            src_dst.crs, "epsg:4326", *src_dst.bounds, densify_pts=21
        )
        center = [(bounds[0] + bounds[2]) / 2, (bounds[1] + bounds[3]) / 2]
        minzoom, maxzoom = get_zooms(src_dst)

    meta = dict(
        bounds=bounds,
        center=center,
        minzoom=minzoom,
        maxzoom=maxzoom,
        name=os.path.basename(url),
        tilejson="2.1.0",
        tiles=[tile_url],
    )
    return ("OK", "application/json", json.dumps(meta)) 

def tilejson_handler(
    event: Dict,
    sceneid: str,
    tile_format: str = "png",
    tile_scale: int = 1,
    **kwargs: Any,
) -> Tuple[str, str, str]:
    """Handle /tilejson.json requests."""
    # HACK
    token = event["multiValueQueryStringParameters"].get("access_token")
    if token:
        kwargs.update(dict(access_token=token[0]))

    qs = urllib.parse.urlencode(list(kwargs.items()))
    tile_url = (
        f"{APP.host}/tiles/{sceneid}/{{z}}/{{x}}/{{y}}@{tile_scale}x.{tile_format}?{qs}"
    )

    scene_params = landsat8._landsat_parse_scene_id(sceneid)
    landsat_address = f"{LANDSAT_BUCKET}/{scene_params['key']}_BQA.TIF"
    with rasterio.open(landsat_address) as src_dst:
        bounds = warp.transform_bounds(
            src_dst.crs, "epsg:4326", *src_dst.bounds, densify_pts=21
        )
        minzoom, maxzoom = get_zooms(src_dst)
        center = [(bounds[0] + bounds[2]) / 2, (bounds[1] + bounds[3]) / 2, minzoom]

    meta = dict(
        bounds=bounds,
        center=center,
        minzoom=minzoom,
        maxzoom=maxzoom,
        name=sceneid,
        tilejson="2.1.0",
        tiles=[tile_url],
    )
    return ("OK", "application/json", json.dumps(meta)) 

def tilejson_handler(
    url: str, tile_format: str = "png", tile_scale: int = 1, **kwargs: Any
) -> Tuple[str, str, str]:
    """Handle /tilejson.json requests."""
    kwargs.update(dict(url=url))
    qs = urllib.parse.urlencode(list(kwargs.items()))
    tile_url = f"{APP.host}/tiles/{{z}}/{{x}}/{{y}}@{tile_scale}x.{tile_format}"
    if qs:
        tile_url += f"?{qs}"

    with rasterio.open(url) as src_dst:
        bounds = warp.transform_bounds(
            src_dst.crs, "epsg:4326", *src_dst.bounds, densify_pts=21
        )
        minzoom, maxzoom = get_zooms(src_dst)
        center = [(bounds[0] + bounds[2]) / 2, (bounds[1] + bounds[3]) / 2, minzoom]

    meta = dict(
        bounds=bounds,
        center=center,
        minzoom=minzoom,
        maxzoom=maxzoom,
        name=os.path.basename(url),
        tilejson="2.1.0",
        tiles=[tile_url],
    )
    return ("OK", "application/json", json.dumps(meta)) 

def tilejson_handler(
    event: Dict,
    scene: str,
    tile_format: str = "png",
    tile_scale: int = 1,
    **kwargs: Any,
) -> Tuple[str, str, str]:
    """Handle /tilejson.json requests."""
    # HACK
    token = event["multiValueQueryStringParameters"].get("access_token")
    if token:
        kwargs.update(dict(access_token=token[0]))

    qs = urllib.parse.urlencode(list(kwargs.items()))
    tile_url = f"{APP.host}/s2/tiles/{scene}/{{z}}/{{x}}/{{y}}@{tile_scale}x.{tile_format}?{qs}"

    scene_params = sentinel2._sentinel_parse_scene_id(scene)
    sentinel_address = "s3://{}/{}/B{}.jp2".format(
        sentinel2.SENTINEL_BUCKET, scene_params["key"], "04"
    )
    with rasterio.open(sentinel_address) as src_dst:
        bounds = warp.transform_bounds(
            *[src_dst.crs, "epsg:4326"] + list(src_dst.bounds), densify_pts=21
        )
        minzoom, maxzoom = get_zooms(src_dst)
        center = [(bounds[0] + bounds[2]) / 2, (bounds[1] + bounds[3]) / 2, minzoom]

    meta = dict(
        bounds=bounds,
        center=center,
        minzoom=minzoom,
        maxzoom=maxzoom,
        name=scene,
        tilejson="2.1.0",
        tiles=[tile_url],
    )
    return ("OK", "application/json", json.dumps(meta)) 

def get_tile_wms(tile, imagery, folder, kwargs):
    """
    Read a WMS endpoint with query parameters corresponding to a TMS tile

    Converts the tile boundaries to the spatial/coordinate reference system
    (SRS or CRS) specified by the WMS query parameter.
    """
    # retrieve the necessary parameters from the query string
    query_dict = parse_qs(imagery.lower())
    image_format = query_dict.get('format')[0].split('/')[1]
    wms_version = query_dict.get('version')[0]
    if wms_version == '1.3.0':
        wms_srs = query_dict.get('crs')[0]
    else:
        wms_srs = query_dict.get('srs')[0]

    # find our tile bounding box
    bound = bounds(*[int(t) for t in tile.split('-')])
    xmin, ymin, xmax, ymax = transform_bounds(WGS84_CRS, CRS.from_string(wms_srs), *bound, densify_pts=21)

    # project the tile bounding box from lat/lng to WMS SRS
    bbox = (
        [ymin, xmin, ymax, xmax] if wms_version == "1.3.0" else [xmin, ymin, xmax, ymax]
    )

    # request the image with the transformed bounding box and save
    wms_url = imagery.replace('{bbox}', ','.join([str(b) for b in bbox]))
    r = requests.get(wms_url, auth=kwargs.get('http_auth'))
    tile_img = op.join(folder, '{}.{}'.format(tile, image_format))
    with open(tile_img, 'wb') as w:
        w.write(r.content)
    return tile_img 

def _make_tiles(bbox, src_crs, minz, maxz):
    """
    Given a bounding box, zoom range, and source crs,
    find all tiles that would intersect

    Parameters
    -----------
    bbox: list
        [w, s, e, n] bounds
    src_crs: str
        the source crs of the input bbox
    minz: int
        minumum zoom to find tiles for
    maxz: int
        maximum zoom to find tiles for

    Returns
    --------
    tiles: generator
        generator of [x, y, z] tiles that intersect
        the provided bounding box
    """
    w, s, e, n = transform_bounds(*[src_crs, "epsg:4326"] + bbox, densify_pts=0)

    EPSILON = 1.0e-10

    w += EPSILON
    s += EPSILON
    e -= EPSILON
    n -= EPSILON

    for z in range(minz, maxz + 1):
        for x, y in _tile_range(mercantile.tile(w, n, z), mercantile.tile(e, s, z)):
            yield [x, y, z] 

def get_dataset_info(src_path: str) -> Dict:
    """Get rasterio dataset meta."""
    with rasterio.open(src_path) as src_dst:
        bounds = transform_bounds(
            src_dst.crs, "epsg:4326", *src_dst.bounds, densify_pts=21
        )
        min_zoom, max_zoom = get_zooms(src_dst, ensure_global_max_zoom=True)
        return {
            "geometry": {
                "type": "Polygon",
                "coordinates": [
                    [
                        [bounds[0], bounds[3]],
                        [bounds[0], bounds[1]],
                        [bounds[2], bounds[1]],
                        [bounds[2], bounds[3]],
                        [bounds[0], bounds[3]],
                    ]
                ],
            },
            "properties": {
                "path": src_path,
                "bounds": bounds,
                "minzoom": min_zoom,
                "maxzoom": max_zoom,
                "datatype": src_dst.meta["dtype"],
            },
            "type": "Feature",
        } 

def get_bounds(geo_obj, crs=None):
    """Get the ``[left, bottom, right, top]`` bounds in any CRS.

    Arguments
    ---------
    geo_obj : a georeferenced raster or vector dataset.
    crs : int, optional
        The EPSG code (or other CRS format supported by rasterio.warp)
        for the CRS the bounds should be returned in. If not provided,
        the bounds will be returned in the same crs as `geo_obj`.

    Returns
    -------
    bounds : list
        ``[left, bottom, right, top]`` bounds in the input crs (if `crs` is
        ``None``) or in `crs` if it was provided.
    """
    input_data, input_type = _parse_geo_data(geo_obj)
    if input_type == 'vector':
        bounds = list(input_data.geometry.total_bounds)
    elif input_type == 'raster':
        if isinstance(input_data, rasterio.DatasetReader):
            bounds = list(input_data.bounds)
        elif isinstance(input_data, gdal.Dataset):
            input_gt = input_data.GetGeoTransform()
            min_x = input_gt[0]
            max_x = min_x + input_gt[1]*input_data.RasterXSize
            max_y = input_gt[3]
            min_y = max_y + input_gt[5]*input_data.RasterYSize

            bounds = [min_x, min_y, max_x, max_y]

    if crs is not None:
        crs = _check_crs(crs)
        src_crs = get_crs(input_data)
        # transform bounds to desired CRS
        bounds = transform_bounds(src_crs.to_wkt("WKT1_GDAL"),
                                  crs.to_wkt("WKT1_GDAL"), *bounds)

    return bounds 

def _wmts(
    mosaicid: str = None,
    url: str = None,
    tile_format: str = "png",
    tile_scale: int = 1,
    title: str = "Cloud Optimizied GeoTIFF Mosaic",
    **kwargs: Any,
) -> Tuple[str, str, str]:
    """Handle /wmts requests."""
    if tile_scale is not None and isinstance(tile_scale, str):
        tile_scale = int(tile_scale)

    kwargs.pop("SERVICE", None)
    kwargs.pop("REQUEST", None)
    kwargs.update(dict(url=url))
    query_string = urllib.parse.urlencode(list(kwargs.items()))
    query_string = query_string.replace(
        "&", "&"
    )  # & is an invalid character in XML

    with rasterio.open(url) as src_dst:
        bounds = warp.transform_bounds(
            src_dst.crs, "epsg:4326", *src_dst.bounds, densify_pts=21
        )
        minzoom, maxzoom = get_zooms(src_dst)

    return (
        "OK",
        "application/xml",
        wmts_template(
            f"{APP.host}",
            os.path.basename(url),
            query_string,
            minzoom=minzoom,
            maxzoom=maxzoom,
            bounds=bounds,
            tile_scale=tile_scale,
            tile_format=tile_format,
            title=title,
        ),
    ) 

def get_zooms(
    src_dst: Union[DatasetReader, DatasetWriter, WarpedVRT],
    ensure_global_max_zoom: bool = False,
    tilesize: int = 256,
) -> Tuple[int, int]:
    """
    Calculate raster min/max mercator zoom level.

    Parameters
    ----------
    src_dst: rasterio.io.DatasetReader
        Rasterio io.DatasetReader object
    ensure_global_max_zoom: bool, optional
        Apply latitude correction factor to ensure max_zoom equality for global
        datasets covering different latitudes (default: False).
    tilesize: int, optional
        Mercator tile size (default: 256).

    Returns
    -------
    min_zoom, max_zoom: Tuple
        Min/Max Mercator zoom levels.

    """
    bounds = transform_bounds(
        src_dst.crs, constants.WGS84_CRS, *src_dst.bounds, densify_pts=21
    )
    center = [(bounds[0] + bounds[2]) / 2, (bounds[1] + bounds[3]) / 2]
    lat = center[1] if ensure_global_max_zoom else 0

    dst_affine, w, h = calculate_default_transform(
        src_dst.crs,
        constants.WEB_MERCATOR_CRS,
        src_dst.width,
        src_dst.height,
        *src_dst.bounds,
    )

    mercator_resolution = max(abs(dst_affine[0]), abs(dst_affine[4]))

    # Correction factor for web-mercator projection latitude scale change
    latitude_correction_factor = math.cos(math.radians(lat))
    adjusted_resolution = mercator_resolution * latitude_correction_factor

    max_zoom = zoom_for_pixelsize(adjusted_resolution, tilesize=tilesize)

    ovr_resolution = adjusted_resolution * max(h, w) / tilesize
    min_zoom = zoom_for_pixelsize(ovr_resolution, tilesize=tilesize)

    return (min_zoom, max_zoom) 

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 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()