Source-tile intersection queries using Affine projection matrices

affine_reproj_tile_intersect.py

def intersect_bounds(bounds1, bounds2):
    """Find intersection of two bounding boxes
    Attributes
    ----------
    bounds1 : list
        bounds (left, bottom, right, top).
    bounds2 : list
        bounds (left, bottom, right, top).
    Returns
    -------
    bounds : list or None
        bounds (left, bottom, right, top).
        Returns None if the boxes do not intersect
    """

    minx = max(bounds1[0], bounds2[0])
    miny = max(bounds1[1], bounds2[1])
    maxx = min(bounds1[2], bounds2[2])
    maxy = min(bounds1[3], bounds2[3])

    if minx >= maxx or miny >= maxy:
        return None

    return [minx, miny, maxx, maxy]


asset = 's3://naip-visualization/ca/2018/60cm/rgb/34118/m_3411861_ne_11_060_20180723_20190208.tif'
src = rasterio.open(asset)
dst_crs = CRS({"init": "EPSG:3857"})
transform, _, _ = calculate_default_transform(
    src.crs, dst_crs, src.width, src.height, *src.bounds)

corners_merc = []
for _x in [0, src.width]:
    for _y in [0, src.height]:
        corners_merc.append(transform * (_x, _y))

merc_asset_bbox = [
    min([t[0] for t in corners_merc]),
    min([t[1] for t in corners_merc]),
    max([t[0] for t in corners_merc]),
    max([t[1] for t in corners_merc]),
]

# Intersecting tile
x = 701
y = 1635
z = 12
tile_bounds = mercantile.xy_bounds(mercantile.Tile(x, y, z))
intersect_bounds(tile_bounds, merc_asset_bbox)
# [-13178966.668816948,
#  4037035.4799059597,
#  -13177246.953488544,
#  4040767.063267556]

# Non-Intersecting tile
x = 600
y = 1635
z = 12
tile_bounds = mercantile.xy_bounds(mercantile.Tile(x, y, z))
intersect_bounds(tile_bounds, merc_asset_bbox)
# None