9_OGR创建Shapefile

OGR创建Shapefile

# -*- coding: utf-8 -*-
# 学习时间： 2022/11/7 17:20

__author__ = 'He XK'

import os
from osgeo import ogr

# 新建shp驱动
driver = ogr.GetDriverByName('ESRI SHAPEFILE')
# 参考的shp文件、图层和参考系
driver_ref = ogr.Open('shp/province_1.shp')
layer_ref = driver_ref.GetLayer(0)
srs_ref = layer_ref.GetSpatialRef()
# 指定shp路径
output_filename = 'shp/demo17.shp'
# 判断是否存在新建文件，否则删除
if os.path.exists(output_filename):
    driver.DeleteDataSource(output_filename)
else:
    print('THIS SHP IS NOT EXIST!')
# 创建数据
new_shp_ds = driver.CreateDataSource(output_filename)
# 创建一个图层
# 当新建图层之后，shp文件才会生成
new_layer = new_shp_ds.CreateLayer('test', srs=srs_ref, geom_type=ogr.wkbPoint)
# 新建字段
new_field_def = ogr.FieldDefn('ID', ogr.OFTString)
new_field_def.SetWidth(4)
new_layer.CreateField(new_field_def)
new_field_def2 = ogr.FieldDefn('Name', ogr.OFTString)
new_field_def2.SetWidth(10)
new_layer.CreateField(new_field_def2)
# 为图层创建feature
feature_defn = new_layer.GetLayerDefn()
feature = ogr.Feature(feature_defn)
# 创建几何
new_point = ogr.Geometry(ogr.wkbPoint)
new_point.SetPoint(0, 500000.822960446, 3000000.31432952)
feature.SetGeometry(new_point)
# 设定要素字段的值
feature.SetField('ID', '1')
feature.SetField('Name', 'TEST1')
# 将要素写入图层
new_layer.CreateFeature(feature)
# 释放内存
new_point.Destroy()
feature.Destroy()
new_shp_ds.Destroy()
driver_ref.Destroy()

使用ogr创建一个shp文件，并且生成一个点坐标，添加两个字段并填入内容。

整个逻辑过程如下：

1. 首先，生成SHP的Driver；

2. 判断当前路径是否已经有同名的SHP文件存在，若存在则通过os模块删除；

3. 通过Driver创建数据集，即创建矢量文件；

4. 为矢量文件新建图层，同时在参数中设定图层名称、空间参考系、几何类型；

5. 设定字段，并在矢量文件的图层中添加字段；

6. 为图层创建feature；

7. 创建几何图形，并设定几何图形的坐标，将其添加到feature中；

8. 为feature设定字段属性；

9. 将feature写入layer中；

10. 释放内存。

设定空间参考系的方式不止一种，可以通过引入等方式新建参考系并设定；

创建几何图形的方式还能通过wkb的方式创建，适合创建不定长度的几何图形；

使用时查阅相关资料，可以封装一个输入任意坐标系、字段、不定长几何坐标的函数。

注意：

new_layer = new_shp_ds.CreateLayer('test', srs=srs_ref, geom_type=ogr.wkbUnknown)

在上述代码创建新图层时，geom_type是一个约束选项，值是用于约束图层中的几何图形必须是该值类型的图形，如果没有约束那么使用wkbUnknown。

创建点要素

创建点要素使用函数SetPoint或者AddPoint

AddPoint(Geometry self, double x, double y, double z=0)

SetPoint(Geometry self, int point, double x, double y, double z=0)

可以添加多个点，但是只会收录最后一次添加的点。

如果给同一个图层中创建多个点要素，那么就需要对每一个点要素，新建一个feature，同时执行一遍设定字段和写入图层。

创建线要素

创建线要素的主体框架和点要素一致，只是layer的geom_type=ogr.wkbLineString，同时ogr.Geometry(ogr.wkbLineString)，添加线要素的折点时，需要使用多次的AddPoint

new_line.AddPoint(500030.822960446, 3000030.31432952)
new_line.AddPoint(506030.822960446, 3005030.31432952)
new_line.AddPoint(503030.822960446, 3001030.31432952)

也可以通过设定一段wkt

wkt = 'LINESTRING(116.4 39.9, 116.42 39.91, 116.39 39.92, 116.38 39.91)' ## 创建线
line = ogr.CreateGeometryFromWkt(wkt) ## 生成线
feature.SetGeometry(line)

也有其他的一些方式，可以在文档中查询

[外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传(img-7n4UQ0MK-1669793750867)(file://I:\MarkDownFiles\GDAL\GDAL_Pic\2022-11-08-15-22-19-image.png?msec=1669793737872)]

其他一些有用函数：

line.GetPointCount()获取线段中的点数量；

line.GetX(0)， line.GetY(0)获取指定点的XY坐标。

wkbLinearRing线环（问题已解决）

同属wkbLineString，是线要素下的分支。

在创建图层时选择普通的线：

new_layer = new_shp_ds.CreateLayer('test', srs=srs_ref, geom_type=ogr.wkbLineString)

创建几何时选择：

ring = ogr.Geometry(ogr.wkbLinearRing)

wkbLinearRing和wkbLineString的最大区别在于，前者可以使用ring.CloseRings()使线环要素闭环，后者不行。

创建面要素

new_layer = new_shp_ds.CreateLayer('test', srs=srs_ref, geom_type=ogr.wkbPolygon)


wkt = 'POLYGON((506565.822960446 3006854.31432952, 516565.822960446 3006854.314329522 , 516565.822960446 3106854.31432952, 506565.822960446 3006854.31432952)) '
polygon = ogr.CreateGeometryFromWkt(wkt)
feature.SetGeometry(polygon)

其中wkt必须是POLYGON((...))两层括号。

如何实现中心挖空的多边形？据说是使用LinearRing，但是LinearRing在这里报错。

不使用LinearRing，可以写两个多边形区域。

wkt = 'POLYGON((1 1,5 1,5 5,1 5,1 1),(2 2,2 3,3 3,3 2,2 2))'

Polygon的创建还可以通过其他方法，比如写两个Ring，然后通过AddGeometry方法形成多边形。

对于POLYGON有几个函数：

polygon.GetGeometryCount()计算当前polygon的图形个数，若是一个没有挖空的多边形，结果是1；挖空一个洞，结果是2，以此类推；

多边形的坐标获取：

print(polygon)  # 直接输出对象

POLYGON ((1 1,5 1,5 5,1 5,1 1),(2 2,2 3,3 3,3 2,2 2))

outring = polygon.GetGeometryRef(0)
inring = polygon.GetGeometryRef(1)
print(outring, inring)

LINEARRING (1 1,5 1,5 5,1 5,1 1) LINEARRING (2 2,2 3,3 3,3 2,2 2) 

创建复合几何形状MultiGeometry

Geometry — Python GDAL/OGR Cookbook 1.0 documentation (pcjericks.github.io)

创建复合的几何形状，比如wkbMultiPoint多点。

在创建图层时选择多点：

new_layer = new_shp_ds.CreateLayer('test', srs=srs_ref, geom_type=ogr.wkbMultiPoint)

创建几何时有些不同：

# 创建几何
multipoint = ogr.Geometry(ogr.wkbMultiPoint)

point1 = ogr.Geometry(ogr.wkbPoint)
point1.AddPoint(30, 100)
multipoint.AddGeometry(point1)

point2 = ogr.Geometry(ogr.wkbPoint)
point2.AddPoint(25, 60)
multipoint.AddGeometry(point2)

print(multipoint.ExportToWkt())

也可以共用一个Point。

创建多线wkbMultiLineString

multiLine = ogr.Geometry(ogr.wkbMultiLineString)

wkt = 'LINESTRING(15 20, 25 60, 41 60, 188 36, 123 656)'
line1 = ogr.CreateGeometryFromWkt(wkt)
multiLine.AddGeometry(line1)

wkt = 'LINESTRING(60 58, 14 58, 69 78, 14 25, 146 35, 54 96)'
line2 = ogr.CreateGeometryFromWkt(wkt)
multiLine.AddGeometry(line2)

feature.SetGeometry(multiLine)

创建一个GeometryCollection

# 创建几何
# Create a geometry collection
geomcol = ogr.Geometry(ogr.wkbGeometryCollection)
# Add a point
point = ogr.Geometry(ogr.wkbPoint)
point.AddPoint(-122.23, 47.09)
geomcol.AddGeometry(point)

# Add a line
line = ogr.Geometry(ogr.wkbLineString)
line.AddPoint(-122.60, 47.14)
line.AddPoint(-122.48, 47.23)
geomcol.AddGeometry(line)
print(geomcol.ExportToWkt())

创建其他内容的基本结构就是这样，但是问题在于对于某些图形，例如这里的geomcol = ogr.Geometry(ogr.wkbGeometryCollection)，并不能够将它添加进feature.SetGeometry(geomcol)中。

如果layer的图层约束设置为wkbUnkonwn，那么会产生ERROR 1: Attempt to write non-linestring (GEOMETRYCOLLECTION) geometry to ARC type shapefile.错误。

使用OGR拷贝方法创建新的Shapefile

使用OGR拷贝创建新的矢量数据集的方法其实就是将数据读取的顺序再重复一遍。

总体来说这个过程就是构建数据源->构建层->构建要素->构建几何形状->关闭数据源。

在OGR的矢量数据模型中，矢量数据分为datasource，layer，feature三个层次，在这三个层次上，OGR均有不同的复制数据的方法。

Datasource

直接复制源文件，重新写入新的文件中。

DataSource的复制，在于打开新文件后，直接用新文件的驱动复制原来的文件。

output_ds = driver.CopyDataSource(src_file, output_file)

---

其中新增了判断文件是否存在并且删除已有文件的过程。

矢量的Datasource具备Release释放的方法。

似乎Release和Destroy的效果一致。

这样复制会出现编码问题

Warning 1: One or several characters couldn't be converted correctly from UTF-8 to ISO-8859-1. This warning will not be emitted anymore.

from osgeo import ogr, gdal
import os

driver = ogr.GetDriverByName('ESRI SHAPEFILE')
gdal.SetConfigOption("GDAL_FILENAME_IS_UTF8", "YES")
gdal.SetConfigOption("SHAPE_ENCODING", "UTF8")
src_file_path = 'shp/province_1.shp'
output_file = 'shp/demo21_dscopy.shp'
if os.path.exists(output_file):
    driver.DeleteDataSource(output_file)
else:
    print('THIS SHP IS NOT EXIST!')

src_file = driver.Open(src_file_path)
output_ds = driver.CopyDataSource(src_file, output_file)
output_ds.Release()

尝试解决：

查询SetConfigOption，找到了ESRI Shapefile / DBF下有关于Encoding的内容，但是不太有用的样子。

• ENCODING=encoding_name: to override the encoding interpretation of the shapefile with any encoding supported by CPLRecode or to “” to avoid any recoding`

已解决：

以上代码在Python中循环查询Feature的属性，能够正确输出内容，但是，在ArcGIS中查看属性却依然乱码，目测可能是ArcGIS的解码出现问题。

在QGIS中打开SHP文件，调整编码格式后可以正常显示。

Layer

Layer的复制在于，打开目标文件，打开新文件，获取目标文件的图层，在新文件的驱动下直接对原始文件的图层进行复制。

output_layer = output_ds.CopyLayer(layer, 'test')

from osgeo import ogr, gdal
import os

driver = ogr.GetDriverByName('ESRI SHAPEFILE')
gdal.SetConfigOption("GDAL_FILENAME_IS_UTF8", "YES")
gdal.SetConfigOption("SHAPE_ENCODING", "UTF8")
src_file_path = 'shp/province_1.shp'
output_file = 'shp/demo21_dscopy.shp'
if os.path.exists(output_file):
    driver.DeleteDataSource(output_file)
else:
    print('THIS SHP IS NOT EXIST!')

src_file = driver.Open(src_file_path)
layer = src_file.GetLayer(0)
output_ds = driver.CreateDataSource(output_file)
output_layer = output_ds.CopyLayer(layer, 'test')

Feature

在feature层面的复制并没有什么方法，而是获取到原有文件的feature后在新文件中创建feature并且写入。

写入的同时不会保留属性字段，同样需要获取后写入。

from osgeo import ogr, gdal
import os

driver = ogr.GetDriverByName('ESRI SHAPEFILE')
gdal.SetConfigOption("GDAL_FILENAME_IS_UTF8", "YES")
gdal.SetConfigOption("SHAPE_ENCODING", "UTF8")
src_file_path = 'shp/province_1.shp'
output_file = 'shp/demo22_featurecopy.shp'
if os.path.exists(output_file):
    driver.DeleteDataSource(output_file)
else:
    print('THIS SHP IS NOT EXIST!')

src_file = driver.Open(src_file_path)
layer = src_file.GetLayer(0)
output_ds = driver.CreateDataSource(output_file)
output_layer = output_ds.CreateLayer('test', srs=layer.GetSpatialRef(), geom_type=ogr.wkbPolygon)
# 创建字段
layer_info = layer.GetLayerDefn()
output_field_def = ogr.FieldDefn('Name', ogr.OFTString)
output_field_def.SetWidth(10)
output_layer.CreateField(output_field_def)

layer.ResetReading()
i = 0
while i < layer.GetFeatureCount():
    feature = layer.GetFeature(i)
    i += 1
    print(i)
    output_feature = output_layer.CreateFeature(feature)
    feature.SetField('Name', feature.GetField('Name'))

OGR创建属性字段

在上述代码中，我们已经使用过了如何创建字段。

layer_info = layer.GetLayerDefn()
output_field_def = ogr.FieldDefn('Name', ogr.OFTString)
output_field_def.SetWidth(10)
output_layer.CreateField(output_field_def)

有一个方法：

layer_info.AddFieldDefn(output_field_def)

当OGRFeatureDefn中无OGRFeature对象时，此方法才可以调用。此方法的优点是保证在不影响调用者响应的同时，输入的OGRFieldDefn也会被拷贝备份。

可以看出，输出的图层属性表中并无内容。可能是这个函数只在内存中虚拟地创建了一个字段，与磁盘并无关系，或者说这个函数只是用来创建一个属性字段的框架，仅仅用来定义，在建立新的数据时可能起到参考与引用的作用。在软件处理过程中，由于软件的问题可能会导致图形对象与属性记录并非一一对应，严重时矢量数据可能会损坏。

该部分暂时不知道有什么用

添加投影Projection

投影可以从DataSource、Layer、Geometry等具备投影信息的对象中获取。

一般投影可以从参考系中获取然后赋予给新的shp文件。

也可以直接创建一个投影。

创建投影

创建投影使用osr模块。

可以根据EPSG的数字代号进行投影创建。

也有导出投影信息的函数。

from osgeo import osr
spatial_ref = osr.SpatialReference()
spatial_ref.ImportFromEPSG(3857)

重投影一个Geometry

__author__ = 'He XK'
from osgeo import osr, ogr
spatial_ref = osr.SpatialReference()
spatial_ref.ImportFromEPSG(3857)

sr_spatial = osr.SpatialReference()
sr_spatial.ImportFromEPSG(2927)

transform = osr.CoordinateTransformation(sr_spatial, spatial_ref)

point = ogr.CreateGeometryFromWkt('POINT (1120351.57 741921.42)')
print(point)  # POINT (1120351.57 741921.42)
point.Transform(transform)
print(point)  # POINT (-13647563.6238682 5999195.26697041)

Geometry对象具备Transform函数。

首先需要定义一个对象transform = osr.CoordinateTransformation(source, target)来存储投影坐标系。之后进行转换，可以实现重投影。

投影的更改只能修改Feature，所以，对一个SHP文件进行重投影时，需要遍历Feature

投影写入文件

对已经获取的投影信息对象进行spatial_ref.MorphToESRI()格式转换，转化成WKT格式。

之后写入文件导出即可。

# 投影写入文件
spatial_ref.MorphToESRI()
file = open('shp/demo23_output_ref.prj', 'w')
file.write(spatial_ref.ExportToWkt())
file.close()