OpenStreetMap数据Qt5分析实战(基于2020数据)
openstreetmap是一种完全开放的地理信息系统,数据由个人、公司免费捐赠、维护。本文在2019年的基础上,利用新的数据样本,展示Qt5作为棒哒哒的C++重量级框架的强悍。OpenStreetMap作为一款开源GIS,尽管数据的准确性、权威性较差,但数据本身还是有一定的规模,PostgreSQL矢量数据可超过450GB,购得上“大数据”的意思。
提到大数据,笔者不认为大数据一定是数据量非常大,一定要上集群。只要是基于数据分析的思想,对一定规模的数据进行分析并得到结论,都可以叫“大数据”。本文,使用Qt5分析openstreetmap数据库(样本为2020-01导入全球数据),获得城市科技指数这个自定义指标。
openSteetMap详细知识见https://planet.openstreetmap.org/。
我导入了一个虚拟机镜像,可从我的raspberryPi静态小网站获得链接。数据库的搭建见本博客专栏。
1 openstreetmap体现城市科技发达程度
在各类官方、非官方的统计学报告中,通过各种角度对城市进行排名。 有根据GDP的,根据绿化的,等等,但依托的一般都是权威地理信息生产商提供的数据。openStreetMap并不完美,其主要数据都来自不受共同利益约束的个人、小公司,因而准确性差,且很多内容存在错误,被正式场合使用的场景很少(想想在正式场合,地图上出现一个“二婶子家鱼塘”的标记会不会很尬(-:)。但是,由于数据来源很多、利益相关性低,从大数据的角度,基于这个样本的模型客观性较好。我们如果仔细观察,会发现OpenStreetMap的地标能够很好的体现一个地方的人类活动,特别是科技发展。
- 地标数量越多,说明为其贡献数据的独立个体越多——尽管小明标记了“二婶子家”这样的不规范地名,但至少代表这个城市有个小明在使用OpenStreetMap。
- 这些个体是城市的原住民(包括外地学生),从城市走出去的人(游子)、热爱这座城市的客人(游客)。
- 由于编辑、维护这些数据需要一定的技术基础,比如要会上网、打字,还要有一些基础的地理知识,我们可以认为贡献者应该代表了同时具备初、高中以上文化,且对计算机技术比较熟练的科技爱好者。
综上,一个地方地标越多,说明与该地相关的具有基础科技知识的人越多——而知识人才,是未来一个城市发展的希望。
2 设计相对公平的统计原理
各个城市的大小不一样,形状也不同。如果直接用外接多边形(行政区)进行统计,显然,面积越大的城市越占优势。解决此弊端,我们用城市中心区域500平方公里圆形范围(12.616千米半径)内为统计区域。
- 大部分城市的市区应该都能填充12千米为半径的圆,即使这个城市的形状非常奇怪(如长条型)。
- OpenStreetMap的城市要素(如place=city)point基本全部安排在最繁华的地点,因此,可以为这个圆找到非常棒的中心。
- 两个城市市区中心、CBD之间距离一般大于20公里。
- 如果两个城市的最繁华地区低于20公里,则构成城市带,如珠三角、长三角的很多县级市,这样统计也没有问题。
基于上述分析,采用最简单的统计的方法为查找各个地级市最繁华的半径为12公里的圆型区域中,具名地标有多少个。
3 SQL语句设计
PostgreSQL支持地理计算,以下语句会直接搜索出距离中心点最近的12.616千米内的地标ID(500平方公里的圆)。PostGIS允许SQL语句中加入地理空间计算,对集合的子交并补支持的非常好。
select distinct osm_id from planet_osm_line where ST_Covers(ST_Transform(ST_SetSRID(ST_MakeBox2D(ST_Point(112.556-0.5, 37.8562-0.5),ST_Point(112.556+0.5, 37.8562+0.5)),4326),3857),way) and ST_Distance(ST_Transform(ST_SetSRID(ST_MakePoint(112.556, 37.8562), 4326),3857),way) < 12616
上述功能需要解释一下:
- 城市的中心位置为 经度112.556,纬度37.8562
- 查询条件分为上下两部分。上部分为ST_Covers开始的一段,主要目的是加快查询速率。下半部分为测距,约束只查询12千米(500平方千米)的圆形区域。
- ST_MakeBox2D 创建一个大范围的矩形区域,包裹住12千米的范围。ST_SetSRID告诉数据库坐标为经纬度。ST_Transform转换经纬度为OpenStreetMap的摩卡托投影系。ST_Covers 计算上述矩形是不是覆盖到了某个地标。使用这大幅度提高GIS索引的查询效率。
- ST_Distance计算距离。这个距离是摩卡托的距离单位,在中低纬度,可近似认为是米。
- 上述SQL语句将迅速搜索400GB量级的矢量数据库,在几秒内返回结果。
4 地理数据获得
可以通过CSDN资源
https://download.csdn.net/download/u012266955/10664500
获得CSV格式的城市经纬度。由于该城市经纬度是行政区域的中心,因而需要和OpenStreetMap的位置进行矫正,以便找到最繁华的市中心。
该数据类似:
我们根据城市的名字,加上行政中心经纬度,使用名称在数据库搜索,并以地标最多的位置作为中心。
select name, St_X(St_Transform(way,4326)),St_Y(St_Transform(way,4326)) ,
ST_Distance(ST_Transform(ST_SetSRID(ST_MakePoint(116.753, 38.2658), 4326),3857),
way) as dis
from planet_osm_point where
ST_Covers(ST_Transform(ST_SetSRID(ST_MakeBox2D(ST_Point(116.753-2, 38.2658-2),ST_Point(116.753+2, 38.2658+2)),4326),3857),
way)
and place in ('city','town','suburb')
and (name like '沧州%' or name like '沧州市%')
order by admin_level, dis asc
5 统计结果
统计结果前50如下:
顺序 | 城市 | 省份 | 命名地点个数 |
---|---|---|---|
1 | 广州 | 广东 | 21907 |
2 | 上海 | 上海 | 17013 |
3 | 北京 | 北京 | 16937 |
4 | 珠海 | 广东 | 15763 |
5 | 深圳 | 广东 | 14001 |
6 | 杭州 | 浙江 | 11977 |
7 | 南昌 | 江西 | 9497 |
8 | 济南 | 山东 | 8282 |
9 | 成都 | 四川 | 7689 |
10 | 长沙 | 湖南 | 7309 |
11 | 合肥 | 安徽 | 7088 |
12 | 南通 | 江苏 | 6730 |
13 | 南京 | 江苏 | 6660 |
14 | 苏州 | 江苏 | 6360 |
15 | 佛山 | 广东 | 6206 |
16 | 武汉 | 湖北 | 5883 |
17 | 临沂 | 山东 | 4726 |
18 | 常州 | 江苏 | 4588 |
19 | 青岛 | 山东 | 4557 |
20 | 天津 | 天津 | 4355 |
21 | 兰州 | 甘肃 | 4350 |
22 | 江门 | 广东 | 4310 |
23 | 无锡 | 江苏 | 3739 |
24 | 潍坊 | 山东 | 3590 |
25 | 福州 | 福建 | 3225 |
26 | 长春 | 吉林 | 3136 |
27 | 昆明 | 云南 | 2959 |
28 | 哈尔滨 | 黑龙江 | 2937 |
29 | 清远 | 广东 | 2917 |
30 | 重庆 | 重庆 | 2768 |
31 | 唐山 | 河北 | 2695 |
32 | 汕头 | 广东 | 2621 |
33 | 沈阳 | 辽宁 | 2570 |
34 | 郑州 | 河南 | 2467 |
35 | 厦门 | 福建 | 2454 |
36 | 九江 | 江西 | 2399 |
37 | 中山 | 广东 | 2308 |
38 | 桂林 | 广西 | 2205 |
39 | 乌鲁木齐 | 新疆 | 2151 |
40 | 洛阳 | 河南 | 2104 |
41 | 湘潭 | 湖南 | 2053 |
42 | 鸡西 | 黑龙江 | 2038 |
43 | 西安 | 陕西 | 2020 |
44 | 泰安 | 山东 | 1996 |
45 | 大连 | 辽宁 | 1994 |
46 | 太原 | 山西 | 1892 |
47 | 宁波 | 浙江 | 1871 |
48 | 石家庄 | 河北 | 1870 |
49 | 贵阳 | 贵州 | 1844 |
50 | 自贡 | 四川 | 1832 |
… |
如果稍加改变,只选择博物馆的数目、大学的数目、厕所的数目,则可以体现出不同的属性,从而绘制出蛛网图。
6 源代码
细节主要包括:
- 使用opensteetmap中的多组地名进行地标中心选取,每个城市都取最多地标的中心作为参照。
- 只统计具有name属性的地标。以防很多自动标注工具批量生成的内容干扰统计。
- 下载的csv文件要添加到Qt5资源里,或者指定文件夹路径。
- Linux下,最好另存作为UTF-8编码,避免乱码。
#include <QCoreApplication>
#include <QFile>
#include <QSqlDatabase>
#include <QSqlQuery>
#include <QSqlError>
#include <QTextStream>
#include <QVector>
#include <cassert>QVector<QVector<double> > correctPosition(const double lon, const double lat,const QString & name1, const QString & name2);
qint64 queryCount(const double lon, const double lat, const QString & filter);int main(int argc, char *argv[]) {QCoreApplication a(argc, argv);QTextStream st_out(stdout,QIODevice::WriteOnly);try {//打开数据库,准备选择QSqlDatabase dbq = QSqlDatabase::addDatabase("QPSQL","TEST");if (!dbq.isValid())throw dbq.lastError();dbq.setHostName("127.0.0.1");dbq.setPort(5433);dbq.setUserName("archosm");dbq.setPassword("archosm");dbq.setDatabaseName("gis");if (!dbq.open())throw dbq.lastError();//读取CSV数据,https://download.csdn.net/download/u012266955/10664500QFile f_city(":/ctz_city_type.csv");QFile f_city_out("ctz_city_type.csv");if (!f_city.open(QIODevice::ReadOnly))throw f_city.errorString();if (!f_city_out.open(QIODevice::WriteOnly))throw f_city_out.errorString();QTextStream csv_in(&f_city);QTextStream csv_out(&f_city_out);csv_out<<"city,provice,citycode,citytype,cityclass,citype,conspeople,""cityname,cityym,lng_gd,lat_gd,spare,corr_lon,corr_lat,osm_count_1000km^2,url\n";//skip first line(titles)csv_in.readLine();while (!csv_in.atEnd()){const QString line = csv_in.readLine();const QStringList lst_elems = line.split(",");if (lst_elems.size()<11)continue;const double lon = lst_elems.at(9).toDouble();const double lat = lst_elems.at(10).toDouble();if (lon < 40 || lon >140 || lat<5 || lat >65)throw QString("Bad Input lat lon.");const QVector<QVector<double> > corrected =correctPosition(lon,lat,lst_elems.first(),lst_elems.at(7));qint64 best_hit = 0;double best_lat = lat, best_lon = lon;foreach(const QVector<double> & cter, corrected){const qint64 count = queryCount(cter[0],cter[1]," name is not null ");if (count>best_hit){best_hit = count;best_lon = cter[0];best_lat = cter[1];}}foreach (const QString & v, lst_elems) {csv_out<<v<<",";}csv_out<<best_lon<<","<<best_lat<<","<<best_hit<<",";csv_out<<QString("\"https://www.openstreetmap.org/#map=10/%2/%1\"").arg(best_lon).arg(best_lat)<<"\n";st_out << lst_elems.at(0)<<":"<<best_hit << "\n";st_out.flush();}f_city.close();f_city_out.close();}catch (QString err) {st_out << err << "\n";st_out.flush();}catch (QSqlError err) {st_out << err.text() << "\n";st_out.flush();}a.processEvents();return 0;
}qint64 queryCount(const double lon, const double lat, const QString & filter)
{//2000 km^2圆形区域内。QTextStream st_out(stdout,QIODevice::WriteOnly);static const QString s("select count(osm_id) from (""select distinct osm_id from (""select distinct osm_id from planet_osm_line where ST_Covers(""ST_Transform(ST_SetSRID(ST_MakeBox2D(ST_Point(%1-0.5, %2-0.5),ST_Point(%1+0.3, %2+0.3)),4326),3857),way) ""and ST_Distance(ST_Transform(ST_SetSRID(ST_MakePoint(%1, %2), 4326),3857),way) < 12616 and (%3) "" union ""select distinct osm_id from planet_osm_polygon where ST_Covers(""ST_Transform(ST_SetSRID(ST_MakeBox2D(ST_Point(%1-0.5, %2-0.5),ST_Point(%1+0.3, %2+0.3)),4326),3857),way) ""and ST_Distance(ST_Transform(ST_SetSRID(ST_MakePoint(%1, %2), 4326),3857),way) < 12616 and (%3) "" union ""select distinct osm_id from planet_osm_point where ST_Covers(""ST_Transform(ST_SetSRID(ST_MakeBox2D(ST_Point(%1-0.5, %2-0.5),ST_Point(%1+0.3, %2+0.3)),4326),3857),way) ""and ST_Distance(ST_Transform(ST_SetSRID(ST_MakePoint(%1, %2), 4326),3857),way) < 12616 and (%3) "" union ""select distinct osm_id from planet_osm_roads where ST_Covers(""ST_Transform(ST_SetSRID(ST_MakeBox2D(ST_Point(%1-0.5, %2-0.5),ST_Point(%1+0.3, %2+0.3)),4326),3857),way) ""and ST_Distance(ST_Transform(ST_SetSRID(ST_MakePoint(%1, %2), 4326),3857),way) < 12616 and (%3) "") all_eles) joined_uniq;");QSqlQuery query(QSqlDatabase::database("TEST"));QString sql = s.arg(lon).arg(lat).arg(filter);//st_out<<sql;if (query.exec(sql)==false)throw query.lastError();qint64 count = 0;if (query.next())count = query.value(0).toLongLong();return count;
}
QVector<QVector<double> > correctPosition(const double lon, const double lat, const QString & name1, const QString & name2)
{QVector<QVector<double> > res;QVector<double> lonlat_raw{lon,lat};//res<<lonlat;static const QString s("select St_X(St_Transform(way,4326)),St_Y(St_Transform(way,4326)) ,""ST_Distance(ST_Transform(ST_SetSRID(ST_MakePoint(%1, %2), 4326),3857),way) as dis,name from planet_osm_point"" where ST_Covers(""ST_Transform(ST_SetSRID(ST_MakeBox2D(ST_Point(%1-2, %2-2),ST_Point(%1+2, %2+2)),4326),3857),way) "" %5 and (name ='%3' or name ='%4') "" order by admin_level, dis asc ");QTextStream st_out(stdout,QIODevice::WriteOnly);QTextStream st_in(stdin,QIODevice::ReadOnly);QSqlQuery query(QSqlDatabase::database("TEST"));QString sql = s.arg(lon).arg(lat).arg(name1).arg(name2).arg(" and place in ('city','town','suburb') ");//st_out<<"\n"<<sql;if (query.exec(sql)==false)throw query.lastError();while (query.next()){QVector<double> lonlat{0,0};lonlat[0] = query.value(0).toDouble();lonlat[1] = query.value(1).toDouble();res<<lonlat;}if (res.size()==0){sql = s.arg(lon).arg(lat).arg(name1).arg(name2).arg(" and place is not null ");if (query.exec(sql)==false)throw query.lastError();while (query.next()){QVector<double> lonlat{0,0};lonlat[0] = query.value(0).toDouble();lonlat[1] = query.value(1).toDouble();res<<lonlat;}}if (res.size()==0){sql = s.arg(lon).arg(lat).arg(name1).arg(name2).arg(" ");if (query.exec(sql)==false)throw query.lastError();while (query.next()){QVector<double> lonlat{0,0};lonlat[0] = query.value(0).toDouble();lonlat[1] = query.value(1).toDouble();res<<lonlat;}}if (res.size()==0){res<<lonlat_raw;st_out<<"\n"<<sql<<" \nNot corrented\n";}return res;
}
代码在Linux Qt5.11 中编译通过。
题外-Qt是最棒的C++框架
无论是字符串处理,还是数据库读写,Qt都能够以静态语言的效率,做到动态语言才有的易用性。Qt5采用元对象系统,很好的模拟出动态语言的诸多特性。其实现的代码,可在Arm处理器、PC以及大型的Unix系统中运行,其具备的可视化模块(Charts、Data Visualization、Qwt)可非常生动地展示数据的全貌。使用Qt5实现功能的代码长度不但显著优于C++的其他框架,也直逼Python等动态语言。
Qt是各类著名软件的框架:
- Adobe Photoshop Album
- Autodesk Maya
- Autodesk 3ds Max
- Google Earth
- WPS
- Oracle VirtualBox
- Wireshark
- Ocatve
- OpenShot Video Editor
- Xilinx ISE
- LaTeX Lyx
- VLC
- Tableau desktop
- OpenPilot
- RStudio
- Skypy
OpenStreetMap数据Qt5分析实战(基于2020数据)相关推荐
- 《Wireshark数据包分析实战(第2版)》目录—导读
版权声明 Wireshark数据包分析实战(第2版) Copyright © 2011 by Chris Sanders. Title of English-language original:Pra ...
- R语言时间序列(time series)分析实战:时序数据加载、绘制时间序列图
R语言时间序列(time series)分析实战:时序数据加载.绘制时间序列图 目录
- 数据中台建设方案-基于大数据平台(下)
数据中台建设方案 -基于大数据平台(下) 1数据中台建设方案 1.1 总体建设方案 1.2大数据集成平台 1.3大数据计算平台 1.3.1数据计算层建设
- 数据中台建设方案-基于大数据平台
数据中台建设方案 -基于大数据平台- 1 数据中台建设方案 1.1 总体建设方案 通过对客户大数据应用平台服务需求的理解,根据建设目标.设计原则的多方面考虑,建议采用星环科技Transwarp Dat ...
- 数据中台建设方案-基于大数据平台(上)
数据中台建设方案 -基于大数据平台(上) 1 数据中台建设方案 1.1 总体建设方案 通过对客户大数据应用平台服务需求的理解,根据建设目标.设计原则的多方面考虑,建议采用星环科技Transwarp D ...
- 【Python】电商用户行为数据可视化分析实战
本文中,云朵君将和大家一起从多个角度使用多个可视化技术,根据各种因素跟踪客户在电子商务网站的花费时间. 关于数据集 数据集来自kaggle -- Machine Hack. 先进电子商务的用户数量激增 ...
- 《Wireshark数据包分析实战》读书笔记
1.OSI参考模型中的特殊功能: 表示层(第六层):进行用来保护数据的多种加密和解密操作. 会话层(第五层):负责以全双工或者半双工的方式来创建会话和关闭连接. 传输层(第四层):提供面向连接和无连接 ...
- 读书笔记 1.数据包分析技术与网络基础 Wireshark数据包分析实战 第3版
1.数据包分析技术与网络基础 1.2.1 协议 发起连接 :是由客户端还是服务器发起连接?在真正通信之前必须要交换哪些信息? 协商连接参数 :通信需要进行协议加密吗?加密密钥如何在通信双方进行传输? ...
- 【金猿产品展】亚信科技“数据探索分析平台”——深挖数据价值,助客户高效管理和经营生产...
亚信科技产品 本产品由亚信科技投递并参与"数据猿年度金猿策划活动--2020大数据产业创新服务产品榜单及奖项"评选. 大数据产业创新服务媒体 --聚焦数据 · 改变商业 亚信科技数 ...
最新文章
- php五只猴子分椰子_tubes五只雪茄_phillies雪茄五只装
- 论文简述 | 融合关键点和标记的基于图优化的可视化SLAM
- js 正则判断字符串是否为字母或数字
- MasterPage技术
- 笔记-计算机网络基础-开放系统互连参考模型OSI
- 冒泡排序的多种写法、逻辑
- mysql修改忘记了root密码忘记了,mysql忘记root密码后,重新设置、修改root密码
- matlab linspace
- mysql中数据定义语言_SQL数据定义语言(DDL)
- 五分钟没有操作自动退出_消防设施操作员 精选练习题10.31
- 【ES】ES 好文档积累
- Linux地图投影Proj4应用,Proj.4简介与使用
- 面试准备每日系列:计算机底层之并发编程(三)JVM-垃圾回收
- 关于: 为什么要写注释----谈一下个人体会
- 2021年4大免费ER图工具
- 局域网、城域网、广域网、国际互联网(internet)
- 4G和4G LTE之间的区别是什么?
- ROOT大师PC版 v1.7.6.7190 绿色免费版
- 金仓数据库 KingbaseES SQL 语言参考手册 (10. 查询和子查询)
- python群控微信_带你用 Python 实现自动化群控(入门篇)