矢栅处理(杂)：OSM数据下载及shapefile读写

发表于 分类于 阅读次数： Waline： 本文字数： 3.3k 阅读时长 ≈ 3 分钟

OSM地图数据是志愿地理信息数据(VGI)的代表，更新快(天)覆盖区域广，在大区域研究中其数据尤其是路网应用广泛。这里简单介绍一下数据下载渠道及osm格式转换即shapefile的通用读写代码。

涉及多个国家的研究可能没有其他可替代的全球数据产品，所以这个数据应用很广。官网下载并不方便，可用通过Geofabrik下载对应区域和国家的shape文件，其中中国的shapefile文件是完全开放的：

  • 下载路径，点击[raw directory index]可以下载历年的数据
  • 但是可能会发现有的国家并没有提供，仅提供了.som.pbf数据或者bz2的数据
  • shapefile的文件数据是按照数据类型（poi、路网、土地利用等等提供，是转换好的），但.osm的数据则是原始格式，pbf直接在qgis打开可以看到5个文件（点、线、面存储）

数据读写

osm文件实质就是多个矢量图层组成的文件，读取之后按照shapefile图层读写方式处理即可。数据量少可以直接使用相关插件完成转换及导出。量大则可以考虑开源的命令行工具或代码完成处理

已知的可用工具：QGIS可以直接打开导出shape，ArcMap插件?

由于是矢量图层，直接使用GDAL中的矢量处理模块ogr完成处理

为了便于理解，读与写分开展示，实际上可以直接读取并完成写出。几何类型WKBGeometryType对应的代码

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
>>> from osgeo import ogr, osr
>>> inDS = ogr.Open('./moldova-latest.osm.pbf')
>>> inDS.GetLayerCount()
5
>>> layer = inDS.GetLayer(1)
>>> layer.GetGeomType()
2 
# 输出可以看到5种类型分别为1 2 5 6 7
# ogr.wkbPoint = 1
# ogr.wkbLineString = 2
# ogr.wkbMultiLineString = 5
# ogr.wkbMultiPolygon = 6
# ogr.wkbGeometryCollection = 7
>>> lydefn = layer.GetLayerDefn() 
>>> lydefn.GetFieldCount()
9
# 9列

读shape的字段，以线图层为例

1
2
3
4
5
6
7
8
9
10
11
12
13
#获取图层 定义及字段
lydefn = layer.GetLayerDefn() 
fieldlist = []
for i in range(lydefn.GetFieldCount()):
    fddefn = lydefn.GetFieldDefn(i)
    fddict = {
        'name':fddefn.GetName(),
        'type':fddefn.GetType(),
        'width':fddefn.GetWidth(),
        'decimal':fddefn.GetPrecision()}
    fieldlist.append(fddict)
fieldlist
# ['osm_id', 'name', 'highway', 'waterway', 'aerialway', 'barrier', 'man_made', 'z_order', 'other_tags']

读shape的要素

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
# 读取要素
geomlist = []
featurelist = []
feature = layer.GetNextFeature() #获得第一个要素
while feature is not None:
    geom = feature.GetGeometryRef()
    geomlist.append(geom.ExportToWkt())
    feature2 = {}
    for fd in fieldlist:
        feature2[fd['name']] = feature.GetField(fd['name'])
    featurelist.append(feature2)
    feature = layer.GetNextFeature()
ds=None # 释放数据源
# 预览数据
len(featurelist)
# 39008
featurelist[0]
"""
{'osm_id': '4529536',
 'name': None,
 'highway': None,
 'waterway': None,
 'aerialway': None,
 'barrier': None,
 'man_made': None,
 'z_order': 0,
 'other_tags': None}
"""

写shape,根据数据类型，就可以找到想要提取的文件并写出，可以通过修改fieldlist清理输出的shp文件的字段：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
# 创建输出文件
outShapefile = "./data_line2.shp"
outDriver = ogr.GetDriverByName("ESRI Shapefile")
outDS = outDriver.CreateDataSource(outShapefile)
srs = osr.SpatialReference()
srs.ImportFromEPSG(4326)
out_layer = outDS.CreateLayer(outShapefile, srs, geom_type=ogr.wkbLineString, options = ['ENCODING=UTF-8'])
# 写入字段
for fd in fieldlist:
    field = ogr.FieldDefn(fd['name'],fd['type'])
    if 'width' in fd:
        field.SetWidth(fd['width'])
    if 'decimal' in fd:
        field.SetPrecision(fd['decimal'])
    out_layer.CreateField(field)
# 写入要素
for i in range(len(featurelist)): 
    geom = ogr.CreateGeometryFromWkt(geomlist[i])
    # 创建空行
    newfeature = ogr.Feature(out_layer.GetLayerDefn())
    newfeature.SetGeometry(geom)
    # 写入字段及值
    for fd in fieldlist: 
        fieldname = fd['name']
        infeature = featurelist[i]
        newfeature.SetField(fieldname, infeature[fieldname]) 
    out_layer.CreateFeature(newfeature)
# 释放数据源并完成写出
inDS=None 
outDS = None

以上代码可以自行简化合并为一个部分，且可以根据需求进一步对shape进行分类导出操作。由于暂无需求暂不展开。