寻找可索引GeoTIFF切片目录并按CRS坐标查询文件的Python模块

以下几个Python工具/模块可以满足你的需求，结合实际场景选择即可：

Rasterio（基础通用） #

• 核心能力：读取GeoTIFF的元数据（边界框、CRS），配合空间几何库实现自定义索引与查询
• 实现逻辑：先遍历目录提取所有TIFF的边界和CRS，存入索引结构；查询时将输入坐标/多边形转换到对应CRS，对比边界框的包含/相交关系，返回匹配文件
• 代码示例片段：

import rasterio
from rasterio.crs import CRS
from shapely.geometry import Point
from shapely.ops import transform
import pyproj
import os

# 构建本地索引
tiff_index = []
tiff_dir = "/path/to/your/tiffs"
for filename in os.listdir(tiff_dir):
    if filename.endswith(".tif"):
        filepath = os.path.join(tiff_dir, filename)
        with rasterio.open(filepath) as src:
            tiff_index.append({
                "path": filepath,
                "bounds": src.bounds,
                "crs": src.crs
            })

# 点查询函数
def query_point(target_point, point_crs):
    matched_files = []
    for tiff in tiff_index:
        # 转换查询点到当前TIFF的CRS
        proj = pyproj.Transformer.from_crs(point_crs, tiff["crs"], always_xy=True)
        transformed_pt = transform(proj.transform, target_point)
        # 判断点是否在TIFF边界内
        b = tiff["bounds"]
        if b.left <= transformed_pt.x <= b.right and b.bottom <= transformed_pt.y <= b.top:
            matched_files.append(tiff["path"])
    return matched_files

# 使用示例：查询WGS84坐标对应的TIFF
target_pt = Point(116.397, 39.908)
result = query_point(target_pt, CRS.from_epsg(4326))

GeoPandas（空间查询更便捷） #

• 核心能力：将TIFF边界转为GeoDataFrame，利用内置的空间查询API快速筛选
• 实现逻辑：遍历目录生成包含TIFF路径、边界多边形的GeoDataFrame；查询时将输入空间对象对齐CRS，通过空间连接或几何判断筛选匹配项
• 代码示例片段：

import geopandas as gpd
from shapely.geometry import box, Polygon
import rasterio
import os

# 构建GeoDataFrame索引
tiff_records = []
tiff_dir = "/path/to/your/tiffs"
for filename in os.listdir(tiff_dir):
    if filename.endswith(".tif"):
        filepath = os.path.join(tiff_dir, filename)
        with rasterio.open(filepath) as src:
            tiff_records.append({
                "path": filepath,
                "geometry": box(*src.bounds)
            })

# 统一CRS为WGS84
gdf = gpd.GeoDataFrame(tiff_records, crs=rasterio.open(filepath).crs).to_crs("EPSG:4326")

# 多边形查询函数
def query_polygon(target_poly):
    query_gdf = gpd.GeoDataFrame({"geometry": [target_poly]}, crs="EPSG:4326")
    # 空间连接筛选相交的TIFF
    matched = gpd.sjoin(gdf, query_gdf, predicate="intersects")
    return matched["path"].tolist()

# 使用示例：查询多边形覆盖的TIFF
target_poly = Polygon([(116, 39), (117, 39), (117, 40), (116, 40)])
result = query_polygon(target_poly)

TileDB（大规模切片场景） #

• 核心能力：专为大规模栅格/空间数据设计的存储引擎，支持高效的空间索引与查询
• 实现逻辑：将目录中的GeoTIFF导入TileDB数组，利用其内置空间索引直接通过坐标范围或多边形查询对应切片
• 代码示例片段：

import tiledb
import tiledb.vector_search as vs
from tiledb.crs import CRS
import os

tiff_dir = "/path/to/your/tiffs"
index_uri = "tiff_spatial_index"

# 初始化索引（仅需执行一次）
if not tiledb.array_exists(index_uri):
    vs.from_geotiffs(index_uri, tiff_dir, crs=CRS.from_epsg(4326))

# 范围查询
with tiledb.open(index_uri) as arr:
    # 查询经纬度边界内的TIFF
    query_result = arr.query(bbox=(116, 39, 117, 40))
    matched_files = query_result.meta["geotiff_paths"]

额外提示 #

• 若切片数量过万，建议结合SQLite+SpatiaLite存储索引，避免每次查询遍历所有文件
• 查询前务必确保输入坐标/多边形的CRS与索引CRS一致，或提前完成坐标转换

内容的提问来源于stack exchange，提问作者cmosig