openstreetmap是一种完全开放的地理信息系统，数据由个人、公司免费捐赠、维护。在这个博客的前文中，我们大多围绕搭建地图环境展开讨论。实际上，它更具价值的是数据本身。今天，我们来看使用Qt5分析openstreetmap数据库（样本为2019-01导入全球数据），获得城市科技指数这个自定义指标。

openSteetMap详细知识、数据、虚拟机见:
http://www.goldenhawking.org:8088/
https://www.openstreetmap.org
https://planet.openstreetmap.org/

数据库的搭建见本博客其他博文。专栏：
https://blog.csdn.net/goldenhawking/column/info/22354

1 openstreetmap体现城市科技发达程度

在各类官方、非官方的统计学报告中，通过各种角度对城市进行排名。有根据GDP的，根据绿化的，等等。openStreetMap作为在国内并不很流行的GIS，因为主要数据都来自不受共同利益约束的个人、小公司，因而更能体现一个地方的科技发达程度。

地标数量越多，说明为其贡献数据的独立个体越多。
这些个体可能是城市的原住民，或者和城市有关的人（游客、学生、社团）。
由于编辑、维护这些数据需要一定的技术基础，比如要会上网、打字，还要有一些基础的地理知识，我们可以认为贡献者应该代表了同时具备初中以上文化，且对计算机技术比较熟练的人。
综上，一个地方地标越多，说明与该地相关的具有基础科技知识的人越多——而人，是未来一个城市发展的希望。

2 设计相对公平的统计原理

各个城市的大小不一样，形状也不同。如果直接用外接多边形（行政区）进行统计，显然，面积越大的城市越占优势。我们用城市中心区域500平方公里圆形范围（12.616千米半径）内的地标个数为统计方法，主要理由：

大部分城市的市区应该都能填充12千米为半径的圆，即使这个城市的形状非常奇怪（如长条型）
OpenStreetMap的城市要素（如place=city）point基本全部安排在最繁华的地点，因此，可以为这个圆找到非常棒的中心。
两个城市市区中心、CBD之间距离一般大于20公里。

因此，统计的方法为查找各个地级市最繁华的半径为12公里的圆型区域中，地标有多少个。

3 SQL语句设计

PostgreSQL支持地理计算，以下语句会直接搜索出距离中心点最近的5.642千米内的地标ID(100平方公里的圆）。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.3, 37.8562+0.3)),4326),3857),way) and ST_Distance(ST_Transform(ST_SetSRID(ST_MakePoint(112.556, 37.8562), 4326),3857),way) < 12616

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 统计结果

统计结果如下:

顺序	城市	省份	命名地点个数
1	广州	广东	19237
2	北京	北京	15733
3	上海	上海	14878
4	珠海	广东	12470
5	深圳	广东	12134
6	杭州	浙江	11178
7	南昌	江西	9208
8	济南	山东	6593
9	合肥	安徽	6463
10	南通	江苏	6457
11	成都	四川	6446
12	南京	江苏	6335
13	苏州	江苏	5849
14	武汉	湖北	5667
15	青岛	山东	4547
16	临沂	山东	4511
17	佛山	广东	4389
18	西安	陕西	4197
19	天津	天津	4056
20	长沙	湖南	3939
21	兰州	甘肃	3697
22	潍坊	山东	3222
23	江门	广东	3162
24	福州	福建	3028
25	长春	吉林	2707
26	无锡	江苏	2683
27	唐山	河北	2597
28	汕头	广东	2593
29	清远	广东	2572
30	昆明	云南	2455
31	常州	江苏	2439
32	厦门	福建	2154
33	中山	广东	2071
34	湘潭	湖南	2004
35	鸡西	黑龙江	1975
36	洛阳	河南	1954
37	大连	辽宁	1832
38	张家口	河北	1785
39	平顶山	河南	1725
40	宁波	浙江	1706
41	泰安	山东	1700
42	贵阳	贵州	1642
43	石家庄	河北	1597
44	哈尔滨	黑龙江	1536
45	九江	江西	1500
46	遵义	贵州	1405
47	宁德	福建	1391
48	蚌埠	安徽	1376
49	镇江	江苏	1374
50	桂林	广西	1374
…

6 源代码

细节主要包括：

使用opensteetmap中的多组地名进行地标中心选取，每个城市都取最多地标的中心作为参照。
只统计具有name属性的地标。以防很多自动标注工具批量生成的内容干扰统计。

#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;
}

openStreetMap数据分析举例-使用Qt统计城市科技指数排名相关推荐

智慧城市发展指数中国第一，深圳数字化转型全景展示
借助新兴的信息技术能够随时随地感知.捕获.传递和处理信息,能够实现对城市的精细化.智能化管理,完善公共服务设施和城市功能,创造安全.便捷.高效和清洁环保的环境,让人们能够在智慧城市中更好生活. 目前, ...
Python数据分析pandas之分组统计透视表
Python数据分析pandas之分组统计透视表数据聚合统计 Padans里的聚合统计即是应用分组的方法对数据框进行聚合统计,常见的有min(最小).max(最大).avg(平均值).sum(求和) ...
中国硬科技城市发展指数正式发布，西安跻身前十
硬科技属于由科技创新构成的物理世界,是需要长期研发投入.持续积累才能形成的原创技术,其具有极高技术门槛和技术壁垒,难以被复制和模仿. 11月7日上午,2017全球硬科技创新大会在西安开幕,吸引了全球众 ...
2008全球城市竞争力最新排名出炉
东莞入围近五年发展最快城市 "东莞奇迹"再次吸引了世界的眼光.第五届城市竞争力国际论坛昨日(27日)在扬州发布了<2007-2008全球城市竞争力报告>(以下简称< ...
数字城市网络安全指数白皮书附下载
随着"智慧城市""新基建"等战略的持续推进,5G.工业互联网.边缘计算等技术的快速发展,物联网与关键信息基础设施深度融合,在提高行业运行效率和便捷性的同时,也面 ...
python数据分析怎么样_Python数据分析必知必会：TGI指数
这是Python数据分析实战的第一个案例,详细解读TGI指数,并用Python代码实现基础的TGI偏好分析. 经常有一些专业的数据分析报告,会提到TGI指数,例如"基于某某TGI指数,我们发 ...
数字政府发展指数排名出炉！上海、浙江、北京位列前三，你的城市排第几？（附报告全文下载）...
本文约3600字,建议阅读6分钟. <2020数字政府发展指数报告>原创性地设计了中国数字政府发展指数评估指标体系. 关注微信公众号"数据派THU",对话框发送关键词& ...
2004中国城市综合竞争力排名
作者:天马行空(anlaiwuyou)个人主页:天马行空(anlaiwuyou) 出处:茶亦醉人何必酒书能香我无须花-敏思博客个人主页:茶亦醉人何必酒书能香我无须花-敏思博客发表于:2004年 ...
全国所有城市人均GDP排名（包含县级市 611 ）
全国所有城市人均GDP排名(包含县级市 611 zz) 排名人均国内生产总值 (元) 1 深圳(粤) 136071.3 2 大庆(黑) 89962.56 3 珠海(粤) 66550.61 4 ...
2014中国城市“鬼城”指数排行榜发布
2014中国城市"鬼城"指数排行榜发布行业动态投资时报宿小庆2014-10-12 07:37 我要分享 452 中国找"鬼"记中国或现50个"鬼城 ...

openStreetMap数据分析举例-使用Qt统计城市科技指数排名