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

from pathlib import Path
from osgeo import ogr, osr

def 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