README.md

Test

GeoTIFF: https://s3-us-west-2.amazonaws.com/remotepixel-pub/data/Int16_nodata9999.tif

$ docker run --name dockerimage -p 8000:8000 --volume $(pwd)/:/local --rm -it lambci/lambda:build-python3.6 bash

$ yum install wget

$ cd /local

$ pip install pip -U
$ pip install rasterio rio-tiler

$ wget https://gist.githubusercontent.com/vincentsarago/0b229741c3261a9a447c7ab5dbd9eb05/raw/e03871160d0a8016534aeb39449833aae86e9a13/get_tile_border.py

# create 4 geotiff in your local directory
$ python get_tile_border.py

get_tile_border.py

"""test."""

import os
from rio_tiler.utils import array_to_image, get_vrt_transform

import mercantile

import rasterio
from rasterio.enums import Resampling
from rasterio.transform import from_bounds
from rasterio.vrt import WarpedVRT

input = "https://s3-us-west-2.amazonaws.com/remotepixel-pub/data/Int16_nodata9999.tif"
tile_z = 12
tile_x = 676
tile_y = 1620


def tile_read(src_path, bounds, tilesize, warp_resampling, read_resampling):
    """~Copy from rio-tiler."""
    with rasterio.open(src_path) as src_dst:
        vrt_transform, vrt_width, vrt_height = get_vrt_transform(src_dst, bounds)
        vrt_params = dict(
            crs="epsg:3857", 
            transform=vrt_transform,
            width=vrt_width,
            height=vrt_height,
            nodata=src_dst.nodata,
            src_nodata=src_dst.nodata,
            add_alpha=False,
            resampling=Resampling[warp_resampling],
        )

        out_shape = (len(src_dst.indexes), tilesize, tilesize)
        with WarpedVRT(src_dst, **vrt_params) as vrt:
            data = vrt.read(
                out_shape=out_shape,
                resampling=Resampling[read_resampling],
            )
            mask = vrt.dataset_mask(out_shape=(tilesize, tilesize))
            return data, mask


def main(warp_resampling, read_resampling):
    mercator_tile = mercantile.Tile(x=tile_x, y=tile_y, z=tile_z)
    bounds = mercantile.xy_bounds(mercator_tile)
    w, s, e, n = bounds
    
    tilesize = 256
    
    dst_crs = {'init': 'EPSG:3857'}
    dst_transform = from_bounds(w, s, e, n, tilesize, tilesize)
 
    tile, mask = tile_read(
        input, bounds, tilesize, warp_resampling, read_resampling
    )
    with open(f"{tile_z}-{tile_x}-{tile_y}_{warp_resampling}_{read_resampling}.tif", "wb") as f:
        f.write(array_to_image(
            tile,
            dtype=tile.dtype,
            img_format="GTiff",
            crs=dst_crs,
            transform=dst_transform
        ))


if __name__ == '__main__':
    main("nearest", "nearest")
    main("nearest", "bilinear")
    main("bilinear", "bilinear")
    main("bilinear", "nearest")

#     fout = f"{tile_z}-{tile_y}-{tile_y}_nearest_gdal.tif"
#     os.remove(fout)
#     cmd = f"gdalwarp {input} " + \
#           f"{fout} -t_srs EPSG:3857 " + \
#           "-te -13423565.159329513 4177742.217954592 -13413781.219709009 4187526.1575750955 " + \
#           "-te_srs EPSG:3857 -r near -ts 256 256"
#     os.system(cmd)

#     fout = f"{tile_z}-{tile_y}-{tile_y}_bilinear_gdal.tif"
#     os.remove(fout)
#     cmd = f"gdalwarp {input} " + \
#           f"{fout} -t_srs EPSG:3857 " + \
#           "-te -13423565.159329513 4177742.217954592 -13413781.219709009 4187526.1575750955 " + \
#           "-te_srs EPSG:3857 -r bilinear -ts 256 256"
#     os.system(cmd)