OGR
简单要素库是地理空间数据抽象库(Geospatial Data Abstraction Library,GDAL)的一部分,它是一个非常流行的用于读写空间数据的开源库
环境
python 2.7
安装
conda install gdal pyproj
读取矢量数据 shapefile
import sys
from osgeo import ogr
fn = r'C:/Users/miko/Documents/Book_code/python地理数据处理_数据集/python地理数据处理_数据集/osgeopy-data/global/ne_50m_populated_places.shp'
ds = ogr.Open(fn,0)
if ds is None:
sys.exit('Could not open {0}.'.format(fn))
lyr = ds.GetLayer(0)
i=0
for feat in lyr:
pt = feat.geometry()
x = pt.GetX()
y = pt.GetY()
name = feat.GetField('NAME')
pop = feat.GetField('POP_MAX')
print(name,pop,x,y)
i += 1
if i == 10:
break
del ds
打开数据源
ogr.open()
通过Open函数传递文件名和一个可选更新标记,打开一个数据源。这是OGR模块中的一个独立函数,需要事先用这个模块名称定义一个函数名,以便Python可以找到它。
如果第二个变量不默认设置为0,将以只读的模式打开文件。可以通过将其设置为1或者True来以编辑的模式打开
如果shapefile文件无法打开,则Open函数返回为空
获取图层信息
ds.GetLayer()
一个数据源由一个或多个图层组成,所以打开数据源后,需要获得具体的图层
可以通过该函数来获取图层索引或者图层名,进而得到这个数据集中对应的图层
0索引
如果不给GetLayer函数提供任何参数,默认返回这个数据源中的第一个图层
获取具体图层中要素的信息
图层由一个或多个要素组成,每一个要素表示一个地理对象。地理对象的几何形状和属性值都存储在这些要素中。
for feat in lyr:
pt = feat.geometry()
x = pt.GetX()
y = pt.GetY()
name = feat.GetField('NAME')
pop = feat.GetField('POP_MAX')
feat.geometry()
获取要素的几何形状,从而获取其坐标信息
feat.GetField()
通过属性名或索引值获取对应字段的值,该函数返回的数据类型与底层数据集中的数据类型一致。补充函数:feat.GetfieldAsString()
访问特定要素
import sys
from osgeo import ogr
fn = r'C:/Users/miko/Documents/Book_code/python地理数据处理_数据集/python地理数据处理_数据集/osgeopy-data/global/ne_50m_populated_places.shp'
ds = ogr.Open(fn,0)
if ds is None:
sys.exit('Could not open {0}.'.format(fn))
lyr = ds.GetLayer(0)
num_features = lyr.GetFeatureCount() #获取要素总数
last_feature = lyr.GetFeature(num_features-1) #获取图层中的最后一个要素
print(num_features)
print(last_feature.name)
lyr.GetFeatureCount()
获取要素的总数
lyr.GetFeature(FIDs)
获取FID号对应的要素
ospybook 模块
环境 python2.7
安装
pip install https://github.com/cgarrard/osgeopy-code/raw/master/ospybook-latest.zip
查看属性
import ospybook as pb
fn = r'C:/Users/miko/Documents/Book_code/python地理数据处理_数据集/python地理数据处理_数据集/osgeopy-data/global/ne_50m_populated_places.shp'
pb.print_attributes(fn,3,['NAME','POP_MAX'])
print_attributes(lyr_or_fn,[n],[fields],[geom],[reset])
该函数将属性值信息输出到屏幕上
这个函数对于查看小数据量的属性信息运行正常,但如果使用它来输出一个大数据的所有属性,你会后悔的
参数 | 含义 |
---|---|
lyr_or_fn | 可以是一个图层,也可以是一个数据源路径 如果是一个数据源,使用的是第一个图层 |
n | 是一个可选值 用于设置输出记录数目,默认输出所有数值 |
fields | 是一个可选值 用于设置输出结果中包含的属性字段列表,默认包括所有字段 |
geom | 是一个可选的布尔值 用于设置是否输出集合要素类型,默认为True |
reset | 是一个可选的布尔值 用于设置在输出数值之前,是否重置到第一条记录,默认为True |
绘制空间数据
需创建VectorPlotter
类的新实例,需为构造函数传递一个布尔变量
import os
os.chdir('C:\Users\miko\Documents\Book_code\osgeopy-data\global') # 更改工作目录,最好是英文路径
from ospybook.vectorplotter import VectorPlotter
vp = VectorPlotter(True)
vp.plot('ne_50m_admin_0_countries.shp',fill=False) # fill参数使shp文件用一个空心多边形表示
vp.plot('ne_50m_populated_places.shp','bo') # bo符号表示蓝色的圆圈
vp.draw()
ArcGIS 桌面版本打开效果:
plot(self,geom_or_lyr,[symbol],[name],[kwargs])
参数 | 含义 |
---|---|
geom_or_lyr | 可以是一个要素对象、图层或者指向一个数据源的路径。 如果是一个数据源,数据源中的第一个图层会被绘制显示 |
symbol | 是一个可选值 用于设置几何元素的符号样式 |
name | 是一个可选值 用于为数据设定名称、以便后期可以访问它 |
kwargs | 是一个可选值 通过关键字进行指定 kwargs经常被用作一个不确定数量的关键字参数(an indeterminate number of keyword arguments)的缩写 |
获取数据的元数据
获取空间范围
图层的空间范围是从上、下、左、右4个方向的最大和最小边界坐标构成的矩形
lyr.GetExtent()
获取一个图层对象的边界坐标,它会返回一个类似(min_x,max_x,min_y,max_y)的元组
import sys
from osgeo import ogr
fn = r'C:/Users/miko/Documents/Book_code/osgeopy-data/Washington/large_cities.geojson'
ds = ogr.Open(fn,0)
if ds is None:
sys.exit('Could not open {0}.'.format(fn))
lyr = ds.GetLayer(0)
extent = lyr.GetExtent()
print(extent)
获取几何类型
lyr.GetGeomType()
返回的是一个整数值,而不是一个可理解的字符串
如果这个图层有多种几何类型,例如点和面,GetGeomType函数就会返回wkbUnknown
在ORG几何要素常量中的WKB前缀为熟知二进制(Well-KnowingBinary,WKB)的缩写,是用于不同软件程序间进行几何要素类型转换的一种二进制表示标准。因为它是二进制格式,所以人们无法直接阅读获得其表示的内容,但是熟知文本格式(Well-Known Text,WKT)可以直接阅读
import sys
from osgeo import ogr
fn = r'C:/Users/miko/Documents/Book_code/osgeopy-data/Washington/large_cities.shp'
ds = ogr.Open(fn,0)
if ds is None:
sys.exit('Could not open {0}.'.format(fn))
lyr = ds.GetLayer(0)
print(lyr.GetGeomType())
print(lyr.GetGeomType() == ogr.wkbPoint) #与常量进行比较,获取数据类型
print(lyr.GetGeomType() == ogr.wkbPolygon)
几何要素类型 | OGR常量 |
---|---|
Point | wkbPoint |
Mulitpoint | wkbMultiPoint |
Line | wkbLineString |
Multiline | wkbMultiLineString |
Polygon | wkbPolygon |
Multipolygon | wkbMultiPolygon |
Unknown geometry type | wkbUnknown |
No geometry | wkbNone |
获取一个几何类型的可读字符串,可以从要素几何类型中进行获取
import sys
from osgeo import ogr
fn = r'C:/Users/miko/Documents/Book_code/osgeopy-data/Washington/large_cities.shp'
ds = ogr.Open(fn,0)
if ds is None:
sys.exit('Could not open {0}.'.format(fn))
lyr = ds.GetLayer(0)
feat = lyr.GetFeature(0)
print(feat.geometry().GetGeometryName())
获取空间参考系统
lyr.GetSpatialRef()
得到用WKT格式描述空间参考系统的字符串
import sys
from osgeo import ogr
fn = r'C:/Users/miko/Documents/Book_code/osgeopy-data/Washington/large_cities.shp'
ds = ogr.Open(fn,0)
if ds is None:
sys.exit('Could not open {0}.'.format(fn))
lyr = ds.GetLayer(0)
print(lyr.GetSpatialRef())
获取图层自身属性字段信息
lyr.schema
import sys
from osgeo import ogr
fn = r'C:/Users/miko/Documents/Book_code/osgeopy-data/Washington/large_cities.shp'
ds = ogr.Open(fn,0)
if ds is None:
sys.exit('Could not open {0}.'.format(fn))
lyr = ds.GetLayer(0)
for field in lyr.schema:
print(field.name,field.GetTypeName())
矢量数据写入
基于已有的shp文件创建新的shapefile文件
import sys
from osgeo import ogr
ds = ogr.Open(r'C:\Users\miko\Documents\Book_code\osgeopy-data\global',1)
# 打开一个文件夹作为数据源,而不是使用shapefile文件
# 如果文件夹中大多数文件是shapefile格式文件,每个shapefile文件会被视为一个图层
# 参数 1,表示可以在文件夹中创建一个新的图层 (shapefile文件)
if ds is None:
sys.exit('Could not open folder.')
in_lyr = ds.GetLayer('ne_50m_populated_places') #获取shapefile文件
#将没有后缀.shp,传递给GetLayer函数来获取居住地shapefile文件作为一个图层
if ds.GetLayer('captial_cities'):
ds.DeleteLayer('captial_cities')
# 因为OGR不会覆盖现有的图层,所以需要检查现有的图层是否已经存在,如果存在,就删掉它
#创建一个点图层
out_lyr = ds.CreateLayer('captial_cities',in_lyr.GetSpatialRef(),ogr.wkbPoint)
out_lyr.CreateFields(in_lyr.schema)
#创建一个空要素
out_defn = out_lyr.GetLayerDefn()
out_feat = ogr.Feature(out_defn)
#复制几何要素和属性
for in_feat in in_lyr:
if in_feat.GetField('FEATURECLA') == 'Admin-0 capital':
geom = in_feat.geometry()
out_feat.SetGeometry(geom)
for i in range(in_feat.GetFieldCount()):
value = in_feat.GetField(i)
out_feat.SetField(i,value)
out_lyr.CreateFeature(out_feat)
del ds
新创建的为图中棕色部分点
创建新图层
ds.CreateLayer(name,[srs],[geom_type],[options])
创建一个新图层,用来存储输出的数据
参数 | 含义 |
---|---|
name | 表示新建的图层名称 |
srs | 表示新建的图层使用的空间参考系统 默认为空,表示没有设置任何空间参考系统 |
geom_type | 用来表示图层中要素的几何类型 默认为wkbUnknown |
options | 是一个图层创建时的选项列表,只适合于特定的矢量数据类型 |
创建空要素
out_defn = out_lyr.GetLayerDefn()
创建一个要素需要获得要素定义,其包含几何类型和所有属性字段的信息,这样创建的空要素才具有相同的属性字段和几何类型
你需要获得待添加要素的图层的要素定义信息,而且要在图层做了添加、删除或更新了任何字段之后进行获取。 如果先获取要素的定义信息,然后属性字段做了任何改变,获得的信息就过时了,这意味着,基于这些过时的信息插入的要素,会与实际不符,程序崩溃
out_feat = ogr.Feature(out_defn)
为新创建的图层添加要素,需要创建一个用于存储几何要素和属性的虚拟要素,然后将其插入到图层中
要素插入图层
out_lyr.CreateFeature()
将要素插入到图层中,这个函数会在图层中保存一个要素的副本,包括添加的所有信息
图层的保存
del ds
删除ds变量强制关闭文件,并将所有的编辑成果写入到磁盘中
ds.SyncToDisk()
保持数据源打开,但编辑的成果写入磁盘中
创建新的数据源
获取驱动
方式一
是从已经打开的数据集中进行获取,这将允许创建一个新的数据源,其和已存在的数据源矢量数据格式一样
from osgeo import ogr
ds = ogr.Open(r'C:/Users/miko/Documents/Book_code/osgeopy-data/global/ne_50m_admin_0_countries.shp')
driver = ds.GetDriver()
print(driver.name)
from osgeo import ogr
ds = ogr.Open(r'C:/Users/miko/Documents/Book_code/osgeopy-data/Washington/large_cities.geojson')
driver = ds.GetDriver()
print(driver.name)
方式二
使用OGR中GetDriverByName()
函数,传递给它驱动程序的简称
这些驱动程序名称都在OGR网站上有介绍
from osgeo import ogr
esri_fn = r'C:/Users/miko/Documents/Book_code/osgeopy-data/test'
esri_driver = ogr.GetDriverByName('ESRI Shapefile')
esri_ds = esri_driver.CreateDataSource(esri_fn)
if esri_ds is None:
sys.exit('Could not create {0}.'.format(esri_fn))
新创建的数据源会打开等待写入,剩下就可以为数据源添加图层和要素信息
创建新的数据源,不能覆盖现有的数据源
from osgeo import ogr
esri_fn = r'C:/Users/miko/Documents/Book_code/osgeopy-data/test'
esri_driver = ogr.GetDriverByName('ESRI Shapefile')
esri_ds = esri_driver.CreateDataSource(esri_fn)
if esri_ds is None:
sys.exit('Could not create {0}.'.format(esri_fn))
ds = ogr.Open(r'C:\Users\miko\Documents\Book_code\osgeopy-data\global',1)
# 打开一个文件夹作为数据源,而不是使用shapefile文件
# 如果文件夹中大多数文件是shapefile格式文件,每个shapefile文件会被视为一个图层
# 参数 1,表示可以在文件夹中创建一个新的图层 (shapefile文件)
if ds is None:
sys.exit('Could not open folder.')
in_lyr = ds.GetLayer('ne_50m_populated_places') #获取shapefile文件
#将没有后缀.shp,传递给GetLayer函数来获取居住地shapefile文件作为一个图层
if esri_ds.GetLayer('captial_cities'):
esri_ds.DeleteLayer('captial_cities')
# 因为OGR不会覆盖现有的图层,所以需要检查现有的图层是否已经存在,如果存在,就删掉它
#创建一个点图层
out_lyr = esri_ds.CreateLayer('captial_cities',in_lyr.GetSpatialRef(),ogr.wkbPoint)
out_lyr.CreateFields(in_lyr.schema)
#创建一个空要素
out_defn = out_lyr.GetLayerDefn()
out_feat = ogr.Feature(out_defn)
#复制几何要素和属性
for in_feat in in_lyr:
if in_feat.GetField('FEATURECLA') == 'Admin-0 capital':
geom = in_feat.geometry()
out_feat.SetGeometry(geom)
for i in range(in_feat.GetFieldCount()):
value = in_feat.GetField(i)
out_feat.SetField(i,value)
out_lyr.CreateFeature(out_feat)
del esri_ds
使用其他驱动创建数据源
from osgeo import ogr
json_fn = r'C:/Users/miko/Documents/Book_code/osgeopy-data/test'
json_driver = ogr.GetDriverByName('GeoJSON')
json_ds = esri_driver.CreateDataSource(json_fn)
if json_ds is None:
sys.exit('Could not create {0}.'.format(json_fn))
from osgeo import ogr
driver = ogr.GetDriverByName('SQLite')
ds = driver.CreateDataSource(r'D:/osgeopy-data/global/earth.sqlite',['SPATIALITE=yes'])
新建属性字段
给新建的数据源、图层,创建属性
from osgeo import ogr
in_ds = ogr.Open(r'D:\WorkSpace\Jupyter notebook\data\osgeopy-data-washington\osgeopy-data\Washington',1)
out_esri_fn = r'D:\WorkSpace\Jupyter notebook\data\test'
out_esri_driver = ogr.GetDriverByName('ESRI Shapefile')
out_esri_ds = out_esri_driver.CreateDataSource(esri_fn)
in_lyr = in_ds.GetLayer('large_cities')
print(in_lyr.GetSpatialRef())
if out_esri_ds.GetLayer("test_create_lyr"):
out_esri_ds.DeleteLayer("test_create_lyr")
out_lyr = out_esri_ds.CreateLayer("test_create_lyr",in_lyr.GetSpatialRef(),ogr.wkbPoint)
coord_fld = ogr.FieldDefn("name",ogr.OFTString)
coord_fld.SetWidth(6)
out_lyr.CreateField(coord_fld)
coord_fld = ogr.FieldDefn("X_default",ogr.OFTReal)
out_lyr.CreateField(coord_fld)
coord_fld = ogr.FieldDefn("Y_default",ogr.OFTReal)
out_lyr.CreateField(coord_fld)
coord_fld = ogr.FieldDefn("X_precision",ogr.OFTReal)
coord_fld.SetWidth(8)
coord_fld.SetPrecision(3)
out_lyr.CreateField(coord_fld)
coord_fld = ogr.FieldDefn("Y_precision",ogr.OFTReal)
coord_fld.SetWidth(8)
coord_fld.SetPrecision(3)
out_lyr.CreateField(coord_fld)
del out_esri_ds
字段类型
字段数据类型 | OGR常量 |
---|---|
Integer | OFTInteger |
List of integers | OFTIntegerList |
Floating point number | OFTReal |
List of floating point numbers | OFTRealList |
String | OFTString |
List of String | OFTStringList |
Date | OFTDate |
Time of day | OFTTime |
Date and time | OFTDateTime |
定义属性字段的简单写法
coord_fld = ogr.FieldDefn('X_default', ogr.OFTReal)
out_lyr.CreateField(coord_fld)
coord_fld.SetName('Y_default')
out_lyr.CreateField(coord_fld)
给属性表赋值
from osgeo import ogr
in_ds = ogr.Open(r'D:\WorkSpace\Jupyter notebook\data\osgeopy-data-washington\osgeopy-data\Washington',1)
out_esri_fn = r'D:\WorkSpace\Jupyter notebook\data\test'
out_esri_driver = ogr.GetDriverByName('ESRI Shapefile')
out_esri_ds = out_esri_driver.CreateDataSource(esri_fn)
in_lyr = in_ds.GetLayer('large_cities')
print(in_lyr.GetSpatialRef())
if out_esri_ds.GetLayer("test_create_lyr"):
out_esri_ds.DeleteLayer("test_create_lyr")
out_lyr = out_esri_ds.CreateLayer("test_create_lyr",in_lyr.GetSpatialRef(),ogr.wkbPoint)
coord_fld = ogr.FieldDefn("name",ogr.OFTString)
coord_fld.SetWidth(6)
out_lyr.CreateField(coord_fld)
coord_fld = ogr.FieldDefn("X_default",ogr.OFTReal)
out_lyr.CreateField(coord_fld)
coord_fld = ogr.FieldDefn("Y_default",ogr.OFTReal)
out_lyr.CreateField(coord_fld)
coord_fld = ogr.FieldDefn("X_Short",ogr.OFTReal)
coord_fld.SetWidth(8)
coord_fld.SetPrecision(3)
out_lyr.CreateField(coord_fld)
coord_fld = ogr.FieldDefn("Y_Short",ogr.OFTReal)
coord_fld.SetWidth(8)
coord_fld.SetPrecision(3)
out_lyr.CreateField(coord_fld)
out_feat = ogr.Feature(out_lyr.GetLayerDefn())
for in_feat in in_lyr:
pt = in_feat.geometry()
name = in_feat.GetField("name")
out_feat.SetGeometry(in_feat.geometry())
out_feat.SetField("name",name)
out_feat.SetField("X_default",pt.GetX())
out_feat.SetField("Y_default",pt.GetY())
out_feat.SetField("X_Short",pt.GetX())
out_feat.SetField("Y_Short",pt.GetY())
out_lyr.CreateFeature(out_feat)
del out_esri_ds
更改图层属性字段定义
需要更改的字段属性 | OGR常量 |
---|---|
Field name only | ALTER_NAME_FLAG |
Field type only | ALTER_TYPE_FLAG |
Field width and/or precision only | ALTER_WIDTH_PRECISION_FLAG |
All of the above | ALTER_ALL_FLAG |
更新字段的单个属性信息
from osgeo import ogr
ds = ogr.Open(r'D:\WorkSpace\Jupyter notebook\data\test',1)
lyr = ds.GetLayer("test_create_lyr")
i = lyr.GetLayerDefn().GetFieldIndex("name")
print(i)
fld_defn = ogr.FieldDefn('City_Name',ogr.OFTString)
lyr.AlterFieldDefn(i,fld_defn,ogr.ALTER_NAME_FLAG)
del ds
更新字段的多个属性信息
from osgeo import ogr
ds = ogr.Open(r'D:\WorkSpace\Jupyter notebook\data\test',1)
lyr = ds.GetLayer("test_create_lyr")
i = lyr.GetLayerDefn().GetFieldIndex("X_default")
print(i)
width = lyr.GetLayerDefn().GetFieldDefn(i).GetWidth()
fld_defn = ogr.FieldDefn('X_coord',ogr.OFTReal)
fld_defn.SetWidth(width)
fld_defn.SetPrecision(4)
flag = ogr.ALTER_NAME_FLAG + ogr.ALTER_WIDTH_PRECISION_FLAG
lyr.AlterFieldDefn(i,fld_defn,flag)
del ds
复制数据源与文件
import os
import ospybook as pb
original_fn = r'D:\WorkSpace\Jupyter notebook\data\osgeopy-data-washington\osgeopy-data\Washington\large_cities.shp'
new_fn = r'D:\WorkSpace\Jupyter notebook\data\test\large_cities2.shp'
pb.copy_datasource(original_fn, new_fn)
查看图层属性
import os
import ospybook as pb
original_fn = r'D:\WorkSpace\Jupyter notebook\data\osgeopy-data-washington\osgeopy-data\Washington\large_cities.shp'
new_fn = r'D:\WorkSpace\Jupyter notebook\data\test\large_cities2.shp'
pb.copy_datasource(original_fn, new_fn)
# Open the copied shapefile for writing.
ds = ogr.Open(new_fn, 1)
if ds is None:
sys.exit('Could not open {0}.'.format(new_fn))
lyr = ds.GetLayer(0)
# Take a look at the attributes before you change anything.
print('Original attributes')
pb.print_attributes(lyr, geom=True)