zoom_for_poly_wgs84.py

_ogc_crs84 = 'urn:ogc:def:crs:OGC:1.3:CRS84'
_web_mercator_epsg = 'EPSG:3857'


def _meters_per_pixel(zoom: int, lat: float=0.0, tile_size:int=256) -> float:
    """
    Return the pixel resolution for a given mercator tile zoom and lattitude.

    :param zoom: Mercator zoom level
    :param lat: Latitude in decimal degree. Optional (default: 0)
    :param tile_size: Mercator tile size (default: 256).
    :type tile_size: int, optional
    :return: Pixel resolution in meters
    :rtype: float
    """
    return (math.cos(lat * math.pi / 180.0) * 2 * math.pi * 6378137) / (
            tile_size * 2 ** zoom
    )


def zoom_for_pixelsize(pixel_size, max_z=24, tile_size=256, lat=0) -> int:
    """
    Get mercator zoom level corresponding to a pixel resolution.

    Freely adapted from
    https://github.com/OSGeo/gdal/blob/b0dfc591929ebdbccd8a0557510c5efdb893b852/gdal/swig/python/scripts/gdal2tiles.py#L294

    :param pixel_size: Pixel size in meters
    :param max_z: Max mercator zoom level allowed. Optional (default: 24)
    :param tile_size: Mercator tile size. Optional (default: 256).
    :param lat: Latitude of pixel size. Optional (default: 0)
    :return: Mercator zoom level corresponding to the pixel resolution
    """
    for z in range(max_z):
        if pixel_size > _meters_per_pixel(z, lat, tile_size=tile_size):
            return max(0, z - 1)  # We don't want to scale up

    return max_z - 1


def zoom_for_poly_wgs84(poly_wgs84: Polygon, map_size: int = 900) -> int:
    """
    Computes the **approximate** zoom value needed to cover the given polygon in a web map.

    Parameters
    :param poly_wgs84: Polygon in wgs84
    :param map_size: approx map widget size in pixels. optional
    :return: zoom to use to show polgyon
    """
    wm_poly = reproject(poly_wgs84, _ogc_crs84, _web_mercator_epsg)
    (xmin, ymin, xmax, ymax) = wm_poly.bounds
    width_m = xmax - xmin
    pixel_size = width_m / map_size
    return zoom_for_pixelsize(pixel_size, lat=poly_wgs84.centroid.y)