Rio Warp + Overview

info.md

the test.py is equivalent to : gdalwarp -t_srs EPSG:3857 -te -8922952.933898335 5244191.636589372 -8883817.175416324 5283327.395071382 -te_srs EPSG:3857 -ts 128 128 -r bilinear /vsis3/landsat-pds/c1/L8/017/030/LC08_L1TP_017030_20180228_20180308_01_T1/LC08_L1TP_017030_20180228_20180308_01_T1_B1.TIF test.tif

VSI_CACHE=FALSE VSI_CACHE_SIZE=0 CPL_DEBUG=ON CPL_CURL_VERBOSE=YES gdalwarp -t_srs EPSG:3857 -te -8922952.933898335 5244191.636589372 -8883817.175416324 5283327.395071382  -te_srs EPSG:3857 -ts 128 128 -r bilinear /vsis3/landsat-pds/c1/L8/017/030/LC08_L1TP_017030_20180228_20180308_01_T1/LC08_L1TP_017030_20180228_20180308_01_T1_B1.TIF test.tif 2>&1 >/dev/null | grep "Content-Length:"
< Content-Length: 16384
< Content-Length: 8700
< Content-Length: 8700
< Content-Length: 16384
< Content-Length: 358435
< Content-Length: 119216

 (py4) ~/Desktop  VSI_CACHE=FALSE VSI_CACHE_SIZE=0 CPL_DEBUG=ON CPL_CURL_VERBOSE=YES python test.py  2>&1 >/dev/null | grep "Content-Length:"
< Content-Length: 16384
< Content-Length: 16384
< Content-Length: 358435
< Content-Length: 119216

test.py

import numpy as np
import math
import rasterio
import mercantile
from rasterio import transform
from rasterio.warp import transform_bounds, calculate_default_transform
from rasterio.enums import Resampling


def fake_boundless_ovr():
    dataset = 's3://landsat-pds/c1/L8/017/030/LC08_L1TP_017030_20180228_20180308_01_T1/LC08_L1TP_017030_20180228_20180308_01_T1_B1.TIF'
    with rasterio.open(dataset) as src:

        wgs_bounds = transform_bounds(*[src.crs, 'epsg:4326'] + list(src.bounds), densify_pts=21)

        # Get list of Z10 tiles covering the data 
        tiles = list(mercantile.tiles(*wgs_bounds, [10]))
        
        # Select on tile to read
        tile_bounds = mercantile.xy_bounds(tiles[6])
        
        dst_transform, _, _ = calculate_default_transform(src.crs,
                                                          'epsg:3857',
                                                          src.width,
                                                          src.height,
                                                          *src.bounds)

        # Get Full resolution Mercator tile size
        w, s, e, n = tile_bounds
        vrt_width = math.ceil((e - w) / dst_transform.a)
        vrt_height = math.ceil((s - n) / dst_transform.e)

        # Get VRT High resolution transform
        vrt_transform = transform.from_bounds(*tile_bounds, vrt_width, vrt_height)

        vrt_params = dict(
            crs='EPSG:3857',
            resampling=Resampling.bilinear,
            transform=dst_transform,
            width=vrt_width,
            height=vrt_height)

        tilesize = 128
        with rasterio.vrt.WarpedVRT(src, **vrt_params) as vrt:
            # Force decimated read
            data = vrt.read(out_shape=(1, tilesize, tilesize), resampling=Resampling.bilinear)

        # Save result to disk
        meta = {
            'driver': 'GTiff',
            'count': 1,
            'dtype': np.uint16,
            'height': tilesize,
            'width': tilesize,
            'crs': 'epsg:3857',
            'transform': transform.from_bounds(*tile_bounds, tilesize, tilesize}

        with rasterio.open('test_rio.tif', 'w', **meta) as dataset:
            dataset.write(data)


if __name__ == '__main__':
    fake_boundless_ovr()