社区首页 >专栏 >玩转地球: 如何利用SAS绘制现代化地图（附代码）

玩转地球: 如何利用SAS绘制现代化地图（附代码）

大数据文摘

发布于 2018-05-25 08:45:48

3.9K0

文章被收录于专栏：大数据文摘大数据文摘

投稿作者|巫银良

大数据文摘欢迎各类优质稿件

请联系tougao@bigdatadigest.cn

移动互联网应用和大规模社交网络催生了海量的数据分析需求，时空数据作为记录用户和设备在现实世界分布和活跃程度的基础数据，一直为各大互联网电子商务平台和商家所关注。地理空间数据结合其他业务数据如何被分析利用，以及如何在分析中可视化呈现一直是现代化分析平台的一个重要方向。一方面各种地图服务越来越多地集成到应用中，成为应用增强交互的组成部分（比如“附近的服务/人”，甚至连支付包红包都需要呈现各种方位关系，来增强乐趣），另一方面在分析行业，如何能够高效方便地绘制各种地图成为一种基本需求。

SAS 语言中提供了能够绘制地图的能力。考虑到 SAS 并不是地图数据的生产者，SAS 只是利用数据。在早些年，尽管SAS提供的地图数据来源多种多样，但SAS花了大量的时间精力来保证用户地图数据的精确性。随着现代卫星和测绘技术的成熟和一些其他原因（比如不再从CIA 获得世界数据），SAS 不再维护既往的地图数据，而是和第三方厂商合作来提供能够定期更新的最新数据，这样就不必考虑不同地理坐标系统和地缘政治格局变化带来的基础地理数据更新。

在传统上，SAS 缺省提供 MAPSSAS 库和 PROC GMAP, PROC GPROJECT, PROC GREMOVE, GEONCODE 等若干过程步来支持地图绘制功能。利用这些基础数据和过程步，用户能很容易绘制世界地图，各大洲地图，各国家地区地图。从SAS 9.30M2 版本开始，SAS 和 GFK GeoMarketing 合作，提供MPASGFK基础库，它能为用户提供超过240个国家和地区精准的数字邮政代码和行政区划地理数据。Gfk GeoMarketing 的数字地图是世界范围内最全最完整的数字地图，坐标系统为 WGS84 并且定期更新。根据笔者的调查，MapGfK 基础库包括2个世界级（其中一个world_cities为世界城市），22个洲级，175个国家级6个美国州县的地理数据与对应属性数据。虽然看起来很全，但也并非十全十美，比如笔者发现有些版本China地图数据没有包括中国台湾岛的内容，也没有反映2010年的北京核心四区合并为两区等变化。

下面，我们举个最简单的例子，来说明如何在SAS 里绘制地图:

proc gmap map=mapsgfk.world data=mapsgfk.world;

id id;

choro id / nolegend;

run;quit;

运行上面几行代码，SAS 会在结果窗口中输出如下结果：

如果你需要绘制亚洲或者中国地图，则只需要将上面world 改为 Asia或China 即可。

proc gmap map=mapsgfk.asia data=mapsgfk.asia;

id id;

choro id / nolegend;

run;quit;

细心的观众会发现，亚洲地图确实按照各个国家进行了准确的绘制，但中国地图则看起来黑压压的一片（…这个，其实反映的是俺们大中华确实是地大物博啊），并没有什么实用价值；九段线是确实包括在内，但其中竟然没有中国宝岛台湾（不要慌，下面我们介绍如何将缺失的中国台湾部分和中国地图合并为大中华地图）。

data mytaiwan;

set mapsgfk.taiwan;

id2=id1; id1='CN-83';

run;

data GreatChina;

set mapsgfk.china mytaiwan; /*合并中国台湾到大中华*/

run;

proc gmap map= GreatChina data= GreatChina;

id id;

choro id / nolegend;

run;quit;

执行上面的代码，输出结果如下。中国台湾出现在地图正中央（四川盆地）位置。原因是各个分区地图有自己的投影基点，我们需要按照中国数据进行投影。

为了将中国台湾岛移到指定位置，需要在调用 PROC GMAP 前执行如下代码，对中国台湾岛的数据根据中国的投影进行变换：

...

proc GPROJECT data=GreatChina out=GreatChina LATLON PARMIN=mapsgfk.projparm PARMENTRY=china;

id id;

run;

proc gmap map=..

在实际制作地图时，并不需要这么多的细节数据。因此我们需要将不必要的地区和县的边界删除，然后再调用 PROC GMAI绘图。代码如下：

proc sort data=GreatChina out=tmpds;

by ID1;

run;

proc gremove data=tmpds out=tmpds;

by ID1;

id id;

run;

data GreatChina(drop=ID1);

set tmpds;

id=ID1;

run;

为了给各省标注上省名，我们可以利用 MAPSGFK库中已有的地图属性数据来绘制标签。此时需要利用系统自带的宏 %annomac 和 %maplabel 来生成描述数据数据。另外，需要对中国台湾岛的描述数据进行特殊处理，统一到大中华地图中来。代码如下：

data mytaiwan_attr;

set mapsgfk.taiwan_attr;

id2=id1; id2name=id1name;

id1='CN-83'; id1name="Taiwan Sheng"; isoname='China';/*增补*/

drop country ;

run;

data GreatChina_data;

set mapsgfk.China_attr mytaiwan_attr; /*合并中国台湾省的描述数据*/

keep id1 id1name;

rename id1=id id1name=idname;

run;

%annomac;

%maplabel (GreatChina, GreatChina_data, anno_label, idname, id, font=%str(SimSun), color=black, size=1.5, hsys=3);

proc gmap map=GreatChina data=GreatChina_data;

id id;

choro id / nolegend anno=anno_label;

run;quit;

上面的地图显示的是英文名称，而我们希望显示中文名称怎么办？很简单，我们只需要在代码中使用 id1nameU 列，并将字符进行转义即可显示正确：

data mytaiwan_attr;

set mapsgfk.taiwan_attr;

id2=id1; id2name=id1name;

id1='CN-83'; id1name="Taiwan Sheng"; isoname='China';

id1nameu= put('中国台湾省',$uesc200.);

drop country ;

run;

data GreatChina_data;

set mapsgfk.China_attr mytaiwan_attr; /*合并中国台湾省的描述数据*/

keep id1 id1nameU;

rename id1=id id1nameU=idname;

run;

%annomac;

%maplabel (GreatChina, GreatChina_data, anno_label, idname, id, font=%str(SimSun), color=black, size=1.5, hsys=3);

data anno_label; set anno_label; text=unicode(text);run;

proc gmap map= GreatChina data= GreatChina_data;

id id;

choro id / nolegend anno=anno_label;

run;quit;

虽然在 MAPSGFK 基础库中有很多基础地理数据，但在现实中依然不够用怎么办？解决方案有两种：第一种是直接利用实际测绘的地理数据创建自定义地图；第二种方法是利用谷歌地球导出地球上任何地区/建筑的 KML 数据，然后再导入到 SAS 系统里创建地图。比如下图就是用第二种方法创建的谷歌总部第40号楼的地理数据。（下图为 Google Earth里的样子）

%MAPIMPORT(DATAFILE="test.kml",out=%str(mymap), ID=201);

data mymap; set mymap;

x=long;y=lat;

run;

data mymap_data;

attrib ID length=$15 label='Districts code';

attrib IDNAME length=$55 label='Districts name';

infile datalines dsd;

input

IDNAME

;

datalines4;

201,Google Building 40

;;;;

data mymap_data; set mymap_data;

length my_html $100;

my_html='title='||quote(trim(left(idname)));

run;

goptions reset=all;

goptions hsize=1024pt vsize=768pt;

ods html;

%let mymap=mymap;

%let mymap_data=mymap_data;

%annomac;

%maplabel (&mymap, &mymap_data, anno_label, idname, id, font=%str(SimSun), color=WHITE, size=1.5, hsys=3);

proc gmap map=mymap data=mymap_data;

id id;

choro id / nolegend anno=anno_label;

run;quit;

ods html close;

以上代码生成结果如下，为Google总部40号楼的精确地理信息，可用于进一步分析处理。

在互联网上，有时听见一些人抱怨 SAS 语言做出的图表不够美观，显得比较粗陋。其实造成这误解的根本是没有掌握 SAS 强大的特性控制功能和实现的灵活性。为了展示 SAS 在绘制地图方面预留的灵活性和控制，下面将展示若干纯粹利用 SAS 代码绘制的各种现代化的复杂地图。SAS语言天生作为面向分析而设计的语言，它保留了非常多的扩展性；笔者甚至发现在 SAS 地图里可以绘制天气云图（见下图3）。正所谓倚天不出，谁与争锋？在分析行业里只有掌握了如何使用SAS这把倚天剑，才能使数据分析结果的展示一切皆有可能！

图1：SAS绘制空白中国省图

图2：SAS绘制的中国各省的卫星地图

图3：SAS 绘制的带有卫星云图的中国分省图

总结：

SAS GMAP 提供 2D (choropleth) 和 3D (block, prism, surface) 地图的绘制和渲染，用来将分析变量和结果显示在地图上。既往的研究表明，SAS 用户可以桥接任何地图服务商的数据，包括 MAPBOX, MAPQUEST, HERE, GOOGLE，ARCGIS和 AutoNavi（高德）的地图和他们的各种变体：卫星图（SATELLITE）, （街道图）STREETS, （地形图）TERRAIN和（交通图）TRAFFIC 等。 PROC GMAP 的所有奥秘其实都藏在它的 MAP和DATA 参数里，至于如何实现，就需要在实际需求中与具体业务数据结合考虑。