ICode9

精准搜索请尝试: 精确搜索
首页 > 其他分享> 文章详细

学习笔记3_新建字段,修改字段,geopandas生成点

2021-11-15 18:06:18  阅读:713  来源: 互联网

标签:shp 新建 geopandas self feature 笔记 field path csv


任务:

对面要素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")

标签:shp,新建,geopandas,self,feature,笔记,field,path,csv
来源: https://blog.csdn.net/qq_43210879/article/details/121339140

本站声明: 1. iCode9 技术分享网(下文简称本站)提供的所有内容,仅供技术学习、探讨和分享;
2. 关于本站的所有留言、评论、转载及引用,纯属内容发起人的个人观点,与本站观点和立场无关;
3. 关于本站的所有言论和文字,纯属内容发起人的个人观点,与本站观点和立场无关;
4. 本站文章均是网友提供,不完全保证技术分享内容的完整性、准确性、时效性、风险性和版权归属;如您发现该文章侵犯了您的权益,可联系我们第一时间进行删除;
5. 本站为非盈利性的个人网站,所有内容不会用来进行牟利,也不会利用任何形式的广告来间接获益,纯粹是为了广大技术爱好者提供技术内容和技术思想的分享性交流网站。

专注分享技术,共同学习,共同进步。侵权联系[81616952@qq.com]

Copyright (C)ICode9.com, All Rights Reserved.

ICode9版权所有