如何判断给定坐标是否在TIFF栅格文件内并提取N×N窗口？

好的，我来帮你一步步解决这两个问题——先判断示例坐标是否在栅格范围内，再提取指定大小的窗口：

第一步：判断坐标是否在栅格内 #

首先要注意：你的示例坐标是WGS84经纬度（EPSG:4326），但栅格文件用的是EPSG:31370（比利时兰伯特投影），两种坐标系不匹配，所以必须先把经纬度坐标转换成栅格的投影坐标系，再和栅格边界对比。

这里我们用rasterio配合pyproj来完成坐标转换，代码如下：

import rasterio
from pyproj import Transformer

# 打开栅格文件
src = rasterio.open('k_01.tif')

# 创建坐标转换器：从WGS84（EPSG:4326）转栅格的EPSG:31370
# always_xy=True 确保输入顺序是 (经度, 纬度)，和坐标系的x/y对应
transformer = Transformer.from_crs("EPSG:4326", src.crs, always_xy=True)

# 你的示例坐标：纬度51.3334198，经度3.2973934
lon, lat = 3.2973934, 51.3334198
# 转换为栅格坐标系下的(x, y)
x, y = transformer.transform(lon, lat)

# 获取栅格的边界范围
bbox = src.bounds

# 判断转换后的坐标是否在栅格内部
is_inside = (bbox.left <= x <= bbox.right) and (bbox.bottom <= y <= bbox.top)
print(f"坐标是否在栅格内：{is_inside}")

运行这段代码后，就能直观知道你的坐标是否落在k_01.tif的范围内了。

第二步：提取N×N窗口（若坐标在栅格内） #

如果坐标确实在栅格里，接下来需要把地理坐标转换成栅格的行列号，然后计算N×N窗口的范围，最后读取窗口数据。需要注意处理边界情况——比如当坐标靠近栅格边缘时，窗口不能超出栅格的行列范围，否则会报错。

完整代码示例：

if is_inside:
    # 定义你需要的窗口大小N
    N = 5  # 可以根据需求修改，比如改成10、20等

    # 将栅格坐标系下的(x, y)转换为像素行列号（注意：返回的是(col, row)）
    col, row = src.index(x, y)

    # 计算窗口的起始/结束行列
    half_N = N // 2
    # 确保起始行列不小于0，结束行列不超过栅格的总高/宽
    row_start = max(0, row - half_N)
    row_end = min(src.height, row + half_N + 1)  # 左闭右开，所以要+1
    col_start = max(0, col - half_N)
    col_end = min(src.width, col + half_N + 1)

    # 创建窗口对象并读取数据
    window = rasterio.windows.Window.from_slices((row_start, row_end), (col_start, col_end))
    window_data = src.read(window=window)

    # 输出结果：window_data的形状是 (波段数, 窗口高度, 窗口宽度)
    print(f"提取的窗口数据形状：{window_data.shape}")
else:
    print("坐标不在栅格范围内，无法提取窗口")

补充小提示 #

• 如果你的栅格是多波段影像（比如RGB），window_data的第一个维度是波段数，比如3波段的话，形状会是(3, N, N)（如果窗口没被边缘裁剪）。
• 如果你只需要读取单个波段，可以用src.read(1, window=window)（1代表第一个波段）。
• 要是还没安装pyproj，可以用pip install pyproj快速安装。

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