R语言中按多边形提取对应阈值以上的栅格像素值

针对你这个需要给38个州（多边形）设置不同阈值，提取栅格中符合条件的像素（设为1，其余为NoData）的需求，我整理了几个实用的方案，你可以根据自己的工具栈选择：

方案一：ArcPy（ArcGIS Python脚本） #

如果你有ArcGIS Desktop/Pro的许可，用ArcPy写脚本可以高效实现自动化处理：

前置准备 #

确保你的shapefile包含一个数值型字段（比如命名为threshold），每个州对应的阈值已经准确填入该字段。

实现代码 #

import arcpy
import os

# 配置工作环境
arcpy.env.workspace = r"C:\Your_Workspace_Path"
arcpy.env.overwriteOutput = True  # 允许覆盖已有输出

# 定义输入输出路径
input_raster = "your_input_raster.tif"
state_shapefile = "states.shp"
output_tiles_folder = r"C:\Your_Workspace_Path\state_results"

# 创建输出文件夹（不存在则新建）
if not os.path.exists(output_tiles_folder):
    os.makedirs(output_tiles_folder)

# 遍历每个州的多边形要素
with arcpy.da.SearchCursor(state_shapefile, ["OID@", "threshold", "SHAPE@"]) as cursor:
    for row in cursor:
        oid = row[0]
        threshold_val = row[1]
        state_polygon = row[2]
        
        # 1. 提取当前州范围内的栅格子集
        masked_raster = arcpy.sa.ExtractByMask(input_raster, state_polygon)
        
        # 2. 按阈值重分类：大于阈值设为1，其余为NoData
        reclassified_raster = arcpy.sa.Con(masked_raster > threshold_val, 1)
        
        # 3. 保存当前州的结果
        output_raster_path = os.path.join(output_tiles_folder, f"state_{oid}_result.tif")
        reclassified_raster.save(output_raster_path)

# 可选：将所有州的结果合并为单个栅格文件
arcpy.MosaicToNewRaster_management(
    input_rasters=[os.path.join(output_tiles_folder, f) for f in os.listdir(output_tiles_folder) if f.endswith(".tif")],
    output_location=arcpy.env.workspace,
    raster_dataset_name_with_extension="final_combined_result.tif",
    pixel_type="1_BIT",  # 仅存1和NoData，用1位类型节省存储空间
    number_of_bands=1
)

注意事项 #

• 确保threshold字段是整数/浮点型，不要是文本型；
• 运行脚本前需激活ArcPy环境（比如在ArcGIS Pro的Python窗口运行，或配置好ArcGIS的Python环境变量）。

方案二：GDAL + Python（开源无许可） #

如果你没有ArcGIS许可，用GDAL这个开源地理数据处理库完全可以实现相同需求：

前置准备 #

• 安装GDAL库：pip install gdal（注意对应你的系统和Python版本，若安装失败可尝试用conda安装：conda install -c conda-forge gdal）；
• 同样确保shapefile带有threshold数值字段。

实现代码 #

from osgeo import gdal, ogr
import numpy as np
import os

# 注册GDAL/OGR驱动
gdal.AllRegister()
ogr.RegisterAll()

# 定义输入输出路径
input_raster_path = "your_input_raster.tif"
state_shapefile_path = "states.shp"
output_tiles_folder = "state_results"

# 创建输出文件夹
os.makedirs(output_tiles_folder, exist_ok=True)

# 打开栅格和矢量数据源
raster_ds = gdal.Open(input_raster_path)
shape_ds = ogr.Open(state_shapefile_path)
state_layer = shape_ds.GetLayer()

# 获取栅格的地理变换和投影信息
geo_transform = raster_ds.GetGeoTransform()
projection = raster_ds.GetProjection()

# 遍历每个州的要素
for state_feature in state_layer:
    fid = state_feature.GetFID()
    threshold_val = state_feature.GetField("threshold")
    state_geom = state_feature.GetGeometryRef()
    
    # 创建临时矢量文件（GDAL Warp需要输入矢量文件而非单个要素）
    temp_shp_path = f"temp_state_{fid}.shp"
    temp_ds = ogr.GetDriverByName("ESRI Shapefile").CreateDataSource(temp_shp_path)
    temp_layer = temp_ds.CreateLayer("temp_state", srs=state_layer.GetSpatialRef())
    temp_layer.CreateFeature(state_feature.Clone())
    temp_ds = None  # 关闭临时文件释放资源
    
    # 1. 裁剪栅格到当前州范围
    masked_raster_path = f"temp_masked_{fid}.tif"
    gdal.Warp(
        masked_raster_path,
        raster_ds,
        cutlineDSName=temp_shp_path,
        cropToCutline=True,
        dstNodata=np.nan
    )
    
    # 读取裁剪后的栅格数据
    masked_ds = gdal.Open(masked_raster_path)
    raster_array = masked_ds.ReadAsArray()
    
    # 2. 按阈值重分类：大于阈值设为1，其余为NaN（对应NoData）
    result_array = np.where(raster_array > threshold_val, 1, np.nan)
    
    # 3. 保存结果栅格
    output_raster_path = os.path.join(output_tiles_folder, f"state_{fid}_result.tif")
    driver = gdal.GetDriverByName("GTiff")
    output_ds = driver.Create(
        output_raster_path,
        masked_ds.RasterXSize,
        masked_ds.RasterYSize,
        1,
        gdal.GDT_Byte,  # 用Byte类型存储，节省空间
        options=["COMPRESS=LZW"]
    )
    output_ds.SetGeoTransform(masked_ds.GetGeoTransform())
    output_ds.SetProjection(masked_ds.GetProjection())
    output_band = output_ds.GetRasterBand(1)
    output_band.SetNoDataValue(np.nan)
    output_band.WriteArray(result_array)
    
    # 释放资源并清理临时文件
    output_ds = None
    masked_ds = None
    for ext in [".shp", ".dbf", ".shx", ".prj"]:
        temp_file = temp_shp_path.replace(".shp", ext)
        if os.path.exists(temp_file):
            os.remove(temp_file)
    os.remove(masked_raster_path)

# 可选：合并所有结果为单个栅格
vrt_path = "combined_result.vrt"
gdal.BuildVRT(vrt_path, [os.path.join(output_tiles_folder, f) for f in os.listdir(output_tiles_folder) if f.endswith(".tif")])
gdal.Translate("final_combined_result.tif", vrt_path, options=["COMPRESS=LZW"])

注意事项 #

• 临时文件会自动清理，若脚本中断可能残留临时文件，手动删除即可；
• 若栅格数据较大，建议分块处理（GDAL默认会分块，无需额外设置）。

方案三：QGIS图形界面操作（无需编程） #

如果你不想写代码，用QGIS的图形界面也能完成批量处理：

1. 加载数据：将shapefile和栅格数据导入QGIS；
2. 拆分矢量图层：打开「Processing Toolbox」→「Vector general」→「Split vector layer」，选择shapefile，按OID或州名称字段拆分，将38个州的多边形保存到单独文件夹；
3. 批量裁剪栅格：打开「Processing Toolbox」→「GDAL」→「Raster extraction」→「Clip raster by mask layer」，点击「Run as batch process」：
   • 输入栅格选择原始栅格；
   • 掩膜图层选择刚才拆分的所有多边形文件；
   • 输出文件夹指定一个位置；
4. 批量重分类：打开「Processing Toolbox」→「GDAL」→「Raster calculator」，点击「Run as batch process」：
   • 对每个裁剪后的栅格，设置表达式为("输入栅格名称@1" > 对应阈值) * 1；
   • 可以通过「Autofill」功能批量匹配每个栅格对应的阈值（从拆分的多边形属性中提取）；
5. 合并结果：打开「Processing Toolbox」→「GDAL」→「Raster miscellaneous」→「Merge」，选择所有重分类后的栅格，合并为单个输出文件。

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