当前位置：首页 > news >正文

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

news 2025/7/23 11:48:44

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

1、arcgis（思路栅格转矢量，先提取栅格外边界或者栅格转矢量）
下图的工具可以提取栅格的外边界数据。
在这里插入图片描述
下图是栅格转矢量（数据量大不建议这个工具）

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

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

import os
from osgeo import gdal, ogr, osr
import 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 = Nonemask_ds = Nonevector_ds = Noneprint(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)}")