Python rasterio.transform.from_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 _requested_tile_aligned_with_internal_tile(
    src_dst: Union[DatasetReader, DatasetWriter, WarpedVRT],
    bounds: Tuple[float, float, float, float],
    height: int,
    width: int,
) -> bool:
    """Check if tile is aligned with internal tiles."""
    if not src_dst.is_tiled:
        return False

    if src_dst.crs != constants.WEB_MERCATOR_CRS:
        return False

    col_off, row_off, w, h = windows.from_bounds(
        *bounds, height=height, transform=src_dst.transform, width=width
    ).flatten()

    if round(w) % 64 and round(h) % 64:
        return False

    if (src_dst.width - round(col_off)) % 64:
        return False

    if (src_dst.height - round(row_off)) % 64:
        return False

    return True 

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 get_overview_level(
    src_dst: Union[DatasetReader, DatasetWriter, WarpedVRT],
    bounds: Tuple[float, float, float, float],
    height: int,
    width: int,
    dst_crs: CRS = constants.WEB_MERCATOR_CRS,
) -> int:
    """
    Return the overview level corresponding to the tile resolution.

    Freely adapted from https://github.com/OSGeo/gdal/blob/41993f127e6e1669fbd9e944744b7c9b2bd6c400/gdal/apps/gdalwarp_lib.cpp#L2293-L2362

    Attributes
    ----------
        src_dst : rasterio.io.DatasetReader
            Rasterio io.DatasetReader object
        bounds : list
            Bounds (left, bottom, right, top) in target crs ("dst_crs").
        height : int
            Output height.
        width : int
            Output width.
        dst_crs: CRS or str, optional
            Target coordinate reference system (default "epsg:3857").

    Returns
    -------
        ovr_idx: Int or None
            Overview level

    """
    dst_transform, _, _ = calculate_default_transform(
        src_dst.crs, dst_crs, src_dst.width, src_dst.height, *src_dst.bounds
    )
    src_res = dst_transform.a

    # Compute what the "natural" output resolution
    # (in pixels) would be for this input dataset
    vrt_transform = from_bounds(*bounds, width, height)
    target_res = vrt_transform.a

    ovr_idx = -1
    if target_res > src_res:
        res = [src_res * decim for decim in src_dst.overviews(1)]

        for ovr_idx in range(ovr_idx, len(res) - 1):
            ovrRes = src_res if ovr_idx < 0 else res[ovr_idx]
            nextRes = res[ovr_idx + 1]
            if (ovrRes < target_res) and (nextRes > target_res):
                break
            if abs(ovrRes - target_res) < 1e-1:
                break
        else:
            ovr_idx = len(res) - 1

    return ovr_idx 

def get_vrt_transform(
    src_dst: Union[DatasetReader, DatasetWriter, WarpedVRT],
    bounds: Tuple[float, float, float, float],
    height: Optional[int] = None,
    width: Optional[int] = None,
    dst_crs: CRS = constants.WEB_MERCATOR_CRS,
) -> Tuple[Affine, int, int]:
    """
    Calculate VRT transform.

    Attributes
    ----------
    src_dst : rasterio.io.DatasetReader
        Rasterio io.DatasetReader object
    bounds : list
        Bounds (left, bottom, right, top) in target crs ("dst_crs").
    height : int, optional
        Desired output height of the array for the bounds.
    width : int, optional
        Desired output width of the array for the bounds.
    dst_crs: CRS or str, optional
        Target coordinate reference system (default "epsg:3857").

    Returns
    -------
    vrt_transform: Affine
        Output affine transformation matrix
    vrt_width, vrt_height: int
        Output dimensions

    """
    dst_transform, _, _ = calculate_default_transform(
        src_dst.crs, dst_crs, src_dst.width, src_dst.height, *src_dst.bounds
    )
    w, s, e, n = bounds

    if not height or not width:
        vrt_width = math.ceil((e - w) / dst_transform.a)
        vrt_height = math.ceil((s - n) / dst_transform.e)
        vrt_transform = from_bounds(w, s, e, n, vrt_width, vrt_height)
        return vrt_transform, vrt_width, vrt_height

    tile_transform = from_bounds(w, s, e, n, width, height)
    w_res = (
        tile_transform.a
        if abs(tile_transform.a) < abs(dst_transform.a)
        else dst_transform.a
    )
    h_res = (
        tile_transform.e
        if abs(tile_transform.e) < abs(dst_transform.e)
        else dst_transform.e
    )

    vrt_width = math.ceil((e - w) / w_res)
    vrt_height = math.ceil((s - n) / h_res)
    vrt_transform = from_bounds(w, s, e, n, vrt_width, vrt_height)

    return vrt_transform, vrt_width, vrt_height 

def _tile_worker(tile):
    """
    For each tile, and given an open rasterio src, plus a`global_args` dictionary
    with attributes of `base_val`, `interval`, and a `writer_func`,
    warp a continous single band raster to a 512 x 512 mercator tile,
    then encode this tile into RGB.

    Parameters
    -----------
    tile: list
        [x, y, z] indices of tile

    Returns
    --------
    tile, buffer
        tuple with the input tile, and a bytearray with the data encoded into
        the format created in the `writer_func`

    """
    x, y, z = tile

    bounds = [
        c
        for i in (
            mercantile.xy(*mercantile.ul(x, y + 1, z)),
            mercantile.xy(*mercantile.ul(x + 1, y, z)),
        )
        for c in i
    ]

    toaffine = transform.from_bounds(*bounds + [512, 512])

    out = np.empty((512, 512), dtype=src.meta["dtype"])

    reproject(
        rasterio.band(src, 1),
        out,
        dst_transform=toaffine,
        dst_crs="epsg:3857",
        resampling=Resampling.bilinear,
    )

    out = data_to_rgb(out, global_args["base_val"], global_args["interval"])

    return tile, global_args["writer_func"](out, global_args["kwargs"].copy(), toaffine) 

def _compute_image_stats(dataset: 'DatasetReader',
                             max_shape: Sequence[int] = None) -> Optional[Dict[str, Any]]:
        """Compute statistics for the given rasterio dataset by reading it into memory."""
        from rasterio import features, warp, transform
        from shapely import geometry

        out_shape = (dataset.height, dataset.width)

        if max_shape is not None:
            out_shape = (
                min(max_shape[0], out_shape[0]),
                min(max_shape[1], out_shape[1])
            )

        data_transform = transform.from_bounds(
            *dataset.bounds, height=out_shape[0], width=out_shape[1]
        )
        raster_data = dataset.read(1, out_shape=out_shape, masked=True)

        if dataset.nodata is not None:
            # nodata values might slip into output array if out_shape < dataset.shape
            raster_data[raster_data == dataset.nodata] = np.ma.masked

        valid_data = raster_data.compressed()

        if valid_data.size == 0:
            return None

        if np.any(raster_data.mask):
            hull_candidates = RasterDriver._hull_candidate_mask(~raster_data.mask)
            hull_shapes = (geometry.shape(s) for s, _ in features.shapes(
                np.ones(hull_candidates.shape, 'uint8'),
                mask=hull_candidates,
                transform=data_transform
            ))
            convex_hull = geometry.MultiPolygon(hull_shapes).convex_hull
        else:
            # no masked entries -> convex hull == dataset bounds
            w, s, e, n = dataset.bounds
            convex_hull = geometry.Polygon([(w, s), (e, s), (e, n), (w, n)])

        convex_hull_wgs = warp.transform_geom(
            dataset.crs, 'epsg:4326', geometry.mapping(convex_hull)
        )

        return {
            'valid_percentage': valid_data.size / raster_data.size * 100,
            'range': (float(valid_data.min()), float(valid_data.max())),
            'mean': float(valid_data.mean()),
            'stdev': float(valid_data.std()),
            'percentiles': np.percentile(valid_data, np.arange(1, 100)),
            'convex_hull': convex_hull_wgs
        } 

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 _rasterio_read(
    input_file=None,
    indexes=None,
    dst_bounds=None,
    dst_shape=None,
    dst_crs=None,
    resampling=None,
    src_nodata=None,
    dst_nodata=None,
):

    def _read(
        src, indexes, dst_bounds, dst_shape, dst_crs, resampling, src_nodata, dst_nodata
    ):
        height, width = dst_shape[-2:]
        if indexes is None:
            dst_shape = (len(src.indexes), height, width)
            indexes = list(src.indexes)
        src_nodata = src.nodata if src_nodata is None else src_nodata
        dst_nodata = src.nodata if dst_nodata is None else dst_nodata
        dst_left, dst_bottom, dst_right, dst_top = dst_bounds
        with WarpedVRT(
            src,
            crs=dst_crs,
            src_nodata=src_nodata,
            nodata=dst_nodata,
            width=width,
            height=height,
            transform=affine_from_bounds(
                dst_left, dst_bottom, dst_right, dst_top, width, height
            ),
            resampling=Resampling[resampling]
        ) as vrt:
            return vrt.read(
                window=vrt.window(*dst_bounds),
                out_shape=dst_shape,
                indexes=indexes,
                masked=True
            )
    if isinstance(input_file, str):
        logger.debug("got file path %s", input_file)
        try:
            with rasterio.open(input_file, "r") as src:
                return _read(
                    src, indexes, dst_bounds, dst_shape, dst_crs, resampling, src_nodata,
                    dst_nodata
                )
        except RasterioIOError as e:
            try:
                if path_exists(input_file):
                    raise e
            except:
                raise e
            raise FileNotFoundError("%s not found" % input_file)
    else:
        logger.debug("assuming file object %s", input_file)
        return _read(
            input_file, indexes, dst_bounds, dst_shape, dst_crs, resampling,
            src_nodata, dst_nodata
        )