Extract the bounding envelope from an ESRI Shapefile using OGR. The extents represent the minimum and maximum longitudes and latitudes that enclose all features.

from osgeo import ogr

vector_source = ogr.Open("area_of_interest.shp", 1)
feature_layer = vector_source.GetLayer()
left, bottom, right, top = feature_layer.GetExtent()
vector_source.Destroy()

A utility function reads a GeoTIFF and derives its coverage bounds from the affine geotransform and raster dimensions. The geotransform encodes the origin, pixel size, and rotation.

GT[0] – X coordinate of the upper left corner
GT[1] – pixel width (east-west resolusion)
GT[2] – rotation about X axis (usually 0)
GT[3] – Y coordinate of the upper left corner
GT[4] – rotation about Y axis (usual 0)
GT[5] – pixel height (north-south resolution, often negative for north-up images)

from osgeo import gdal

def raster_bbox(image_path):
    handle = gdal.Open(image_path, gdal.GA_ReadOnly)
    gt = handle.GetGeoTransform()
    cols = handle.RasterXSize
    rows = handle.RasterYSize

    west = gt[0]
    north = gt[3]
    east = west + gt[1] * cols
    south = north + gt[5] * rows

    handle = None
    return west, south, east, north

With both extents38known, intersection can be tested by checking weather rectangles overlap. The logic reverses the condition for no overlap: one shape lies completely to the left, right, above, or below the other.

import os

matches = []
directory = "./satellite_scenes"

candidates = [
    entry for entry in os.listdir(directory)
    if entry.lower().endswith(".tif")
]

for filename in candidates:
    path = os.path.join(directory, filename)
    w, s, e, n = raster_bbox(path)

    disjoint = (
        left > e or
        right < w or
        bottom > n or
        top < s
    )
    if not disjoint:
        matches.append(filename)

print(matches)

The full script brings together reading the shapefile, iterating over TIF files, computing bounding boxes, and collecting intersecting images.

from osgeo import ogr, gdal
import os

data_source = ogr.Open("area_of_interest.shp", 1)
layer = data_source.GetLayer()
min_x, min_y, max_x, max_y = layer.GetExtent()
data_source.Destroy()

def get_raster_extent(tif_path):
    dataset = gdal.Open(tif_path, gdal.GA_ReadOnly)
    transform = dataset.GetGeoTransform()
    cols = dataset.RasterXSize
    rows = dataset.RasterYSize
    west = transform[0]
    north = transform[3]
    east = west + transform[1] * cols
    south = north + transform[5] * rows
    dataset = None
    return west, south, east, north

overlapping_scenes = []
folder = "./imagery"

for file in os.listdir(folder):
    if not file.lower().endswith(".tif"):
        continue
    tif_file = os.path.join(folder, file)
    west, south, east, north = get_raster_extent(tif_file)

    # overlap test
    if not (min_x > east or max_x < west or
            min_y > north or max_y < south):
        overlapping_scenes.append(file)

print(overlapping_scenes)