python：在多景遥感影像中根据shp的范围找到所覆盖的遥感影像

这段时间有一个这种需求：给你一个shp文件，让你在一个有很多幅遥感影像的文件夹里面快速找到与shp相交的那几幅影像。

方法如下：

1、首先我们需要读取shp文件的经纬度范围：

from osgeo import ogr

# 打开SHP文件
dataSource = ogr.Open('path_to_your_shapefile.shp', 1)  # 1表示以读写模式打开
layer = dataSource.GetLayer() #获取SHP文件的图层（Layer）对象
shpminx, shpminy, shpmaxx, shpmaxy = layer.GetExtent() 

dataSource.Destroy() #关闭数据源

2、然后获取遥感影像tif的经纬度范围：

from osgeo import gdal

def get_geo_bounds(tif_path):
    # 打开TIF文件
    dataset = gdal.Open(tif_path, gdal.GA_ReadOnly)
    
    # 获取地理变换参数(这个参数详解在下面)
    geo_transform = dataset.GetGeoTransform()
    
    # 计算经纬度范围
    tifminx = geo_transform[0]
    tifmaxy = geo_transform[3] + geo_transform[5] * dataset.RasterYSize #起始纬度+分辨率*行数
    tifmaxx = minx + geo_transform[1] * dataset.RasterXSize
    tifminy = maxy - geo_transform[5] * dataset.RasterYSize
    
    # 关闭数据集
    dataset = None
    
    return (tifminx, tifminy, tifmaxx, tifmaxy)

gdal地理变换参数是什么？

地理变换参数通常包含六个值，分别表示：

• GT(0): 左上角像素的X坐标（经度）。
• GT(1): 东西方向上的像素分辨率。
• GT(2): 像素在东西方向上的旋转参数，在大多数情况下，这个值是0。
• GT(3): 左上角像素的Y坐标（纬度）。
• GT(4): 南北方向上的旋转参数，通常为0。
• GT(5): 北南方向上的像素分辨率（对于北半球的图像通常是负值）。

3、根据shp范围判断有哪些遥感影像与shp相交（可以覆盖shp）

# 存储符合条件的影像路径的列表
selected_images = []

folder_path = '包含TIF文件的文件夹路径'
# 获取文件夹中所有的TIF文件
tif_files = [f for f in os.listdir(folder_path) if f.lower().endswith('.tif')]
# 遍历TIF文件列表并获取每个文件的经纬度范围
for tif_file in tif_files:
    tif_path = os.path.join(folder_path, tif_file)
    tifminx, tifminy, tifmaxx, tifmaxy = get_geo_bounds(tif_path)
    # 判断GeoJSON文件的范围是否与当前遥感影像的范围有交集
    # 这一块是代码的核心，有点绕，可以仔细思考一下
    if not (shpminx > tifmaxx or  # shp的左边界在影像右边界的右侧
        shpmaxx < tifminx or      # shp的右边界在影像左边界的左侧
        shpminy > tifmaxy or      # shp的下边界在影像上边界的上侧
        shpmaxy < shpminy):       # shp的上边界在影像下边界的下侧
    selected_images.append(tif_file)
print(selected_images) #输出结果

完整代码：

from osgeo import ogr
from osgeo import gdal


# 打开SHP文件
dataSource = ogr.Open('path_to_your_shapefile.shp', 1)  # 1表示以读写模式打开
layer = dataSource.GetLayer() #获取SHP文件的图层（Layer）对象
shpminx, shpminy, shpmaxx, shpmaxy = layer.GetExtent() 

dataSource.Destroy() #关闭数据源



def get_geo_bounds(tif_path):
    # 打开TIF文件
    dataset = gdal.Open(tif_path, gdal.GA_ReadOnly)
    
    # 获取地理变换参数(这个参数详解在下面)
    geo_transform = dataset.GetGeoTransform()
    
    # 计算经纬度范围
    tifminx = geo_transform[0]
    tifmaxy = geo_transform[3] + geo_transform[5] * dataset.RasterYSize #起始纬度+分辨率*行数
    tifmaxx = minx + geo_transform[1] * dataset.RasterXSize
    tifminy = maxy - geo_transform[5] * dataset.RasterYSize
    
    # 关闭数据集
    dataset = None
    
    return (tifminx, tifminy, tifmaxx, tifmaxy)

# 存储符合条件的影像路径的列表
selected_images = []

folder_path = '包含TIF文件的文件夹路径'
# 获取文件夹中所有的TIF文件
tif_files = [f for f in os.listdir(folder_path) if f.lower().endswith('.tif')]
# 遍历TIF文件列表并获取每个文件的经纬度范围
for tif_file in tif_files:
    tif_path = os.path.join(folder_path, tif_file)
    tifminx, tifminy, tifmaxx, tifmaxy = get_geo_bounds(tif_path)
    # 判断GeoJSON文件的范围是否与当前遥感影像的范围有交集
    # 这一块是代码的核心，有点绕，可以仔细思考一下
    if not (shpminx > tifmaxx or  # shp的左边界在影像右边界的右侧
        shpmaxx < tifminx or      # shp的右边界在影像左边界的左侧
        shpminy > tifmaxy or      # shp的下边界在影像上边界的上侧
        shpmaxy < shpminy):       # shp的上边界在影像下边界的下侧
    selected_images.append(tif_file)
print(selected_images) #输出结果