Get the closest web mercator zoom for a raster input

get_zoom.py

import math
import sys

import rasterio
from rasterio.warp import (calculate_default_transform, transform)


CHUNK_SIZE = 256


def get_zoom(input, dst_crs="EPSG:3857"):
    with rasterio.drivers():
        with rasterio.open(input) as src:
            # Compute the geographic bounding box of the dataset.
            (west, east), (south, north) = transform(
                src.crs, "EPSG:4326", src.bounds[::2], src.bounds[1::2])

            affine, _, _ = calculate_default_transform(src.crs, dst_crs,
                src.width, src.height, *src.bounds, resolution=None)

            # grab the lowest resolution dimension
            resolution = max(abs(affine[0]), abs(affine[4]))

            return int(round(math.log((2 * math.pi * 6378137) /
                                      (resolution * CHUNK_SIZE)) / math.log(2)))


if __name__ == "__main__":
    print(get_zoom(sys.argv[1]))