python计算shp中每个区域的面积

简介

如果我们拿到一个区域的shp后想对其增加一个面积字段并写入对应小区域的面积，可以在arcgis中用字段计算器完成，但是步骤相对比较繁琐，而且如果需要批量处理多个shp的话会显得更加不方便。这里提供一个python代码可以一站式解决这个问题。

具体代码

from osgeo import ogr, osr
import os# 输入shapefile路径
input_shp = r"D:\your.shp"# 打开数据源
driver = ogr.GetDriverByName("ESRI Shapefile")
datasource = driver.Open(input_shp, 1)  # 1表示可写入if datasource is None:print("无法打开输入的shapefile文件")exit(1)# 获取图层
layer = datasource.GetLayer()# 获取图层的空间参考
source_srs = layer.GetSpatialRef()
print(f"源坐标系统: {source_srs.ExportToWkt()}")# 创建等面积投影坐标系统（适合面积计算）
# 使用Lambert Azimuthal Equal Area投影，中心点设置为数据范围的中心
# 获取数据的范围
extent = layer.GetExtent()
x_center = (extent[0] + extent[1]) / 2  # (minX + maxX) / 2
y_center = (extent[2] + extent[3]) / 2  # (minY + maxY) / 2target_srs = osr.SpatialReference()
target_srs.SetProjCS("Lambert Azimuthal Equal Area")
target_srs.SetWellKnownGeogCS("WGS84")
target_srs.SetLAEA(y_center, x_center, 0, 0)  # 纬度, 经度, 假东, 假北# 创建坐标转换
transform = osr.CoordinateTransformation(source_srs, target_srs)# 检查字段是否存在
layer_defn = layer.GetLayerDefn()
field_index = -1
for i in range(layer_defn.GetFieldCount()):if layer_defn.GetFieldDefn(i).GetName() == "AREA_KM2":field_index = ibreak# 如果字段不存在，则添加
if field_index == -1:field_defn = ogr.FieldDefn("AREA_KM2", ogr.OFTReal)field_defn.SetWidth(12)field_defn.SetPrecision(4)layer.CreateField(field_defn)# 遍历所有要素，计算面积
feature_count = layer.GetFeatureCount()
print(f"开始计算 {feature_count} 个要素的面积...")for i in range(feature_count):feature = layer.GetFeature(i)geometry = feature.GetGeometryRef()# 克隆几何对象并转换到等面积投影projected_geometry = geometry.Clone()projected_geometry.Transform(transform)# 计算面积（平方米）并转换为平方公里area_m2 = projected_geometry.GetArea()area_km2 = area_m2 / 1000000.0# 更新字段feature.SetField("AREA_KM2", area_km2)layer.SetFeature(feature)# 清理feature = Nonegeometry = Noneprojected_geometry = None# 显示进度if (i + 1) % 100 == 0 or i + 1 == feature_count:print(f"已处理: {i + 1}/{feature_count}")# 保存并关闭
datasource = Noneprint("面积计算完成，字段 'AREA_KM2' 已添加，单位为平方公里。")