解决栅格数据裁剪矢量数据问题两种方法，ArcGIS解决与PYTHON解决_arcgis 使用栅格裁剪shape

需求：需要根据栅格范围裁剪矢量数据。但是查到的资料都是矢量数据裁剪栅格。
两种解决方法，arcgis/python

1、arcgis（思路栅格转矢量，先提取栅格外边界或者栅格转矢量）
下图的工具可以提取栅格的外边界数据。

下图是栅格转矢量（数据量大不建议这个工具）

转换以后就可以采用相交工具提取相交部分的矢量

有时候不知道什么原因上面两个工具转换不成功，python解决（获取的不仅有每个格子边界，还有栅格边界，获取以后再arcgis上选中导出即可同步骤1）

import osfrom osgeo import gdal, ogr, osrimport numpy as npdef raster_to_vector_boundary(input_raster_path, output_vector_path, threshold=0, band_index=1, simplify_tolerance=None): \"\"\" 将栅格数据的外边界提取为矢量格式 参数: input_raster_path (str): 输入栅格文件路径 output_vector_path (str): 输出矢量文件路径 threshold (float): 二值化阈值，大于此值的像素被视为有效区域 band_index (int): 要处理的波段索引，从1开始 simplify_tolerance (float): 边界简化容差，用于减少矢量顶点数量 \"\"\" # 打开栅格数据集 raster_ds = gdal.Open(input_raster_path) if raster_ds is None: raise FileNotFoundError(f\"无法打开栅格文件: {input_raster_path}\") # 获取栅格波段 band = raster_ds.GetRasterBand(band_index) if band is None: raise ValueError(f\"栅格中不存在索引为 {band_index} 的波段\") # 读取栅格数据为numpy数组 raster_array = band.ReadAsArray() # 创建二值化掩膜（大于阈值的像素为1，否则为0） mask_array = np.where(raster_array > threshold, 1, 0).astype(np.uint8) # 获取栅格投影和地理变换信息 projection = raster_ds.GetProjection() geotransform = raster_ds.GetGeoTransform() # 创建临时内存栅格用于存储掩膜 driver = gdal.GetDriverByName(\'MEM\') mask_ds = driver.Create(\'\', raster_ds.RasterXSize, raster_ds.RasterYSize, 1, gdal.GDT_Byte) mask_ds.SetProjection(projection) mask_ds.SetGeoTransform(geotransform) mask_band = mask_ds.GetRasterBand(1) mask_band.WriteArray(mask_array) mask_band.FlushCache() # 创建输出矢量数据集 vector_driver = ogr.GetDriverByName(\'GeoJSON\') if output_vector_path.lower().endswith( \'.geojson\') else ogr.GetDriverByName(\'ESRI Shapefile\') # 如果输出文件已存在，删除它 if os.path.exists(output_vector_path): vector_driver.DeleteDataSource(output_vector_path) vector_ds = vector_driver.CreateDataSource(output_vector_path) if vector_ds is None: raise RuntimeError(f\"无法创建矢量文件: {output_vector_path}\") # 创建图层 srs = osr.SpatialReference() srs.ImportFromWkt(projection) layer = vector_ds.CreateLayer(\'boundary\', srs, ogr.wkbPolygon) # 创建字段用于存储面积 area_field = ogr.FieldDefn(\'area\', ogr.OFTReal) layer.CreateField(area_field) # 执行栅格到矢量的转换 gdal.Polygonize(mask_band, None, layer, 0, [], callback=None) # 遍历所有要素并简化边界（如果指定了容差） if simplify_tolerance is not None and simplify_tolerance > 0: for feature in layer: geom = feature.GetGeometryRef() simplified_geom = geom.SimplifyPreserveTopology(simplify_tolerance) feature.SetGeometry(simplified_geom) layer.SetFeature(feature) # 计算并设置面积字段 for feature in layer: geom = feature.GetGeometryRef() area = geom.GetArea() feature.SetField(\'area\', area) layer.SetFeature(feature) # 释放资源 raster_ds = None mask_ds = None vector_ds = None print(f\"成功将栅格外边界提取为矢量数据并保存至: {output_vector_path}\")if __name__ == \"__main__\": # 栅格文件路径（根据你的提供设置） input_raster_path = r\"G:\\**\\***\\result.tif\" # 输出矢量文件路径（默认为栅格所在目录下的boundary.shp） output_dir = r\"G:\\***\\****\" output_vector_path = os.path.join(output_dir, \"boundary.shp\") # 处理参数 threshold = 0 # 二值化阈值，根据实际情况调整 band_index = 1 # 处理第一个波段 simplify_tolerance = 10 # 边界简化容差，可根据需要调整或设为None不简化 try: # 执行栅格到矢量的转换 raster_to_vector_boundary( input_raster_path, output_vector_path, threshold=threshold, band_index=band_index, simplify_tolerance=simplify_tolerance ) except Exception as e: print(f\"处理过程中发生错误: {str(e)}\")