Calculate which city points are within each tif in a directory

index.py

# this script goes through a directory of tif files and a shapefile of cities
# and reports which city points are contained within each tif

import os
import fiona
import rasterio
import shapely
from pyproj import Transformer

tifpath = "path/to/tiffs"

d = os.listdir(tifpath)
for filename in d:
    if not filename.endswith(".tif"):
        continue

    image = rasterio.open(os.path.join(tifpath, filename))
    rasterextent = shapely.geometry.box(*image.bounds)

    # the city points are in coordinate system EPSG:4269 (latitude, longitude)
    # so we'll need to transform into the CRS of the raster
    t = Transformer.from_crs(4269, image.crs, always_xy=True)

    with fiona.open("cities.shp") as cities:
        for feature in cities:
            # transform the city coordinates to match the CRS of the current raster
            x, y = feature.geometry.coordinates
            newx, newy = t.transform(x, y)
            newpoint = shapely.geometry.Point(newx, newy)

            if shapely.contains(rasterextent, newpoint):
                # current city matches the raster
                print(filename, feature.properties["NAME10"])