学习日记：获取属性表值并修改（判断点是否在面内）

需求：

判断arcgis的字段计算器所计算的面要素的质心是否在面内，如果不在，则用PointOnSurface()方法生成一个在面内的点，并读取该点的经纬度，将其定义为属性表内的center_X,center_Y两个字段的值。

重点1：features.Contains(point)，判断点是否在面内

重点2：   wkt = "POINT(%f %f)" % (Point_X,Point_Y)  ## 创建点
                point = ogr.CreateGeometryFromWkt(wkt)  ## 生成点

from osgeo import ogr,osr
import os,re
import numpy as np
def cal_center(shp_path):
driver = ogr.GetDriverByName("ESRI Shapefile")
datasource=driver.Open(shp_path,1)
cen_layer=datasource.GetLayer()
feature_defn=cen_layer.GetLayerDefn()#获取图层定义
feature_count=feature_defn.GetFieldCount()#获取属性表中字段个数
Num_feature=cen_layer.GetFeatureCount()
# for i in range(feature_count):
#
field_defn=feature_defn.GetFieldDefn(i)
#
print("字段名："+field_defn.GetNameRef()+'t'+"字段类型："+field_defn.GetFieldTypeName(field_defn.GetType())+'n')
count=0
id=[]
for feature in cen_layer:
FID = int(feature.GetFID())
Point_X = feature.GetFieldAsDouble('center_X')
Point_Y = feature.GetFieldAsDouble('center_Y')
wkt = "POINT(%f %f)" % (Point_X,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)
feature.SetField('center_X', po_x[0])
feature.SetField('center_Y', po_y[0])
cen_layer.SetFeature(feature)
id.append(FID)
count+=1
area=feature.GetFieldAsDouble('Lake_area')
if area<100:
# count+=1
cen_layer.DeleteFeature(FID)
print("有"+str(count)+"个湖泊中心点不在面内！")
print("总共有"+str(Num_feature)+"个湖库！")
print(id)
# print("面积小于100平方公里的湖库有"+str(count)+"个！")
def GetPoint(shp_path):
txt_path=shp_path[0:-4]+".txt"
driver = ogr.GetDriverByName("ESRI Shapefile")
datasource = driver.Open(shp_path, 1)
cen_layer = datasource.GetLayer()
feature_defn = cen_layer.GetLayerDefn()
# 获取图层定义
feature_count = feature_defn.GetFieldCount()
# 获取属性表中字段个数
count = 0
with open(txt_path,'w') as f:
for feature in cen_layer:
Point_X = feature.GetFieldAsDouble('center_X')
Point_Y = feature.GetFieldAsDouble('center_Y')
f.write(str(Point_X) + "," + str(Point_Y)+'n')
count += 1
f.close()
if __name__ == '__main__':
shp_path=r'C:UsersAdministratorDesktopriversHydroLAKES_polys_v10_shpHydroLAKES_polys_v10.shp'
cal_center(shp_path)
# shp_path=r"C:UsersAdministratorDesktoprivers监测范围-一带一路非洲南部South_Africa_river.shp"
# GetPoint(shp_path)

最后

以上就是无聊自行车最近收集整理的关于学习日记：获取属性表值并修改（判断点是否在面内）的全部内容，更多相关学习日记内容请搜索靠谱客的其他文章。