如何用SAS绘制现代化数字地图

文 | 巫银良

移动互联网应用和大规模社交网络催生了海量的数据分析需求，时空数据作为记录用户和设备在现实世界分布和活跃程度的基础数据，一直为各大互联网电子商务平台和商家所关注。地理空间数据结合其他业务数据如何被分析利用，以及如何在分析中可视化呈现一直是现代化分析平台的一个重要方向。

一方面各种地图服务越来越多地集成到应用中，成为应用增强交互的组成部分(比如“附近的服务/人”，甚至连支付宝红包都需要呈现各种方位关系，来增强乐趣)，另一方面在分析行业，如何能更高效方便地绘制各种地图成为一种基本需求。

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并且定期更新。

图1

据调查，MapGfK基础库包括2个世界级(其中一个world_cities为世界城市)，22个洲级，175个国家级6个美国州县的地理数据与对应属性数据。虽然看起来很全，但也并非十全十美，比如我在操作中发现台湾岛没有包括在中国数据内，2010年的北京核心四区合并为两区等也没有被更新数据。

下面，我们举个最简单的例子，来说明如何用SAS绘制现代化数字地图：

1、如何绘制世界地图

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

id id;

choro id / nolegend;

run;quit;

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

图2

2、如何绘制亚洲和中国地图

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

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

id id;

choro id / nolegend;

run;quit;

输出如图3结果：

图3

细心的你会发现，亚洲地图确实按照各个国家进行了准确的绘制，但中国地图则看起来黑压压的一片(这其实反映的是俺们大中华确实是地大物博啊~~)，九段线确实包括在内，但地图中竟然没有宝岛台湾。

3、如何将丢失的台湾部分合并到大中华地图

不要慌，下面我们介绍如何将缺失的台湾部分和中国地图合并为大中华地图。

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;

执行上面的代码，输出如图4结果。注意一下：台湾竟然出现在屏幕正中央(四川盆地)位置。

图4

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

…

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

id id;

run;

proc gmap map=…

输出如图5结果：

图5

4、如何优化地图

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

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;

输出如图6结果：

图6

5、如何给地图标记省份

为了给各省标注上省名，我们可以利用MAPSGFK库中已有的地图属性数据来绘制标签。此时需要利用系统自带的宏 %annomac 和 %maplabel 来生成描述数据数据。

data GreatChina_data;

set mapsgfk.China_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;

输出如图7结果：

图7

我们可以看到，图7的地图显示的是英文名称，而且台湾也没有包括在内?原因是 mapsgfk.China_attr 没有包含台湾的属性数据，需要我们自己创建。另外，我们要使用中文去显示各省的标注怎么办?只需要使用 id1nameU 并且转码即可显示正确的中文。

data mytaiwan_attr;

set mapsgfk.taiwan_attr;

id2=id1; id2name=id1name;

id1=’CN-83′; id1name=”Taiwan Sheng”; id1nameu= put(‘台湾省’,$uesc200.);

isoname=’China’;

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;

输出如图8结果：

图8

6、如何利用谷歌地球数据绘制地图

虽然在 MAPSGFK 基础库中有很多基础地理数据，但在现实中依然不够用怎么办? 解决方案有两种：

直接利用实际测量的地理数据创建自定义地图;

利用谷歌地球导出地球上任何地区/建筑的 KML 数据，然后导入到 SAS 里创建地图。

例如图9就是利用谷歌地球数据创建的谷歌总部第40号楼的地理数据。(图10为 Google Earth里的样子)

%MAPIMPORTX(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

ID

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;

执行上述代码生成的结果如图9：

图9

图10 Google Earth里的样子

在互联网上，有时听见一些人抱怨SAS语言做出的图表不够美观，显得比较粗陋。其实造成这误解的根本是没有掌握SAS强大的特性控制功能和实现的灵活性。

为了展示SAS在绘制地图方面预留的灵活性和特性控制，下面我将展示若干纯粹利用SAS代码绘制的各种现代化的复杂地图。

图11 SAS绘制空白中国省图

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

SAS语言天生作为面向分析而设计的语言，它保留了非常多的扩展性;我甚至发现在 SAS 地图里可以绘制天气云图(如图13所示)。正所谓倚天不出，谁与争锋?在分析行业里只有掌握了如何使用SAS这把倚天剑，才能使数据分析结果的展示一切皆有可能!

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

SAS GMAP 提供2D(choropleth)和3D(block, prism, surface)地图的绘制和渲染，用来将分析变量和结果显示在地图上。

既往的研究表明，SAS用户可以桥接任何地图服务商的数据，包括 MAPBOX，MAPQUEST, HERE，GOOGLE，ARCGIS和 AutoNavi(高德)的地图和他们的各种变体：SATELLITE(卫星图)，STREETS(街道图)，TERRAIN(地形图)和TRAFFIC(交通图)等。

PROC GMAP的所有奥秘其实都藏在它的MAP和DATA参数里，至于如何实现，就需要在实际需求中与具体业务数据结合考虑。

图14 利用SAS代码绘制现代化地图