玩轉地球: 如何利用SAS繪製現代化地圖(附代碼)

投稿作者|巫銀良

大數據文摘歡迎各類優質稿件

請聯繫[e-mail]

移動網際網路應用和大規模社交網絡催生了海量的數據分析需求，時空數據作為記錄用戶和設備在現實世界分布和活躍程度的基礎數據，一直為各大網際網路電子商務平臺和商家所關注。地理空間數據結合其他業務數據如何被分析利用，以及如何在分析中可視化呈現一直是現代化分析平臺的一個重要方向。一方面各種地圖服務越來越多地集成到應用中，成為應用增強交互的組成部分（比如「附近的服務/人」，甚至連支付包紅包都需要呈現各種方位關係，來增強樂趣），另一方面在分析行業，如何能夠高效方便地繪製各種地圖成為一種基本需求。

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

      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;

以上代碼生成結果如下，為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 參數裡，至於如何實現，就需要在實際需求中與具體業務數據結合考慮。

賽仕軟體研究開發（北京）有限公司

商業智能和可視化分析產品部 技術總監

如需轉載，請在開篇顯著位置註明作者和出處（轉自：大數據文摘 |bigdatadigest），並在文章結尾放置大數據文摘醒目二維碼。無原創標識文章請按照轉載要求編輯，可直接轉載，轉載後請將轉載連結發送給我們；有原創標識文章，請發送【文章名稱-待授權公眾號名稱及ID】給我們申請白名單授權。未經許可的轉載以及改編者，我們將依法追究其法律責任。聯繫郵箱：[e-mail]。
點擊圖片閱讀文章
Kaggle | 使用Python和R繪製數據地圖的十七個經典案例(附資源)