[Python知识库]Python读写矢量数据（2）矢量数据写入（属性数据）——Python地理数据处理学习分享

? ? ? ? 这一节主要介绍矢量数据的写入（只有属性数据，无几何），如果有读者没有读取的基础建议先看一下上一篇文章，需要对矢量数据读取有一定的了解才能继续学习本节。在这里我们用到的数据仍为goble文件夹下的数据，数据获取可参考下方链接中的附录获取。Python读写矢量数据（1）针对读取矢量数据——Python地理数据处理学习分享_Cishazhu_ooook的博客-CSDN博客https://blog.csdn.net/remote_giser/article/details/127351424?spm=1001.2014.3001.5502

1.矢量数据的复制

? ? ? ? 当我们想新建矢量文件（shapefile）比复制以前的矢量数据时，可以采用以下方法。

# 导入ogr模块
from osgeo import ogr
# 设置路径为global文件夹，该目录下存在很多shapefile文件，并打开文件夹
path="D:\\Python_study\\osgeopy-data\\global"
doc=ogr.Open(path,1)
# 打开图层
in_layer=doc.GetLayer("ne_50m_populated_places")
# 查询图层，存在删除图层
if doc.GetLayer('capital_cities'):
	doc.DeleteLayer('capital_cities')
# 创建一个新图层，空间参考信息与in_layer一致，并继承其属性字段
our_layer=doc.CreateLayer('capital_cities',in_layer.GetSpatialRef(),ogr.wkbPoint)
our_layer.CreateFields(in_layer.schema)
# 创建一个用于存储属性和几何的虚拟存储要素，插入图层中
our_defn=our_layer.GetLayerDefn()
our_feat=ogr.Feature(our_defn)
# 给要素赋予属性
for in_feat in in_layer:
	if in_feat.GetField('FEATURECLA')=='Admin-0 capital':
		geom=in_feat.geometry()
		our_feat.SetGeometry(geom)
		for i in range(in_feat.GetFieldCount()):
			value=in_feat.GetField(i)
			our_feat.SetField(i,value)
		our_layer.CreateFeature(our_feat)
del doc

? ? ? ? 这个程序打开方式就与众不同，首先我们打开的不是某个shapefile文件，而是将路径定义在某个文件夹下（这个文件夹包含大量shapefile文件），并用编辑的方式打开（Open(x,1/0)设为1为编辑模式，0为只读）。打开这个文件夹后，这个文件夹中的所有shapefile文件都可被看作是一个单独的图层。因此，在我们可以创建一个新图层，就会为我们默认创建一个shapefile文件。接着我们使用GetLayer（）函数获取了"ne_50m_populated_places"文件当作图层in_layer。

? ? ? ? 在创建新图层capital_cities（shapefile文件）之前，我们要确认文件夹中有无这个图层，如果有将其删除，重新创建，这一步的目的是后面我对程序多次运行对图层进行修改（没有这一步，多次运行将会持续报错有同名文件）。用if语句来判断文件夹时候能获取到（GetLayer（））这个图层，如果存在将用DeleteLayer()函数将其删除。

? ? ? ? 随后我们创建了我们的图层our_layer。在这里我们使用的是CreateLayer()函数在文件夹下创建了新的图层our_layer。这里需要注意CreateLayer()函数有四个参数。CreateLayer(name,[srs],[geo_type],[options])如下：

1. name，代表新建的图层名称。在这里我们使用的是'capital_cities'。
2. srs,代表图层使用的空间参考，默认为空（如不输入将表示没有任何空间参考）。这里继承了in_layer的空间参考，采用的是GetSpatialRef()函数（上节介绍过）。
3. geom_type代表一个几何常量，意思是什么类似的几何要素。这里我们使用了点要素ogr.wkbPoint，上节也介绍过了（wkb你可以理解为给电脑看的格式，wkt是给人看的直白的话就是点point\线line\面polygon，在wkb中点对应wkbPoint、线对应wkbLineString、面对应wkbPolygon）想要更多了解可看上节的对应表。
4. options，一个图层创建的的选项列表，只适用于特定的矢量数据类型。这里我未作设置，只有前三个参数（个人没搞懂这个参数，大概不进行输入的以options=value形式存储字符串）。

????????创建图层完成后，我们接下来要给图层创建属性，利用CreateFields（）的方法继承in_layer的属性（schema属性会获得图层FieldDefn对象列表）。到这里我们创建图层属性的任务基本完成。如果这时候运行程序会发现会获得继承文件的属性，可在GIS软件中查看。

? ? ? ? 从代码的14行开始，我们将进行继承每一个要素对应的属性（几何和文本）添加到文件中。首先我们要创建一个用于存储属性和几何的虚拟要素（变量our_defn），将图层对象属性信息传递给这个虚拟要素，这里采用GetLayerDefn()的方法来传递，让out_defn接受in_layer的图层对象属性。?

? ? ? ? 随后我们需要使用ogr.Feature把这个虚拟要素包含的属性真正给一个变量，在这里我们将虚拟要素our_defn的特征给了our_feat（这里称变量吧）。

? ? ? ? 接下来我们将进行要素对应的属性赋值，这也是最繁琐的一步。在被继承的图层in_layer中遍历每个要素（for循环中的feat）。如果某个要素的'FEATURECLA'等于'Admin-0 capital'，建立一个继承要素的遍量geom，使用geometry()方法继承in_feat要素（循环中指每个要素）的几何。再使用SetGeometry（）的方法让我们变量our_feat继承原来图层要素的几何。

? ? ? ? 几何继承完后我们就可以继承属性值，当然两者继承没有先后（只是我觉着有了几个位置再向内部插值比较直观）。就一次循环过程中而言（最外层的for），使用for循环遍历要素in_feat的属性数量，让变量value获取in_field的属性(有for循环，采用的函数为GetField())，最后用变量our_feat使用SetField（）逐一获取单一value，直到读完整个要素的value。这样我们的特征都被存入再our_feat中了，接下来使用CreateFeatyre()让我们图层全部继承其几何和属性。

2.创建新的数据源

2.1创建正确的驱动

????????在创建新的数据源之前我们要了解到，新的数据源创建需要正确驱动。例如，GeoJSON驱动程序并不能创建一个shapefile文件，即使创建后扩展名为.shp，其核心还是一个GeoJSON文件。

? ? ? ? 下面有两者方法来创建所需的驱动。

????????1.继承从打开的数据中进行获取驱动（该方法创建的数据源将会和打开的数据源一致）。

# 导入模块，打开文件
from osgeo import ogr
doc=ogr.Open('D:\\Python_study\\osgeopy-data\\global\\ne_50m_populated_places.shp')
# 获取驱动
driver=doc.GetDriver()

? ? ? ? 在这里使用.GetDriver（）的方式继承了文件doc的驱动。

? ? ? ? 2.使用OGR的GetDriverByName（）函数，传递给他驱动名字。

# 导入模块
from osgeo import ogr
# 设置驱动
new_driver=ogr.GetDriverByName("ESRI Shapefile")

? ? ? ? 这里我们使用是驱动名为ESRI Shapefile,其他驱动的访问可以参考下代码（该代码与创建驱动无关，但我们可能需要特定驱动才能正确创建数据源）：

# 导入模块
from osgeo import ogr
# 查看驱动数量
num=ogr.GetDriverCount()
print(num)
# 查看驱动
for i in range(int(num)):
	print(ogr.GetDriver(i).name)

????????从图可以看到这里可以获取我们常用的格式 。如果不清楚对应文件可取OGR网站查询每种格式对应的文件。

2.2创建新的数据源

? ? ? ? 只有在正确的驱动下我们才可以进行创建正确的数据源。创建数据源代码如下：

# 导入模块
from osgeo import ogr
# 设置驱动
new_driver=ogr.GetDriverByName("ESRI Shapefile")
nwe_doc=new_driver.CreateDataSource('D:\\Python_study\\osgeopy-data\\global\\new_shp')

?

? ? ? ? ?从图片可以看到，我们在文件夹目录下创建了一个new_shp文件（数据源是个文件夹）。注意，创建新的数据源时，不能覆盖现有的数据源，那么我们就需要在创建新数据源之前删除旧的数据源，目前采用的方法是Python的os.path.exists()函数处理。

# 导入模块
from osgeo import ogr
import os
# 设置驱动
new_driver=ogr.GetDriverByName("ESRI Shapefile")
# 检查是否已经存在
if os.path.exists("D:\\Python_study\\osgeopy-data\\global\\new_shp"):
	new_driver.DeleteDataSource("D:\\Python_study\\osgeopy-data\\global\\new_shp")
# 创建新数据源
nwe_doc=new_driver.CreateDataSource("D:\\Python_study\\osgeopy-data\\global\\new_shp")

2.3新建属性字段

????????OGR定义字段与数据库类型时，首先需查看有哪些字段类型。这里统计了一些常见的字段数据类型在OGR中的对应。

? ? ? ? 要将一个字段添加到图层中，首先需要一个包含字段名称、数据类型、字段宽度、精度等重要信息的对象（使用ogr.FieldDefn()函数）。然后基于这个对象将其添加到图层。

# 导入模块
from osgeo import ogr
import os
# 设置驱动
new_driver=ogr.GetDriverByName("ESRI Shapefile")
# 检查是否已经存在
if os.path.exists("D:\\Python_study\\osgeopy-data\\global\\new_shp"):
	new_driver.DeleteDataSource("D:\\Python_study\\osgeopy-data\\global\\new_shp")
# 创建新数据源
new_doc=new_driver.CreateDataSource("D:\\Python_study\\osgeopy-data\\global\\new_shp")
# 创建新图层
layer=new_doc.CreateLayer("new_layer")
# 定义字段
add_1=ogr.FieldDefn("add_the_1",ogr.OFTInteger)
add_1.SetWidth(8)
add_1.SetPrecision(3)
# 将定义字段添加到图层
layer.CreateField(add_1)

? ? ? ? 别忘了，我们在这里创建完数据源后记得添加新的图层（即shape file文件，在这里我忘了添加坐标系和类型。），然后用变量add_1指向新建的字段并输入字段名称和类型（在对应表中寻找）。随后我们设置了新字段的宽度和精度（分别用到了SetWidth()、SetPrecision()）函数。最后使用CreatField（）添加到了图层中，并在GIS软件中可以到如下结果。

? ? ? ? ?当然有时我们可能想要创建仅仅是字段名不同的两个字段，这里会提供SetName()函数来改变字段名而其他不发生改变，但需要重新将第二个字段添加到属性表中（实际修改了变量然后增进行了第二次增加）。

# 导入模块
from osgeo import ogr
import os
# 设置驱动
new_driver=ogr.GetDriverByName("ESRI Shapefile")
# 检查是否已经存在
if os.path.exists("D:\\Python_study\\osgeopy-data\\global\\new_shp"):
	new_driver.DeleteDataSource("D:\\Python_study\\osgeopy-data\\global\\new_shp")
# 创建新数据源
new_doc=new_driver.CreateDataSource("D:\\Python_study\\osgeopy-data\\global\\new_shp")
# 创建新图层
layer=new_doc.CreateLayer("new_layer")
# 定义字段
add_1=ogr.FieldDefn("add_the_1",ogr.OFTInteger)
add_1.SetWidth(8)
add_1.SetPrecision(3)
# 将定义字段添加到图层
layer.CreateField(add_1)
# 仅名称不一样的字段
add_1.SetName("add_the_2")
layer.CreateField(add_1)

? ? ? ? 查看结果如下：

2.4修改字段信息?

? ? ? ? 对于实例中的add_the_1名称不满意，在这里我们可以对字段就行修改（包括精度、类型等）。采用AlterFieldDefn(iField,field_def,nFlags)进行修改。

1. iField，索引的字段名，本例去索引的为add_the_1
2. field_def，新属性字段定义
3. nFlags，改变的属性，下表代表属性对应的代码

?????????

# 导入模块
from osgeo import ogr
import os
# 设置驱动
new_driver=ogr.GetDriverByName("ESRI Shapefile")
# 检查是否已经存在
if os.path.exists("D:\\Python_study\\osgeopy-data\\global\\new_shp"):
	new_driver.DeleteDataSource("D:\\Python_study\\osgeopy-data\\global\\new_shp")
# 创建新数据源
new_doc=new_driver.CreateDataSource("D:\\Python_study\\osgeopy-data\\global\\new_shp")
# 创建新图层
layer=new_doc.CreateLayer("new_layer")
# 定义字段
add_1=ogr.FieldDefn("add_the_1",ogr.OFTInteger)
add_1.SetWidth(8)
add_1.SetPrecision(3)
# 将定义字段添加到图层
layer.CreateField(add_1)
# 修改字段信息
i=layer.GetLayerDefn().GetFieldIndex('add_the_1')
fld_defn=ogr.FieldDefn("change_1",ogr.OFTString)
fld_defn.SetWidth(10)
fld_defn.SetPrecision(6)
layer.AlterFieldDefn(i,fld_defn,ogr.ALTER_ALL_FLAG)

? ? ? ? 其实这串代码也十分好解释，首先让 i 获取到要修改的字段（寻找字段名），重新定义一个字段信息并设置精度和宽度。最后使用上述提过的函数来完成字段的改变，当然我们也可以针对单一进行改变。

?要素添加、更新和删除

? ? ? ? 经过上方的学习，我们已经能修改数据源、图层、属性表的字段了，但我们现在还需要对属性表的要素进行更改。（下方代码用到的数据是50m人口）

# 导入模块
from osgeo import ogr
# 读取文件路径
fn="D:\\Python_study\\osgeopy-data\\global\\ne_50m_populated_places.shp"
# 打开文件
ds=ogr.Open(fn,1)
# 打开图层
lyr=ds.GetLayer(0)
# 添加新字段
add_1=ogr.FieldDefn("ID_add",ogr.OFTInteger)
add_1.SetWidth(8)
add_1.SetPrecision(3)
lyr.CreateField(add_1)
# 要素添加
n=1
for feat in lyr:
	feat.SetField("ID_add",n)
	lyr.SetFeature(feat)
	n+=1
	if n==10:
		break
# 关闭文件
del ds

? ? ? ? ?首先我们要找到我们需要添加的字段名称，然后使用for循环来改变要素。SetField("字段名“，更改内容)是用来确定某一要素对应的字段和要改变的内容。SetFeature()是对这个要素对应字段进行改变。

改变单一要素对应值：

# 导入模块
import sys
from osgeo import ogr
# 读取文件路径
fn="D:\\Python_study\\osgeopy-data\\global\\ne_50m_populated_places.shp"
# 打开文件
ds=ogr.Open(fn,1)
if ds is None:
	sys.exit('Could not open{0}.'.format(fn))
# 打开图层
lyr=ds.GetLayer(0)
# 改变单一要素字段
for feat in lyr:
	if feat.GetField('ID_add')==4:
		feat.SetField("ID_add",99)
		lyr.SetFeature(feat)
# 关闭文件
del ds

? ? ? ? 从上面我们可以知道，要获取某要素时我们往往通过for来进行遍历，然后让其等于某字段的内容。（常采用FID或者名称）这是有一种常用的矢量写入的方法。

? ? ? ? 对于要素的删除，则更简单。我们只需要知道要素某个字段的值，便可删除。?这里使用DeleteFeature()的方法将循环参数对应FID的要素删除。（FID的获取使用GetFID()）

# 导入模块
import sys
from osgeo import ogr
# 读取文件路径
fn="D:\\Python_study\\osgeopy-data\\global\\ne_50m_populated_places.shp"
# 打开文件
ds=ogr.Open(fn,1)
if ds is None:
	sys.exit('Could not open{0}.'.format(fn))
# 打开图层
lyr=ds.GetLayer(0)
# 要素删除
for feat in lyr:
	if feat.GetField('ID_add')==5:
		lyr.DeleteFeature(feat.GetFID())
# 关闭文件
del ds