我们以Shapefile为例,使用以下函数创建新的数据源(DataSource):
driver.CreateDataSource(<filename>)
需注意的是,应先判断这个要被创建的文件(filename)不能存在,否则会出错。 通常需要先判断一下文件是否已存在,若存在的话,需使用其他文件名,或将已存在的矢量数据直接删除。
删除矢量数据
大多数的GIS数据格式,如shapefile,mapinfo tab等,都不是单一文件。
像 Shapefile,除了最基本的shp文件外,还有保存属性的dbf文件。
因此,在删除GIS数据时,可使用 os
模块删除,不过需要查找相关的文件进行全部删除,较之OGR的删除数据函数 DeleteDataSource()
要烦琐得多。
import os
from osgeo import ogr
driver=ogr.GetDriverByName('ESRI Shapefile')
out_shp = '/tmp/xx_world_borders.shp'
if os.path.exists(out_shp):
driver.DeleteDataSource(out_shp)
dir(out_shp)[:6] + ['... ...'] + dir(out_shp)[-3:]
['__add__', '__class__', '__contains__', '__delattr__', '__dir__', '__doc__', '... ...', 'translate', 'upper', 'zfill']
然后,使用下面的语句创建数据:
ds =driver.CreateDataSource(out_shp)
layer=ds.CreateLayer('test',geom_type=ogr.wkbPoint)
/opt/conda/lib/python3.12/site-packages/osgeo/gdal.py:311: FutureWarning: Neither gdal.UseExceptions() nor gdal.DontUseExceptions() has been explicitly called. In GDAL 4.0, exceptions will be enabled by default. warnings.warn(
若想添加一个新字段,只能在layer里加,而且还得保证layer里没有数据。
若添加的字段是字符串,还需设定宽度。
fieldDefn = ogr.FieldDefn( 'id',ogr.OFTString)
fieldDefn.SetWidth(4)
layer.CreateField(fieldDefn)
0
若添加一个新的feature,首先得完成上一步,把字段field都添加完整。
然后从layer中读取相应的feature类型,并创建feature:
featureDefn = layer.GetLayerDefn()
feature = ogr.Feature(featureDefn)
设定几何形状:
point = ogr.Geometry(ogr.wkbPoint)
point.SetPoint(0,123,123)
feature.SetGeometry(point)
0
设定某字段数值:
feature.SetField('id', 23)
将feature写入layer:
layer.CreateFeature(feature)
0
最后,从内存中清除 ds
,将数据写入磁盘中。
ds.Destroy()
OGR中定义的几何类型
OGR关键字 | 名称 | 编码值 |
---|---|---|
wkbUnknown | 未知类型 | 0 |
wkbPoint | 点 | 1 |
wkbLineString | 线 | 2 |
wkbPolygon | 多边形 | 3 |
wkbMultiPoint | 多点 | 4 |
wkbMultiLineString | 多线 | 5 |
wkbMultiPolygon | 多部分多边形 | 6 |
wkbGeometryCollection | 几何形状集合 | 7 |
wkbNone | 无 | 100 |
wkbLinearRing | 线环 | 101 |
wkbPoint25D | 2.5维点 | 0x80000001 |
wkbLineString25D | 2.5维线 | 0x80000002 |
wkbPolygon25D | 2.5维多边形 | 0x80000003 |
wkbMultiPoint25D | 2.5维多点 | 0x80000004 |
wkbMultiLineString25D | 2.5维多线 | 0x80000005 |
wkbMultiPolygon25D | 2.5维多部分多边形 | 0x80000006 |
wkbGeometryCollection25D | 2.5维几何形状集合 | 0x80000007 |
from osgeo import ogr
point = ogr.Geometry(ogr.wkbPoint)
point.AddPoint(10,20)
print(point)
POINT Z (10 20 0)
from osgeo import ogr
line = ogr.Geometry(ogr.wkbLineString)
line.AddPoint(10,10)
line.AddPoint(20,20)
print (line)
LINESTRING Z (10 10 0,20 20 0)
使用线对象的SetPoint修改点坐标,例如:
print (line)
LINESTRING Z (10 10 0,20 20 0)
另外,还有一些其他有用的函数,如统计所有点数量:
print (line.GetPointCount())
2
又如读取0号点的x坐标和y坐标:
print (line.GetX(0))
10.0
print (line.GetY(0))
10.0
ring = ogr.Geometry(ogr.wkbLinearRing)
ring.AddPoint(0,0)
ring.AddPoint(100,0)
ring.AddPoint(100,100)
ring.AddPoint(0,100)
在结尾处,用命令CloseRings
闭合ring,或者设定最后一个点坐标与第一个点相同,也可以闭合ring。
ring.CloseRings()
print(ring)
LINEARRING Z (0 0 0,100 0 0,100 100 0,0 100 0,0 0 0)
outring = ogr.Geometry(ogr.wkbLinearRing)
outring.AddPoint(0,0)
outring.AddPoint(100,0)
outring.AddPoint(100,100)
outring.AddPoint(0,100)
outring.AddPoint(0,0)
这里将最后点坐标设置与起始点相同,来闭合ring:
inring = ogr.Geometry(ogr.wkbLinearRing)
inring.AddPoint(25,25)
inring.AddPoint(75,25)
inring.AddPoint(75,75)
inring.AddPoint(25,75)
inring.CloseRings()
polygon = ogr.Geometry(ogr.wkbPolygon)
polygon.AddGeometry(outring)
0
polygon.AddGeometry(inring)
0
最后三行命令比较重要:先建立一个polygon对象,然后再添加外层ring和内层ring。
若查看polygon有几个ring,可运行下列代码:
print (polygon.GetGeometryCount())
2
对于polygon对象,无法直接输出坐标,而是先获取几何对象,也就是ring。从polygon中读取ring、index的顺序与创建polygon的顺序相同。
outring2 = polygon.GetGeometryRef(0)
inring2 = polygon.GetGeometryRef(1)
print(outring2)
LINEARRING Z (0 0 0,100 0 0,100 100 0,0 100 0,0 0 0)
print(inring2)
LINEARRING Z (25 25 0,75 25 0,75 75 0,25 75 0,25 25 0)
multipoint = ogr.Geometry(ogr.wkbMultiPoint)
point = ogr.Geometry(ogr.wkbPoint)
point.AddPoint(10,10)
multipoint.AddGeometry(point)
0
point.AddPoint(20,20)
multipoint.AddGeometry(point)
0
读取MultiGeometry的Geometry,其方法与polygon读取ring是一样的,因此Polygon是一种内置的MultiGeometry。
注意:不要删除一个已存在Feature geometry,可能会导致python崩溃,可以删除脚本运行期间创建的Geometry,可能是手动创建,或者调用其他函数自动创建的。就算这个Geometry已经创建了其它的Feature,你还是可以删除它:
例如:Polygon.Destroy()
使用OGR创建数据集的几何形状
对于创建Geometry来讲,wkt是比较直观的。wkt是最简单的字符串格式,而且Python提供了大量的函数来处理字符串。使用wkt创建矢量数据集与拷贝创建数据集的基本步骤是一致的。首先是创建Driver,按照顺序依次创建:dataSource,Layer,Feature。需要分别创建点、线、与多边形。使用extent
变量来存储多边形的四个角点,点与线坐标根据extent
变量来生成,如:
extent = [400, 1100, 300, 600]
创建点数据集
一般情况不会生成单独的点,需要使用For循环语句来对数组参数进行遍历。
创建线数据集
注意坐标的格式:不同点坐标之间用“,”隔断,坐标的不同维度用空格隔断。
创建多边形数据集
首先使用wkt来定义一个矩形。矩形的范围由extent
变量定义,矩形四个角点坐标根据extend生成,存储到wkt变量中。通过CreateGeometryFromWkt()
创建一个矩形形状,并放入Feature中,最后调用Destroy(),将数据写入磁盘,这样就构成了一个矩形的数据集。
from osgeo import ogr
import os,math
inshp = '/data/gdata/GSHHS_c.shp'
ds = ogr.Open(inshp)
driver = ogr.GetDriverByName("ESRI Shapefile")
outputfile = '/tmp/xx_world_borders_copy.shp'
if os.access( outputfile, os.F_OK ):
driver.DeleteDataSource(outputfile)
pt_cp = driver.CopyDataSource(ds, outputfile)
pt_cp.Release()
dir(pt_cp)[:6] + ['... ...'] + dir(pt_cp)[-3:]
Warning 1: Value 50654050.694499999 of field area of feature 0 not successfully written. Possibly due to too larger number with respect to field width Warning 1: Value 50654050.694499999 of field area of feature 1 not successfully written. Possibly due to too larger number with respect to field width Warning 1: Value 29220969.727000002 of field area of feature 2 not successfully written. Possibly due to too larger number with respect to field width Warning 1: Value 20154740.09 of field area of feature 3 not successfully written. Possibly due to too larger number with respect to field width Warning 1: Value 17534164.066799998 of field area of feature 4 not successfully written. Possibly due to too larger number with respect to field width Warning 1: Value 13919140.565400001 of field area of feature 1736 not successfully written. Possibly due to too larger number with respect to field width
['AbortSQL', 'AddBand', 'AddFieldDomain', 'AddRelationship', 'AdviseRead', 'BeginAsyncReader', '... ...', '_invalidate_children', 'this', 'thisown']
from osgeo import ogr
import os,math
inshp = '/data/gdata/GSHHS_c.shp'
ds = ogr.Open(inshp)
driver = ogr.GetDriverByName("ESRI Shapefile")
outputfile ='/tmp/xx_world_borders_copy2.shp'
if os.access( outputfile, os.F_OK ):
driver.DeleteDataSource( outputfile )
layer = ds.GetLayer()
newds = driver.CreateDataSource(outputfile)
pt_layer = newds.CopyLayer(layer, 'abcd')
newds.Destroy()
dir(pt_layer)[:6] + ['... ...'] + dir(pt_layer)[-3:]
Warning 1: Value 50654050.694499999 of field area of feature 0 not successfully written. Possibly due to too larger number with respect to field width Warning 1: Value 50654050.694499999 of field area of feature 1 not successfully written. Possibly due to too larger number with respect to field width Warning 1: Value 29220969.727000002 of field area of feature 2 not successfully written. Possibly due to too larger number with respect to field width Warning 1: Value 20154740.09 of field area of feature 3 not successfully written. Possibly due to too larger number with respect to field width Warning 1: Value 17534164.066799998 of field area of feature 4 not successfully written. Possibly due to too larger number with respect to field width Warning 1: Value 13919140.565400001 of field area of feature 1736 not successfully written. Possibly due to too larger number with respect to field width
['AlterFieldDefn', 'AlterGeomFieldDefn', 'Clip', 'CommitTransaction', 'CreateFeature', 'CreateField', '... ...', 'schema', 'this', 'thisown']
from osgeo import ogr
import os,math
inshp = '/data/gdata/GSHHS_c.shp'
ds = ogr.Open(inshp)
driver = ogr.GetDriverByName("ESRI Shapefile")
outputfile ='/tmp/xx_world_borders_copy3.shp'
if os.access( outputfile, os.F_OK ):
driver.DeleteDataSource( outputfile )
newds = driver.CreateDataSource(outputfile)
layernew = newds.CreateLayer('worldcopy',None,ogr.wkbLineString)
layer = ds.GetLayer()
extent = layer.GetExtent()
extent
(-180.0, 180.00000000000023, -89.99999999999994, 83.53036100000006)
from osgeo import ogr
import os,math
inshp = '/data/gdata/GSHHS_c.shp'
ds = ogr.Open(inshp)
driver = ogr.GetDriverByName("ESRI Shapefile")
outputfile ='/tmp/xx_world_borders_copy3.shp'
if os.access( outputfile, os.F_OK ):
driver.DeleteDataSource( outputfile )
newds = driver.CreateDataSource(outputfile)
layernew = newds.CreateLayer('worldcopy',None,ogr.wkbLineString)
layer = ds.GetLayer()
feature = layer.GetNextFeature()
while feature is not None:
layernew.CreateFeature(feature)
feature = layer.GetNextFeature()
newds.Destroy()
dir(feature)[:6] + ['... ...'] + dir(feature)[-3:]
['__bool__', '__class__', '__delattr__', '__dir__', '__doc__', '__eq__', '... ...', '__sizeof__', '__str__', '__subclasshook__']
特别注意:Destroy()
不能省略,Destroy()
除了销毁数据还有数据flush到磁盘的作用。如果没有此命令,那么刚才创建的一系列要素就不会被写入磁盘的文件中。
OGR属性字段的定义与使用
GIS数据除了图形要素之外,更重要的就是属性数据。其实GIS数据的属性可以理解成关系型数据库。只不过在这个关系数据库中,有一列特定的索引值,且数据库中的记录与图形对象也是一一对应的。
在OGR中定义属性字段与赋值
OGR定义字段与数据库类型时,首先需查看有哪些字段类型(下表)。
OGR中定义的字段类型
枚举值 | 类型 | 位数 | 编码值 |
---|---|---|---|
OFTInteger | 整型 | 32位整型 | 0 |
OFTIntegerList | 整型列表 | 32位整型 | 1 |
OFTReal | 实数型 | 双精度浮点型 | 2 |
OFTRealList | 实数型列表 | 双精度浮点型 | 3 |
OFTString | 字符串 | ASCII字符 | 4 |
OFTStringList | 字符串列表 | ASCII字符串数组 | 5 |
OFTWideString | 宽字符串 | 不再使用 | 6 |
OFTWideStringList | 宽字符串列表 | 不再使用 | 7 |
OFTBinary | 二进制型 | 原始的二进制数据 | 8 |
OFTDate | 日期型 | 9 | |
OFTTime | 时间型 | 10 | |
OFTDateTime | 日期时间型 | 11 |
我们还可以给这个矩形添加属性数据,首先定义表头,也就是定义字段,需强调的是定义一个字段需要定义其字段描述(包括字段的名称与类型),然后将这个字段描述创建到layer中。有了表头,就可以使用SetField
方法在输入Feature时添加属性表内容,这时打开QGIS ShapeFile时就可以看到结果了。若使用中文的字段值,请注意字符串字段的宽度,若超出宽度,显示内容就会出现问题。
定义字段的位置
当OGRFeatureDefn中无OGRFeature对象时,此方法才可以调用。此方法的优点是保证在不影响调用者响应的同时,输入的OGRFieldDefn也会被拷贝备份。
import os
from osgeo import ogr
driver = ogr.GetDriverByName("ESRI Shapefile")
extfile = '/tmp/rect3.shp'
if os.access( extfile, os.F_OK ):
driver.DeleteDataSource( extfile )
newds = driver.CreateDataSource(extfile)
layernew = newds.CreateLayer('rect3',None,ogr.wkbPolygon)
fieldcnstr = ogr.FieldDefn("关于",ogr.OFTString)
fieldcnstr.SetWidth(36)
fieldf = ogr.FieldDefn("f",ogr.OFTReal)
laydef = layernew.GetLayerDefn()
laydef.AddFieldDefn(fieldcnstr)
laydef.AddFieldDefn(fieldf)
laydef.GetFieldCount()
ERROR 1: OGRFeatureDefn::AddFieldDefn() not allowed on a sealed object ERROR 1: OGRFeatureDefn::AddFieldDefn() not allowed on a sealed object
0
newds.Destroy()
可以看出,输出的图层属性表中并无内容。可能是这个函数只在内存中虚拟地创建了一个字段,与磁盘并无关系,或者说这个函数只是用来创建一个属性字段的框架,仅仅用来定义,在建立新的数据时可能起到参考与引用的作用。在软件处理过程中,由于软件的问题可能会导致图形对象与属性记录并非一一对应,严重时矢量数据可能会损坏。
创建MapInfo
文件时可能会遇到的问题
以TAB格式创建的文件,在写入第一个要素前,需要在相应的驱动上为每个投影设置合理的范围,不过目前尚未有运用OGRDataSource
为新文件设置缺省框架的相应机制,因此在创建新的Layer时,先暂时在mapinfo
驱动上用下面的默认语句来设置框架范围:
- 一个经纬坐标系的框架范围BOUNDS(-180,-90)(180,90)
- 其他的用BOUNDS(-30000000,-15000000)(30000000,15000000)
如果在创建Layer时没有提供坐标系,可以使用以上的投影实例。但如果是地理坐标,那么此例子提供的也只是低精度的地理坐标系。你可以添加-a\_srs wgs84
到ogr2ogr
中,强制地转换成地理模型。
Mapinfo有很多属性限制:
- 只可以创建整数(Integer),实数(Real),字符串(String)字段,其他的类型均不能创建。
- 通常用字符串字段的宽来设置dat文件中的存储大小,因此比字段宽数长的字符串就会被限制。
- 若字符串字段没有分配宽时,其默认值是254。
但是BOUNDS (-30000000,-15000000) (30000000,15000000)这个框架的限制其实并不合理。相对于北京坐标系而言,中国西安的投影坐标系得东移500公里,再加之带号,会有超过30000000的情况,又如福建在3度分带,加之北京坐标系,加之带号差不多在39-40位,加上东移的500000m,最后的坐标值变成39XXXXXX,超出了(-30000000,30000000)的范围。那么用OGR创建的Mapinfo
的福建部分就是错的。但是mif
却可以正常创建,其缺点就是其ACSII字符文件比较大。
from osgeo import ogr
ds = ogr.Open('/data/gdata/GSHHS_c.shp')
layer = ds.GetLayer(0)
spatialRef = layer.GetSpatialRef()
print(spatialRef)
GEOGCS["WGS 84", DATUM["WGS_1984", SPHEROID["WGS 84",6378137,298.257223563, AUTHORITY["EPSG","7030"]], AUTHORITY["EPSG","6326"]], PRIMEM["Greenwich",0, AUTHORITY["EPSG","8901"]], UNIT["degree",0.0174532925199433, AUTHORITY["EPSG","9122"]], AXIS["Latitude",NORTH], AXIS["Longitude",EAST], AUTHORITY["EPSG","4326"]]
Shapefile的投影信息一般存储在.prj文件中,如果没有这个文件,上述函数生成结果为None。
创建投影
建立一个新的Projection的方法如下:首先导入osr库,使用osr.SpatialReference()
创建SpatialReference对象,再向SpatialReference对象导入投影信息:
ImportFromWkt(<wkt>)
ImportFromEPSG(<epsg>)
ImportFromProj4(<proj4>)
ImportFromESRI(<proj_lines>)
ImportFromPCI(<proj>, <units>, <parms>)
ImportFromUSGS(<proj_code>, <zone>)
ImportFromXML(<xml>)
使用下列语句导出Projection字符串:
ExportToWkt()
ExportToPrettyWkt()
ExportToProj4()
ExportToPCI()
ExportToUSGS()
ExportToXML()
根据已有的投影创建新投影
创建空间参考最简单的办法是用ImportFromWkt
导入Wkt中。Wkt获得最简便的方式是安装一个PostgreSQL,再安装PostGIS。在Windows系统中,有安装PostGIS的选项,而Linux系统需要自己编译。PostGIS的spatial_reference
表中包含全套的空间参考,几乎所有的投影坐标信息的Wkt表达全在里面。按需拷贝,开箱即用。
假设已有Wkt,创建一个空间参考:
from osgeo import osr
wkt = spatialRef.ExportToWkt()
spatial = osr.SpatialReference()
spatial.ImportFromWkt(wkt)
0
targetSR = osr.SpatialReference()
targetSR.ImportFromEPSG(4326) #Geo WGS84
coordTrans = osr.CoordinateTransformation(spatial, targetSR)
feature = layer.GetFeature(0)
geom = feature.GetGeometryRef()
geom.ExportToWkt()
geom.Transform(coordTrans)
geom.ExportToWkt()[:120] + ' ... ...'
'POLYGON ((107.544833 76.914028,106.533222 76.499556,111.093417 76.7653890000001,113.934917 75.835,113.583278 75.53119400 ... ...'
targetSR.MorphToESRI()
file = open('/tmp/test.prj', 'w')
file.write(targetSR.ExportToWkt())
file.close()