对面要素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)
注意：
①：利用外接矩形获取中心点：
使用GetEnvelope()获取外接矩形时（结果是包含Xmax,Xmin,Ymax,Ymin的列表），一定要先获取feature的空间属性：
geometry = feature.GetGeometryRef() extent = geometry.GetEnvelope()
②：对字段进行赋值：
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:UsersAdministratorDesktopriversHydroLAKES_polys_v10_shpHydroLAKES_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")
