How to use the rasterio.warp.transform_bounds function in rasterio

To help you get started, we’ve selected a few rasterio examples, based on popular ways it is used in public projects.

Secure your code as it's written. Use Snyk Code to scan source code in minutes - no build needed - and fix issues immediately.

github mapbox / rasterio / tests / test_warp.py View on Github external
def test_transform_bounds():
    with rasterio.open('tests/data/RGB.byte.tif') as src:
        l, b, r, t = src.bounds
        assert np.allclose(
            transform_bounds(src.crs, {'init': 'EPSG:4326'}, l, b, r, t),
            (
                -78.95864996545055, 23.564991210854686,
                -76.57492370013823, 25.550873767433984
            )
github mapbox / rasterio / tests / test_warp.py View on Github external
def test_transform_bounds_src_crs_none():
    with pytest.raises(CRSError):
        transform_bounds(None, WGS84_crs, 0, 0, 0, 0)
github cogeotiff / rio-cogeo / tests / test_web.py View on Github external
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"] + list(src_dst.bounds),
                        densify_pts=21
                    )
                )

                leftTile = mercantile.tile(bounds[0], bounds[3], max_zoom)
                tbounds = mercantile.xy_bounds(leftTile)
                west, north = tbounds.left, tbounds.top
                assert out_dst.transform.xoff == west
                assert out_dst.transform.yoff == north

                rightTile = mercantile.tile(bounds[2], bounds[1], max_zoom)
                tbounds = mercantile.xy_bounds(rightTile)
                east, south = tbounds.right, tbounds.bottom

                lrx = round(
github RemotePixel / remotepixel-tiler / remotepixel_tiler / landsat.py View on Github external
) -> 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))
github opendatacube / datacube-core / datacube / api.py View on Github external
def geospatial_warp_bounds(input_coord, input_crs='EPSG:4326', output_crs='EPSG:4326', tolerance=0.):
    '''
    Converts coordinates, adding tolerance if they are the same point for index searching
    :param input_crs: EPSG or other GDAL string
    :param output_crs: EPSG or other GDAL string (Search is 'EPSG:4326')
    :param input_coord: {'left': float, 'bottom': float, 'right': float, 'top': float}
    :return: {'lat':Range,'lon':Range}
    '''
    if any(v is None for v in input_coord.values()):
        raise ValueError('Missing coordinate in input_coord {}'.format(input_coord))
    left, bottom, right, top = rasterio.warp.transform_bounds(input_crs, output_crs, **input_coord)
    output_coords = {'left': left, 'bottom': bottom, 'right': right, 'top': top}

    if bottom == top:
        output_coords['top'] = top - tolerance
        output_coords['bottom'] = bottom + tolerance
    if left == right:
        output_coords['left'] = left - tolerance
        output_coords['right'] = right + tolerance
    return output_coords
github opendatacube / datacube-core / datacube / api / _conversion.py View on Github external
def geospatial_warp_bounds(input_coord, input_crs='EPSG:4326', output_crs='EPSG:4326', tolerance=0.):
    '''
    Converts coordinates, adding tolerance if they are the same point for index searching
    :param input_crs: EPSG or other GDAL string
    :param output_crs: EPSG or other GDAL string (Search is 'EPSG:4326')
    :param input_coord: {'left': float, 'bottom': float, 'right': float, 'top': float}
    :return: {'lat':Range,'lon':Range}
    '''
    if any(v is None for v in input_coord.values()):
        raise ValueError('Unable to convert spatial bounds: Missing coordinate in input_coord {}'.format(input_coord))
    left, bottom, right, top = rasterio.warp.transform_bounds(input_crs, output_crs, **input_coord)
    output_coords = {'left': left, 'bottom': bottom, 'right': right, 'top': top}

    if bottom == top:
        output_coords['top'] = top + tolerance
        output_coords['bottom'] = bottom - tolerance
    if left == right:
        output_coords['left'] = left - tolerance
        output_coords['right'] = right + tolerance
    return output_coords
github cogeotiff / rio-tiler / rio_tiler / aws.py View on Github external
tileformat : str
        Image format to return (Accepted: "jpg" or "png")
    rgb : tuple, int, optional (default: (1, 2, 3))
        Bands index for the RGB combination.
    tilesize : int, optional (default: 256)
        Output image size.

    Returns
    -------
    out : numpy ndarray
    """

    source_address = '{}/{}/{}'.format(prefix, bucket, key)

    with rasterio.open(source_address) as src:
        wgs_bounds = transform_bounds(
            *[src.crs, 'epsg:4326'] + list(src.bounds), densify_pts=21)
        nodata = src.nodata

    if not utils.tile_exists(wgs_bounds, tile_z, tile_x, tile_y):
        raise TileOutsideBounds(
            'Tile {}/{}/{} is outside image bounds'.format(
                tile_z, tile_x, tile_y))

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

    w, s, e, n = tile_bounds

    out = utils.tile_band_worker(source_address,
                                 tile_bounds,
                                 tilesize,
github cogeotiff / rio-tiler / rio_tiler / sentinel2.py View on Github external
kwargs : optional
        These are passed to 'rio_tiler.sentinel2._sentinel_stats'
        e.g: histogram_bins=20'

    Returns
    -------
    out : dict
        Dictionary with image bounds and bands statistics.

    """
    scene_params = _sentinel_parse_scene_id(sceneid)
    sentinel_address = "{}/{}".format(SENTINEL_BUCKET, scene_params["key"])

    dst_crs = CRS({"init": "EPSG:4326"})
    with rasterio.open("{}/preview.jp2".format(sentinel_address)) as src:
        bounds = transform_bounds(
            *[src.crs, dst_crs] + list(src.bounds), densify_pts=21
        )

    info = {"sceneid": sceneid}
    info["bounds"] = {"value": bounds, "crs": dst_crs.to_string()}

    addresses = [
        "{}/preview/B{}.jp2".format(sentinel_address, band) for band in SENTINEL_BANDS
    ]

    _stats_worker = partial(_sentinel_stats, percentiles=(pmin, pmax), **kwargs)
    with futures.ThreadPoolExecutor(max_workers=MAX_THREADS) as executor:
        responses = executor.map(_stats_worker, addresses)
    info["statistics"] = {
        b: v for b, d in zip(SENTINEL_BANDS, responses) for k, v in d.items()
    }
github mapbox / rasterio / rasterio / rio / features.py View on Github external
def __call__(self):
            for i, path in enumerate(input):
                with rasterio.open(path) as src:
                    bounds = src.bounds
                    if dst_crs:
                        bbox = transform_bounds(src.crs,
                                                dst_crs, *bounds)
                    elif projection == 'mercator':
                        bbox = transform_bounds(src.crs,
                                                {'init': 'epsg:3857'}, *bounds)
                    elif projection == 'geographic':
                        bbox = transform_bounds(src.crs,
                                                {'init': 'epsg:4326'}, *bounds)
                    else:
                        bbox = bounds

                if precision >= 0:
                    bbox = [round(b, precision) for b in bbox]

                yield {
                    'type': 'Feature',
                    'bbox': bbox,
                    'geometry': {
                        'type': 'Polygon',
                        'coordinates': [[
                            [bbox[0], bbox[1]],
                            [bbox[2], bbox[1]],
                            [bbox[2], bbox[3]],