Python通过OSM获取全国屋顶数据（经纬度、轮廓、面积、地址等）

引文

本文实现了使用Python通过OSM获取全国屋顶数据（经纬度、轮廓、面积、地址等）的效果，详细过程和效果可参考以下内容。

1. 导入相关依赖库

首先通过pip install pyrosm安装pyrosm，如果是Windows系统在安装中出现了问题，可参考https://pyrosm.readthedocs.io/en/latest/installation.html 尝试解决。

import requests
import pandas as pd
import os
from pyrosm import OSM

2. 下载各省的pbf文件

由于全国的pbf文件过大，而且在后续处理过程中容易爆内存，因此选择分各个省份进行处理后再合并。

url = 'http://download.openstreetmap.fr/extracts/asia/china/'
osm_dic = {
    'anhui': {'name_zh': '安徽','url': url + 'anhui-latest.osm.pbf'},
    'beijing': {'name_zh': '北京','url': url + 'beijing-latest.osm.pbf'},
    'chongqing': {'name_zh': '重庆','url': url + 'chongqing-latest.osm.pbf'},
    'fujian': {'name_zh': '福建','url': url + 'fujian-latest.osm.pbf'},
    'gansu': {'name_zh': '甘肃','url': url + 'gansu-latest.osm.pbf'},
    'guangdong': {'name_zh': '广东','url': url + 'guangdong-latest.osm.pbf'},
    'guangxi': {'name_zh': '广西','url': url + 'guangxi-latest.osm.pbf'},
    'guizhou': {'name_zh': '贵州','url': url + 'guizhou-latest.osm.pbf'},
    'hainan': {'name_zh': '海南','url': url + 'hainan-latest.osm.pbf'},
    'hebei': {'name_zh': '河北','url': url + 'hebei-latest.osm.pbf'},
    'heilongjiang': {'name_zh': '黑龙江','url': url + 'heilongjiang-latest.osm.pbf'},
    'henan': {'name_zh': '河南','url': url + 'henan-latest.osm.pbf'},
    'hong_kong': {'name_zh': '香港','url': url + 'hong_kong-latest.osm.pbf'},
    'hubei': {'name_zh': '湖北','url': url + 'hubei-latest.osm.pbf'},
    'hunan': {'name_zh': '湖南','url': url + 'hunan-latest.osm.pbf'},
    'inner_mongolia': {'name_zh': '内蒙古','url': url + 'inner_mongolia-latest.osm.pbf'},
    'jiangsu': {'name_zh': '江苏','url': url + 'jiangsu-latest.osm.pbf'},
    'jiangxi': {'name_zh': '江西','url': url + 'jiangxi-latest.osm.pbf'},
    'jilin': {'name_zh': '吉林','url': url + 'jilin-latest.osm.pbf'},
    'liaoning': {'name_zh': '辽宁','url': url + 'liaoning-latest.osm.pbf'},
    'macau': {'name_zh': '澳门','url': url + 'macau-latest.osm.pbf'},
    'ningxia': {'name_zh': '宁夏','url': url + 'ningxia-latest.osm.pbf'},
    'qinghai': {'name_zh': '青海','url': url + 'qinghai-latest.osm.pbf'},
    'shaanxi': {'name_zh': '陕西','url': url + 'shaanxi-latest.osm.pbf'},
    'shandong': {'name_zh': '山东','url': url + 'shandong-latest.osm.pbf'},
    'shanghai': {'name_zh': '上海','url': url + 'shanghai-latest.osm.pbf'},
    'shanxi': {'name_zh': '山西','url': url + 'shanxi-latest.osm.pbf'},
    'sichuan': {'name_zh': '四川','url': url + 'sichuan-latest.osm.pbf'},
    'tianjin': {'name_zh': '天津','url': url + 'tianjin-latest.osm.pbf'},
    'tibet': {'name_zh': '西藏','url': url + 'tibet-latest.osm.pbf'},
    'xinjiang': {'name_zh': '新疆','url': url + 'xinjiang-latest.osm.pbf'},
    'yunnan': {'name_zh': '云南','url': url + 'yunnan-latest.osm.pbf'},
    'zhejiang': {'name_zh': '浙江','url': url + 'zhejiang-latest.osm.pbf'},
    'taiwan': {'name_zh': '台湾','url': 'https://download.geofabrik.de/asia/taiwan-latest.osm.pbf'}
}

tip：由于目标文件在外网，所以下载速度会比较慢，使用科学上网工具可显著提高文件的下载速度！

path = 'osm_data/'
if not os.path.exists(path):
    os.makedirs(path)

def get_all_pbf():
    def download_pbf(file_path, url):
        headers={'User-Agent':'Mozilla/5.0 (X11; Linux x86_64) AppleWebKit/537.36 (KHTML, like Gecko) Chrome/94.0.4606.81 Safari/537.36'}
        response = requests.get(url=url,headers=headers).content
        with open(file_path, "wb") as f:
            f.write(response)

    for location in osm_dic:
        file_path = path + location + '.osm.pbf'
        url = osm_dic[location]['url']
        try:
            download_pbf(file_path, url)
            print(location+'.osm.pbf下载成功！文件目录为：' + file_path)
        except:
            print(location+'.osm.pbf下载失败！')

get_all_pbf()

3. 读取osm.pbf的建筑物数据，保存为csv文件

def pbf2csv(location):
    try:
        osm = OSM(path + location + '.osm.pbf')
        buildings = osm.get_buildings()
        buildings = buildings[buildings['geometry'].apply(lambda x: 'POLYGON' in str(x))][['addr:city', 'addr:postcode', 'addr:street', 'name', 
        'building','amenity', 'building:levels', 'height', 'office', 'shop', 'id', 'timestamp', 'version', 'geometry', 'tags','osm_type']]
        buildings.to_csv(path + location + '.csv', index=0, encoding='utf_8_sig')
        print(location+'.csv 提取成功！')
    except:
        print(location+'.csv 提取失败！！！')

for location in osm_dic:
    pbf2csv(location)

4. 简单展示数据结果

import pandas as pd
import geopandas as gpd
from shapely import wkt

df = pd.read_csv(path  + 'macau.csv')
df['geometry'] = df['geometry'].apply(wkt.loads)
gdf=gpd.GeoDataFrame(df,geometry=df.geometry)
gdf.plot()

5. 全国屋顶数据

由于篇幅所限，后续还有通过高德API解析地址、wgs_84坐标系和gcj_02的坐标系转换等步骤没有详细记录。
且由于全国屋顶数据量较大，需要最终结果文件的朋友可以私信我获取。