如何用Python从葡萄牙地形坡度TIFF中筛选特定坡度区域并生成Shapefile？

1. 安装依赖库 #

先确保安装所需工具库：

pip install rasterio geopandas matplotlib numpy

2. 核心代码实现 #

import rasterio
import geopandas as gpd
import numpy as np
import matplotlib.pyplot as plt
from rasterio.mask import mask
from rasterio.features import shapes

# 读取葡萄牙Shapefile和坡度TIFF
portugal_shp = gpd.read_file("path/to/portugal.shp")
slope_tiff = rasterio.open("path/to/slope.tif")

# 确保CRS一致，若不一致则转换Shapefile的CRS匹配TIFF
if portugal_shp.crs != slope_tiff.crs:
    portugal_shp = portugal_shp.to_crs(slope_tiff.crs)

# 裁剪TIFF到葡萄牙边界范围（减少计算量）
cropped_slope, cropped_transform = mask(
    slope_tiff,
    portugal_shp.geometry,
    crop=True,
    nodata=slope_tiff.nodata
)

# 创建坡度筛选掩码（筛选坡度<20%的区域）
slope_data = cropped_slope[0]  # 取单波段数据（坡度TIFF通常为单波段）
mask = (slope_data < 20) & (slope_data != slope_tiff.nodata)  # 排除无数据区域

# 将掩码转换为矢量多边形
results = (
    {"geometry": shape, "properties": {"slope_class": "<20%"}}
    for shape, value in shapes(slope_data, mask=mask, transform=cropped_transform)
    if value == 1
)

# 转为GeoDataFrame
filtered_gdf = gpd.GeoDataFrame.from_features(list(results), crs=slope_tiff.crs)

# 只保留葡萄牙境内的区域（避免裁剪后仍有境外无效像素）
filtered_gdf = gpd.overlay(filtered_gdf, portugal_shp, how="intersection")

# 绘制结果
fig, ax = plt.subplots(figsize=(10, 10))
portugal_shp.plot(ax=ax, facecolor="white", edgecolor="black", linewidth=1)
filtered_gdf.plot(ax=ax, facecolor="#4CAF50", alpha=0.7)
plt.title("葡萄牙坡度小于20%的区域")
plt.axis("off")
plt.show()

# 可选：保存筛选后的Shapefile
filtered_gdf.to_file("portugal_slope_less_20.shp")

关键注意事项 #

• CRS匹配：必须保证Shapefile和TIFF的坐标参考系统一致，否则裁剪、交集操作会出错。
• 无数据处理：一定要排除TIFF中的无数据（nodata）区域，避免生成无效多边形。
• 掩码逻辑：可根据需求调整条件，比如slope_data > 20筛选坡度大于20%的区域。
• 性能优化：若TIFF文件过大，先裁剪到Shapefile范围再处理，能大幅减少内存占用。

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