基于gdal的面矢量面积和中心点计算（python）

基于gdal的面矢量面积和中心点计算（python）

from pathlib import Pathfrom osgeo import ogr, osrdef compute_some_metrics(inShpPath): '''计算面积''' driver = ogr.GetDriverByName("ESRI Shapefile") dataSource = driver.Open(inShpPath, 1) layer = dataSource.GetLayer() src_srs = layer.GetSpatialRef() tgt_srs = osr.SpatialReference() tgt_srs.ImportFromEPSG(32649) transform = osr.CoordinateTransformation(src_srs, tgt_srs) # geosr.SetWellKnownGeogCS("WGS_1984_UTM_Zone_49N") new_field1 = ogr.FieldDefn("Area", ogr.OFTReal) new_field1.SetWidth(32) new_field1.SetPrecision(16) layer.CreateField(new_field1) new_field2 = ogr.FieldDefn("X", ogr.OFTReal) new_field2.SetWidth(32) new_field2.SetPrecision(16) layer.CreateField(new_field2) new_field3 = ogr.FieldDefn("Y", ogr.OFTReal) new_field3.SetWidth(32) new_field3.SetPrecision(16) layer.CreateField(new_field3) for feature in layer: geom = feature.GetGeometryRef() geom2 = geom.Clone() geom2.Transform(transform) xmin, xmax, ymin, ymax = geom2.GetEnvelope() x = (xmin + xmax) / 2 y = (ymin + ymax) / 2 area_in_sq_m = geom2.GetArea() # area_in_sq_km = area_in_sq_m / 1000000 feature.SetField("Area", area_in_sq_m) layer.SetFeature(feature) feature.SetField("X", x) layer.SetFeature(feature) feature.SetField("Y", y) layer.SetFeature(feature) del dataSource

版权声明：本文内容由网络用户投稿，版权归原作者所有，本站不拥有其著作权，亦不承担相应法律责任。如果您发现本站中有涉嫌抄袭或描述失实的内容，请联系我们jiasou666@gmail.com 处理，核实后本网站将在24小时内删除侵权内容。