任务：

对面要素shp文件进行处理，获取每个要素的中心点，（判断中心点在不在面上），存放到shp的属性表的新建字段中。在根据这些点坐标利用geopandas库生成点shp文件。

代码部分：

对类进行初始化操作

from osgeo import ogr,osr
import os
import re
import csv
import time
import numpy as np
import pandas as pd
import geopandas as gpd
from shapely.geometry import Point

class rivers:
    def __init__(self,shp_path):
        driver = ogr.GetDriverByName("ESRI Shapefile")
        self.datasource=driver.Open(shp_path,1)
        self.cen_layer=self.datasource.GetLayer()
        self.feature_defn=self.cen_layer.GetLayerDefn()#获取图层定义
        self.feature_count=self.feature_defn.GetFieldCount()#获取属性表中字段个数
        self.Num_feature=self.cen_layer.GetFeatureCount()
        #生成的txt,csv,shp文件的路径
        self.txt_path = shp_path[0:-4] + ".txt"
        self.csv_path = shp_path[0:-4] + ".csv"
        self.shp_path = shp_path[0:-4] + "_point.shp"

        # for i in range(self.feature_count):
        #     field_defn=self.feature_defn.GetFieldDefn(i)
        #     print("字段名："+field_defn.GetNameRef()+'\t'+"字段类型："+field_defn.GetFieldTypeName(field_defn.GetType())+'\n')#遍历所有字段其类型

创建字段：

    def create_field(self,field_name):
        fieldIndex = self.feature_defn.GetFieldIndex(field_name)
        if fieldIndex < 0:
            # 添加字段
            fieldDefn = ogr.FieldDefn(field_name, ogr.OFTReal)#添加字段名及设置其类型,Real为浮点类型
            fieldDefn.SetPrecision(12)#设置精度
            self.cen_layer.CreateField(fieldDefn, 1)
            fieldIndex2 = self.feature_defn.GetFieldIndex(field_name)
            if fieldIndex2 > 0:
                print("字段创建成功：", field_name)
        else:
            print(field_name+"字段已存在，无需创建！")

注意：

fieldDefn = ogr.FieldDefn(field_name, ogr.OFTReal)#添加字段名及设置其类型,Real为浮点类型 

这一句中，OGR创建字段的类型有以下几种：

详情见链接[Python] GDAL/OGR操作矢量数据（shp、GeoJSON）_GeoDoer-CSDN博客 

 将中点坐标写入属性表：

    def write_center(self):
        for feature in self.cen_layer:
            geometry = feature.GetGeometryRef()
            extent = geometry.GetEnvelope()
            po_x=(extent[0]+extent[1])/2
            po_y = (extent[2] + extent[3]) / 2
            # print("X坐标："+str(po_x)+"  Y坐标："+str(po_y))
            feature.SetField('center_X', float(po_x))
            feature.SetField('center_Y', float(po_y))
            self.cen_layer.SetFeature(feature)

 注意：

feature.SetField('center_X', float(po_x)) #对字段进行赋值

  对字段进行赋值之后，一定要将其写入属性表，也可以理解为更新到属性表中：

self.cen_layer.SetFeature(feature) #写入属性表

检查部分：

检查点是否在面内

features.Contains(point)

注意point的生成方式以及features的表示内容：

wkt = "POINT(%f %f)" % (float(Point_X),float(Point_Y))  ## 创建点
point = ogr.CreateGeometryFromWkt(wkt)  ## 生成点
features=feature.GetGeometryRef()#获取Geometry属性

如果不在，则使用

PointOnSurface()

方法找到一定在面内的一点，再将该面要素对应的属性表数据修改为点坐标

 def check_center(self):
        count=0
        lake_num=0
        id=[]
        for feature in self.cen_layer:
            FID = int(feature.GetFID())
            Point_X = feature.GetFieldAsDouble('center_X')
            Point_Y = feature.GetFieldAsDouble('center_Y')
            wkt = "POINT(%f %f)" % (float(Point_X),float(Point_Y))  ## 创建点
            point = ogr.CreateGeometryFromWkt(wkt)  ## 生成点
            features=feature.GetGeometryRef()#获取Geometry属性

            if not features.Contains(point):#判断点是否在面上
                point_t = features.PointOnSurface()#获取面内的一点坐标
                data = str(point_t).split(' ', 1)[1]
                po_x = re.findall("[(](.*?)[ ]", data)
                po_y = re.findall("[ ](.*?)[)]", data)
                # print("X坐标：" + str(po_x[0]) + "  Y坐标：" + str(po_y[0]))
                feature.SetField('center_X', float(po_x[0]))
                feature.SetField('center_Y', float(po_y[0]))
                self.cen_layer.SetFeature(feature)#将修改的字段数据写入图层
                id.append(FID)
                count+=1
            area=feature.GetFieldAsDouble('Lake_area')
            if area<100:
                lake_num+=1
                cen_layer.DeleteFeature(FID)
        print("总共有"+str(self.Num_feature)+"个湖库！")
        print("有" + str(count) + "个湖泊中心点不在面内！")
        if id:
            print("中心点不在面内的面ID为："+id)
        print("面积小于100平方公里的湖库有"+str(lake_num)+"个！")

根据属性表中的字段值生成txt文件 

 

    def GetPoint(self,shp_path):
        # driver = ogr.GetDriverByName("ESRI Shapefile")
        # datasource = driver.Open(shp_path, 1)
        # cen_layer = datasource.GetLayer()
        # self.feature_defn = self.cen_layer.GetLayerDefn()  # 获取图层定义
        # feature_count = self.feature_defn.GetFieldCount()  # 获取属性表中字段个数
        with open(self.txt_path,'w') as f:
            for feature in self.cen_layer:
                Point_X = feature.GetFieldAsDouble('center_X')#以double的数据格式读取字段值
                Point_Y = feature.GetFieldAsDouble('center_Y')
                f.write(str(Point_X) + "," + str(Point_Y)+'\n')
        f.close()

txt转为csv

 def txt_to_csv(self):
        fh = open(self.csv_path, "w+", newline='')
        writer = csv.writer(fh)
        writer.writerow(["lon", "lat"])#列名
        file = open(self.txt_path, encoding='utf-8')
        data = file.readlines()
        res = []
        for i in range(len(data)):
            points = data[i].strip().split(',')
            d = [points[0], points[1]]
            res.append(d)
        writer.writerows(res)
        file.close()
        fh.close()

通过csv生成 point点矢量

geopandas生成shp数据，详情参考Python GeoPandas 文本经纬度转换为点要素、线要素 - 简书

    def csv_to_points(self):
        df = pd.read_csv(self.csv_path, header=0, encoding='gbk')
        geometry=[Point(xy) for xy in zip(df.lon, df.lat)]#以Point格式读取csv
        gdf = gpd.GeoDataFrame(df, crs="EPSG:4326",geometry=geometry)#将其坐标系设为WGS_84
        gdf.to_file(self.shp_path, encoding='gbk')#将其写入shp文件中

 主函数

 

if __name__ == '__main__':
    start=time.time()
    shp_path=r'C:\Users\Administrator\Desktop\rivers\HydroLAKES_polys_v10_shp\HydroLAKES_polys_v10.shp'
    river_shp=rivers(shp_path)#实例化对象
    field_names=['center_X','center_Y']#需要创建的字段名列表
    for field_name in field_names:
        river_shp.create_field(field_name)
    # river_shp.write_center()#将中心点坐标填充进字段表
    river_shp.check_center()#检查中心点在不在面内，并进行修正
    river_shp.GetPoint(shp_path)#生成保存中心点坐标的txt文件
    river_shp.txt_to_csv()#将txt文件转为csv
    river_shp.csv_to_points()#根据csv生成对应点要素shp文件
    end=time.time()
    print("用时:"+str(end-start)+"s")