kylebarron / usgs-topo-tiler

Python package to read Web Mercator map tiles from USGS Historical Topographic Maps
https://kylebarron.dev/usgs-topo-mosaic
MIT License
77 stars 3 forks source link

Find min/max zoom based on metadata #9

Closed kylebarron closed 4 years ago

kylebarron commented 4 years ago

rio-tiler has a nice helper function get_zooms:

def get_zooms(src_dst, ensure_global_max_zoom=False, tilesize=256):
    """
    Calculate raster min/max mercator zoom level.
    Parameters
    ----------
    src_dst: rasterio.io.DatasetReader
        Rasterio io.DatasetReader object
    ensure_global_max_zoom: bool, optional
        Apply latitude correction factor to ensure max_zoom equality for global
        datasets covering different latitudes (default: False).
    tilesize: int, optional
        Mercator tile size (default: 256).
    Returns
    -------
    min_zoom, max_zoom: Tuple
        Min/Max Mercator zoom levels.
    """
    bounds = transform_bounds(src_dst.crs, "epsg:4326", *src_dst.bounds, densify_pts=21)
    center = [(bounds[0] + bounds[2]) / 2, (bounds[1] + bounds[3]) / 2]
    lat = center[1] if ensure_global_max_zoom else 0

    dst_affine, w, h = calculate_default_transform(
        src_dst.crs, "epsg:3857", src_dst.width, src_dst.height, *src_dst.bounds
    )

    mercator_resolution = max(abs(dst_affine[0]), abs(dst_affine[4]))

    # Correction factor for web-mercator projection latitude scale change
    latitude_correction_factor = math.cos(math.radians(lat))
    adjusted_resolution = mercator_resolution * latitude_correction_factor

    max_zoom = zoom_for_pixelsize(adjusted_resolution, tilesize=tilesize)

    ovr_resolution = adjusted_resolution * max(h, w) / tilesize
    min_zoom = zoom_for_pixelsize(ovr_resolution, tilesize=tilesize)

    return (min_zoom, max_zoom)

Essentially the only things you need are the map's CRS and bounds. Bounds are clear in the bulk csv metadata. There's also a projection column with the following unique values:

Polyconic                        115709
Lambert Conformal Conic           26067
Transverse Mercator               17647
Universal Transverse Mercator     17118
Unstated                           9648

So you should see if you can easily get a CRS from those

kylebarron commented 4 years ago

Should be able to get the pixel resolution (i.e. meters per pixel) from the DPI + scale: https://gis.stackexchange.com/questions/85307/how-to-calculate-size-of-a-pixel-for-a-scanned-map

Then pass that pixel resolution to zoom_for_pixelsize with tilesize=512 to get the max zoom.

rio-tiler gets the minzoom from

ovr_resolution = adjusted_resolution * max(h, w) / tilesize
min_zoom = zoom_for_pixelsize(ovr_resolution, tilesize=tilesize)

Or you could assume that all COGs have 6-7 levels of overviews? And set minzoom = maxzoom - 7?