目录
1. 获取数据的元数据
地理空间元数据,即关于数据的数据,它在地理信息中用于描述地理数据集的内容、质量、表示方式、空间参照系、管理方式以及数据集的其他特征等,帮助和促进人们有效地定位、评价、比较、获取和使用地理相关数据,它是实现地理空间数据集共享的核心内容之一。
可以使用 GetFeatureCount() 函数获取一个图层的要素数量。但只适用于图层,不适用于数据源,因为数据源中的每个图层可以有不同数量的要素、集合类型、空间范围及属性等。
图层的空间范围:上、下、左、右4个方向的最大和最小边界坐标构成的矩形,可以使用GetExtent()函数来获取图层的边界坐标,返回 (min_x, max_x, min_y, max_y) 的数字元组。
例如:
>>> import os
>>> from osgeo import ogr
>>> ds = ogr.Open(r'E:\Google chrome\Download\GIS with python\osgeopy data\osgeopy-data-washington\large_cities.geojson')
>>> lyr = ds.GetLayer(0)
>>> extent = lyr.GetExtent()
>>> print(extent)
(-122.66148376464844, -117.4260482788086, 45.638729095458984, 48.759552001953125)
>>> lyr = ds.GetLayer(0)
>>> extent = lyr.GetExtent()
>>> print(extent)
(-175.22056447761656, 179.21664709402887, -89.99999981438727, 78.21668438639699)
>>> print('Upper left corner: {}, {}'.format(extent[0], extent[3]))
Upper left corner: -175.22056447761656, 78.21668438639699
>>> print('Lower right corner: {}, {}'.format(extent[1], extent[2]))
Lower right corner: 179.21664709402887, -89.99999981438727
获取图层的几何类型,用 GetGeoType() 函数返回一个整数值,不是可理解的字符串。
OGR模块具有多个常量, 可以通过 GetGeoType() 函数将获取的数值与表中的数值进行对比,来检查它属于哪种数据类型。
几何要素类型 | OGR常量 |
---|---|
Point | wkbPoint |
Mulitpoint | wkbMultiPoint |
Line | wkbLineString |
Multiline | wkbMultiLineString |
Polygon | wkbPolygon |
Multipolygon | wkbMultiPolygon |
Unknown geometry type | wkbUnknown |
No geometry | wkbNone |
检查文件类型:
# 交互式
>>> lyr = ds.GetLayer(0)
>>> print(lyr.GetGeomType())
1
>>> print(lyr.GetGeomType() == ogr.wkbPoint) # 与常量进行比较,点数据
True
>>> print(lyr.GetGeomType() == ogr.wkbPolygon) # 与常量进行比较,面数据
False
若图层有多种几何类型,则函数返回 wkbUnknown。
import sys
from osgeo import ogr
fn = r'E:\Google chrome\Download\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)
feat = lyr.GetFeature(0) # 获取图层第一个要素
print(feat.geometry().GetGeometryName()) # 获取几何类型,并返回名字
POINT
1.1 获取空间参考系统
lyr.GetSpatialRef()
>>> print(lyr.GetSpatialRef())
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"]]
1.2 获取图层自身属性字段信息
lyr.schema()
>>> for field in lyr.schema:
print(field.name, field.GetTypeName())
SCALERANK Integer
NATSCALE Integer
LABELRANK Integer
FEATURECLA String
NAME String
NAMEPAR String
<snip>
1.3 OGR库(补充)
2. 矢量数据写入
import sys
from osgeo import ogr
# 打开写入的数据源
ds = ogr.Open(r'E:\Google chrome\Download\global', 1) # 传入参数1,创建新的文件(.shp文件)
if ds is None:
sys.exit('Could not open folder.')
# 获取shaefile文件
in_lyr = ds.GetLayer('ne_50m_populated_places') # 将.shp命名,可以不带后缀
# 若图层存在,则删除该图层
if ds.GetLayer('capital_cities'):
ds.DeleteLayer('capital_cities')
# 创建一个点图层
out_lyr = ds.CreateLayer('capital_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
CreateLayer(name,[srs], [geom_type], [option]):
- name:新建图层名。
- srs:新建图层空间参考系统,默认为空,表示未设置任何空间参考系统。数据如果没有空间参考信息,将难以确定要素在地球上的位置。 如KML只支持WGS84基准面。
- geom_type:几何常量类型,表示图层中要素的几何类型,默认为wkbUnknown。
- options:只适用于特定的矢量数据类型。
基于已有.shp文件,创建一个新的.shp文件,并使用相同的空间参考系统。
import sys
from osgeo import ogr
ds = ogr.Open(r'E:\Google chrome\Download\global',1)
if ds is None:
sys.exit('Could not open folder.')
in_lyr = ds.GetLayer('ne_50m_populated_places') # 获取shapefile文件
if ds.GetLayer('captial_cities'):
ds.DeleteLayer('captial_cities')
# 创建一个点图层
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.SyncToDisk
2.1 创建新数据源
创建新数据源的关键是正确使用驱动程序,且一种驱动程序只能处理一种类型的矢量数据。如使用GeoJSON驱动程序创建一个shp文件任然是创建一个GeoJSON文件。
获取驱动程序
方法一: 从已打开的数据集中获取
from osgeo import ogr
ds = ogr.Open(r'E:\Google chrome\Download\california\california_50m_wind_albers.shp')
driver = ds.GetDriver()
print(driver.name)
ESRI Shapefile
方法二:使用 OGR中GetDriverByName()函数,传递给它驱动程序的简称。可以使用:ogrinfo --formats 获得驱动程序。
import sys
from osgeo import ogr
esri_fn = r'E:\Google chrome\Download\california'
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'E:\Google chrome\Download\global',1)
if ds is None:
sys.exit('Could not open folder.')
in_lyr = ds.GetLayer('ne_50m_populated_places')
if esri_ds.GetLayer('captial_cities'):
esri_ds.DeleteLayer('captial_cities')
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
数据源和图层创建选项是完全不同的,使用数据源选项创建一个SpatialLite数据源,而不使用SQLite:
import sys
from osgeo import ogr
json_fn = r'E:\Google chrome\Download\global'
json_driver = ogr.GetDriverByName('GeoJSON')
json_ds = json_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'E:\Google chrome\earth.sqlite',['SPATIALITE=yes'])
注意:创建鑫都数据源时,不能覆盖现有的数据源。如果可能会覆盖,需要在创建之前删除就数据。所以要使用驱动程序删除现有数据源,而不能使用python内置函数。因为驱动程序可以确保所有必要的文件都被删除掉,如shp、dbf、shx等。
方法:
if os.path.exists(json_fn):
json_driver.DeleteDataSource(json_fn)
json_ds = json_driver.CreateDataSource(json_fn)
if json_ds is None:
sys.exit('Could not create {0}.'.format(json_fn))
OGR例外,ogr.UseExceptions() 将在程序错误的地方自动报错。
import sys
import os
from osgeo import ogr
data_dir = r'E:\Google chrome\Download\california'
# 开启例外
ogr.UseExceptions()
fn = os.path.join(data_dir, 'output', 'natural_earth_50m.sqlite')
driver = ogr.GetDriverByName('SQLite')
print('Doing some preliminary analysis...')
try:
# 如果文件已经存在,这里将会报错
ds = driver.CreateDataSource(fn)
lyr = ds.CreateLayer('layer')
# 试着保存一部分数据
except RuntimeError as e:
# 输出错误信息,并继续
print(e)
print('Doing some more analysis...')
2.2 新建属性字段
通过FieldDefn构造函数提供名字和数据类型来创建一个字段。
字段数据类型 | OGR常量 |
---|---|
Integer | OFTInteger |
List of integers | OFTIntegerList |
Floating point number | OFTReal |
String | OFTString |
List of String | OFTStringList |
Date | OFTDate |
Date and time | OFTDateTime |
创建一个基本字段定义后,使用它添加一个字段到图层前,并且可以添加一些其他的限制,比如浮点精度或者字段宽度。
在shapefile文件中设置字段精度,必须设置字段宽度。
import sys
import os
from osgeo import ogr
data_dir = r'E:\Google chrome\Download\GIS with python\osgeopy data\osgeopy-data-washington\osgeopy-data'
in_fn = os.path.join(data_dir, 'Washington', 'large_cities.shp')
in_ds = ogr.Open(in_fn, 0)
if in_ds is None:
sys.exit('Could not open {0}.'.format(in_fn))
in_lyr = in_ds.GetLayer(0)
# 创建输出shp文件
driver = in_ds.GetDriver()
out_fn = os.path.join(data_dir, 'output', 'precision_test.shp')
if os.path.exists(out_fn):
driver.DeleteDataSource(out_fn)
out_ds = driver.CreateDataSource(out_fn)
if out_ds is None:
sys.exit('Could not create {0}.'.format(out_fn))
# 创建一个shp图层
out_lyr = out_ds.CreateLayer('precision_test',
in_lyr.GetSpatialRef(),
ogr.wkbPoint)
# 将name字段设置为宽度为6,但它将被展开。
name_fld = ogr.FieldDefn('name', ogr.OFTString)
name_fld.SetWidth(6)
out_lyr.CreateField(name_fld)
# 使用默认精度创建两个属性字段
coord_fld = ogr.FieldDefn('X_default', ogr.OFTReal)
out_lyr.CreateField(coord_fld)
coord_fld.SetName('Y_default')
out_lyr.CreateField(coord_fld)
# 使用更小的精度创建两个属性字段
# 创建并添加第一个属性字段
coord_fld = ogr.FieldDefn('X_short', ogr.OFTReal)
coord_fld.SetWidth(8) # 设置字段宽度
coord_fld.SetPrecision(3) # 设置字段精度
# 重用FieldDefn来创建第二个字段
# 这里创建两个相同的字段定义对象
# 使创建更加容易
out_lyr.CreateField(coord_fld)
coord_fld.SetName('Y_short')
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.SetName('Y_short')
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)
3. 更新数据
3.1 改变图层定义
根据数据类型的不同,可以添加、删除、更改字段来编辑图层定义。对于一个满意的字段定义,可以使用AlterFiledDefn函数用新字段替换现有字段:
AlterFiledDefn(iField,filed_def, nFlags)
- AlterFiledDefn是改变的属性字段对应的索引
- filed_def是新属性字段的定义对象
- nFlags是一个整数,是表中一个或多个常量的总和
需要更改的字段属性 | 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 |
步骤 :创建一个包含新的属性的字段定义,找到现有字段索引,确定要使用的OGR常量。如将“Name”改为“CITY_NAME”:
首先,查看原有的数据集。
import sys
import os
import ospybook as pb
from osgeo import ogr
data_dir = r'E:\Google chrome\Download\GIS with python\osgeopy data\osgeopy-data-washington\osgeopy-data'
original_fn = os.path.join(data_dir, 'Washington', 'large_cities.shp')
new_fn = os.path.join(data_dir, 'Washington', 'large_cities2.shp')
# 复制.shp文件
pb.copy_datasource(original_fn, new_fn)
# 打开复制的shp文件准备写入
ds = ogr.Open(new_fn, 1)
if ds is None:
sys.exit('Could not open {0}.'.format(new_fn))
lyr = ds.GetLayer(0)
# 在做改变之前,先查看数据集
print('Original attributes')
pb.print_attributes(lyr, geom=False)
然后,改变字段属性
# 获取图层中Name的所有相关属性信息
# 创建新的文件
# 定义并使用改变新的文件属性
i = lyr.GetLayerDefn().GetFieldIndex('Name')
fld_defn = ogr.FieldDefn('City_Name', ogr.OFTString)
lyr.AlterFieldDefn(i, fld_defn, ogr.ALTER_NAME_FLAG)
# 改变多个字段的多个属性信息
# 改变POINT_X为X_coord,并设置精度为4
# 创建新的字段定义时,使用的是原始字段宽度。如果设置宽度不是足够大,结果将错误
# 所以需要使用与原始数据一样的字段宽度,这样会比较安全
# 所以字段必须重写
lyr_defn = lyr.GetLayerDefn()
i = lyr_defn.GetFieldIndex('X')
width = lyr_defn.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)
3.2 添加、更新和删除要素
添加要素:首先基于字段创建一个空要素,填充它,然后,把它插入到该图层中去。
# 为图层中的每一个要素添加一个唯一的ID
lyr.ResetReading()
lyr.CreateField(ogr.FieldDefn('ID', ogr.OFTInteger))
n = 1
for feat in lyr:
feat.SetField('ID', n)
lyr.SetFeature(feat)
n += 1
print('\nID has been added and precision has taken effect')
pb.print_attributes(lyr, geom=False)
删除要素:
lyr.ResetReading()
for feat in lyr:
if feat.GetField('City_Name') == 'Seattle':
lyr.DeleteFeature(feat.GetFID())
print('\nSeattle deleted')
pb.print_attributes(lyr, geom=False)
回收空间:
ds.ExecuteSQL('REPACK ' + lyr.GetName())
ds.ExecuteSQL('RECOMPUTE EXTENT ON ' + lyr.GetName())
print('\nDatabase packed')
pb.print_attributes(lyr, geom=False)
4. 例子
1.使用默认精度创建json文件
import sys
import os
from osgeo import ogr
data_dir = r'E:\Google chrome\Download'
shp_fn = os.path.join(data_dir, 'global', 'ne_50m_admin_0_countries.shp')
shp_ds = ogr.Open(shp_fn, 0)
if shp_ds is None:
sys.exit('Could not open {0}'.format(shp_fn))
shp_lyr = shp_ds.GetLayer(0)
json_driver = ogr.GetDriverByName('GeoJSON')
# 使用默认精度创建一个json文件
# 创建数据源
json_fn = os.path.join(data_dir, 'GIS with python', 'africa-default.geojson')
if os.path.exists(json_fn):
json_driver.DeleteDataSource(json_fn)
json_ds = json_driver.CreateDataSource(json_fn)
if json_ds is None:
sys.exit('Could not create {0}.'.format(json_fn))
# 创建无选择项的图层
json_lyr = json_ds.CreateLayer('africa',
shp_lyr.GetSpatialRef(),
ogr.wkbMultiPolygon)
# 写入数据
shp_lyr.ResetReading()
json_feat = ogr.Feature(json_lyr.GetLayerDefn())
for shp_feat in shp_lyr:
if shp_feat.GetField('CONTINENT') == 'Africa':
json_feat.SetGeometry(shp_feat.geometry())
json_lyr.CreateFeature(json_feat)
del json_ds
结果:
2.使用精度为6,创建json文件
import sys
import os
from osgeo import ogr
data_dir = r'E:\Google chrome\Download'
shp_fn = os.path.join(data_dir, 'global', 'ne_50m_admin_0_countries.shp')
shp_ds = ogr.Open(shp_fn, 0)
if shp_ds is None:
sys.exit('Could not open {0}'.format(shp_fn))
shp_lyr = shp_ds.GetLayer(0)
json_driver = ogr.GetDriverByName('GeoJSON')
# 使用精度为6,创建一个json文件
# 用选择项COORDINATE_PRECISION创建一个json文件,并且设置精度为为6
# 创建数据源
json_fn = os.path.join(data_dir, 'GIS with python', 'africa-6digit.geojson')
if os.path.exists(json_fn):
json_driver.DeleteDataSource(json_fn)
json_ds = json_driver.CreateDataSource(json_fn)
if json_ds is None:
sys.exit('Could not create {0}.'.format(json_fn))
lyr_options = ['COORDINATE_PRECISION=6']
json_lyr = json_ds.CreateLayer('africa',
shp_lyr.GetSpatialRef(),
ogr.wkbMultiPolygon,
lyr_options)
# 写入数据
shp_lyr.ResetReading()
json_feat = ogr.Feature(json_lyr.GetLayerDefn())
for shp_feat in shp_lyr:
if shp_feat.GetField('CONTINENT') == 'Africa':
json_feat.SetGeometry(shp_feat.geometry())
json_lyr.CreateFeature(json_feat)
json_ds
结果:
3.Bounding box
import sys
import os
from osgeo import ogr
data_dir = r'E:\Google chrome\Download'
shp_fn = os.path.join(data_dir, 'global', 'ne_50m_admin_0_countries.shp')
shp_ds = ogr.Open(shp_fn, 0)
if shp_ds is None:
sys.exit('Could not open {0}'.format(shp_fn))
shp_lyr = shp_ds.GetLayer(0)
json_driver = ogr.GetDriverByName('GeoJSON')
# 使用选择项COORDINATE_PRECISION 和 WRITE_BBOX 创建一个json文件
json_fn = os.path.join(data_dir, 'GIS with python', 'africa-bbox.geojson')
if os.path.exists(json_fn):
json_driver.DeleteDataSource(json_fn)
json_ds = json_driver.CreateDataSource(json_fn)
if json_ds is None:
sys.exit('Could not create {0}.'.format(json_fn))
lyr_options = ['COORDINATE_PRECISION=6', 'WRITE_BBOX=YES']
json_lyr = json_ds.CreateLayer('africa',
shp_lyr.GetSpatialRef(),
ogr.wkbMultiPolygon,
lyr_options)
# Write some data.
shp_lyr.ResetReading()
json_feat = ogr.Feature(json_lyr.GetLayerDefn())
for shp_feat in shp_lyr:
if shp_feat.GetField('CONTINENT') == 'Africa':
json_feat.SetGeometry(shp_feat.geometry())
json_lyr.CreateFeature(json_feat)
del json_ds
结果:
4.创建一个新的例子
import sys
import os
from osgeo import ogr
data_dir = r'E:\Google chrome\Download'
shp_fn = os.path.join(data_dir, 'global', 'ne_50m_admin_0_countries.shp')
json_fn = os.path.join(data_dir, 'GIS with python', 'africa.geojson')
# 打开输入文件
shp_ds = ogr.Open(shp_fn, 0)
if shp_ds is None:
sys.exit('Could not open {0}'.format(shp_fn))
shp_lyr = shp_ds.GetLayer(0)
# 创建输出文件
json_driver = ogr.GetDriverByName('GeoJSON')
if os.path.exists(json_fn):
json_driver.DeleteDataSource(json_fn)
json_ds = json_driver.CreateDataSource(json_fn)
if json_ds is None:
sys.exit('Could not create {0}.'.format(json_fn))
lyr_options = ['COORDINATE_PRECISION=6']
json_lyr = json_ds.CreateLayer('africa',
shp_lyr.GetSpatialRef(),
ogr.wkbMultiPolygon,
lyr_options)
# 在图层中添加属性文件
name_fld = ogr.FieldDefn('Name', ogr.OFTString)
json_lyr.CreateField(name_fld)
pop_fld = ogr.FieldDefn('Population', ogr.OFTInteger)
json_lyr.CreateField(pop_fld)
# 添加整型文件,及错误的添加字符串
test_fld = ogr.FieldDefn('Test_field', ogr.OFTInteger)
json_lyr.CreateField(test_fld)
# 获取属性定义
feat_defn = json_lyr.GetLayerDefn()
# 创建输出特征并反复使用
json_feat = ogr.Feature(feat_defn)
for shp_feat in shp_lyr:
if shp_feat.GetField('CONTINENT') == 'Africa':
# 复制属性值到文件中
name = shp_feat.GetField('NAME')
pop = shp_feat.GetField('POP_EST')
json_feat.SetField('Name', name)
json_feat.SetField('Population', pop)
# 在整型文件中放置字符串
json_feat.SetField('Test_field', name)
# 赋值几何形状
json_feat.SetGeometry(shp_feat.geometry())
# 将数据输入到geojson文件中
json_lyr.CreateFeature(json_feat)
del json_ds, shp_ds
结果: