公告

[公告]
2014/01/17
由於已經是faculty的關係,不太有足夠時間寫部落格。因此更新的速度會相當緩慢。再加上近幾年來SAS GLOBAL FORUM沒有出現讓我覺得驚艷的技術文件,所以能分享的文章相對也減少許多。若有人推薦值得分享的SAS技術文件,請利用『問題討論區』告知。

2013/07/19
臉書留言板的功能因為有不明原因故障,因此特此移除。而intensedebate的留言板因管理不易,也一併移除。目前已經開啟內建的 G+ 留言系統,所以請有需要留言的朋友,可直接至『問題討論區』裡面留言。


2011年8月16日 星期二

Wow! You Did That Map with SAS/GRAPH®?

原文載點:http://www.nesug.info/Proceedings/nesug07/np/np01.pdf

身為一個專門從事空間時間序列分析的學者我來講,每天跟不同的地圖奮鬥是很稀鬆平常的是。目前市面上最熱門的地圖分析軟體莫過於ArcGIS,但這套軟體僅提供有限的分析方法,而其主要的功能僅僅就是繪製地圖而已。如果你用了別的軟體進行分析,就必須要將結果另存新檔,然後把結果跟地圖檔一起輸入到ArcGIS裡面才能開始繪圖。想當然爾,這裡面也沒有太多強大的資料整理功能,所以一切的一切都還是要經過另外的軟體將資料結果整理好才行。如果分析和繪製地圖的動作可以一次在SAS裡面完成的話將會省下許多時間和精力。NESUG 2007年有一篇教學文件,提供一些在SAS繪製地圖的資訊。

PROC MAPIMPORT

這是SAS專門用來輸入地圖資料的語法,不過必須先知道去哪裡可以找到地圖資料。以下是作者列出的幾個不錯的網頁資源:

(1) ESRI (a GIS software company)
http://arcdata.esri.com/data/tiger2000/tiger_download.cfm

(2) US Census Bureau
http://www.census.gov/geo/www/cob/index.html

(3) CDC (Centers for Disease Control)
http://www.cdc.gov/epiinfo/shape.htm

(4) CIP (International Potato Center, part of the Consultative Group on International Agricultural Research)
http://research.cip.cgiar.org/gis/index.php

其中,ESRI 就是 ArcGIS 的母公司,而美國地圖大都可以在 US Census Bureau 的網站上免費下載。CDC 的網站則是可以下載一個名叫Epi Info的軟體,裡面除了有地圖外還有其他的流病資料。CIP 的網站則是除了地圖以外還提供氣候資料。個人的經驗是 US Census Bureau 網站就足夠提供所有需要使用美國資料來進行研究的地圖。如果想進行台灣的研究,則可以到交通部的運研所下載:

http://www.iot.gov.tw/ct.asp?xItem=154948&ctNode=1091

如果要別的國家的地圖資料,有兩個網站可以下載:

http://www.diva-gis.org/gData
http://www.vdstech.com/map_data.htm

接下來的範例是用北卡(North Carolina)的資料畫地圖。先去 US Census Bureau 找北卡的 shape 檔:

North Carolina - zt37_d00_shp.zip (1,672,705 bytes)

解壓縮後存在自定目錄,然後執行下列程式:
/*Version 9.1 */
proc mapimport 
       datafile="path/zt37_d00.shp" 
       out=mymap91;
run;
/* Version 9.2 */
proc mapimport;
       datafile="path/zt37_d00.shp" 
       out=mymap92;
       id zcta;
run;
記得把上面的 path 改為自己的路徑。如此一來 SAS 就可以把 shapefile 裡面的資料都輸入進去。

MAPSONLINE 

這是 SAS 官網提供的線上工具,裡面也有地圖資料和一些其他資訊。網址如下:

http://support.sas.com/rnd/datavisualization/mapsonline/html/home.html

重點是如果你不小心把 sashelp 資料庫裡面一個叫做 ZIPCODE 的檔案刪掉的話,可以去上面這個網址重新下載,因為這個資料裡面包含美國所有的地理資訊,要畫美國地圖一定得連結到這個檔案。裡面還有一個很威的資料名為 WORLDCTS,包含全世界各大城市的基本資料,可以到這個連結直接下載:

http://support.sas.com/rnd/datavisualization/mapsonline/html/sampledata.html

CHOOSING MAP COLORS

在繪製地圖時,定義顏色是相當重要的一件事情。如何讓地圖上的顏色不會眼花撩亂並且能夠使讀者一目了然,這已經無關統計分析能力,而是一門藝術了。SAS 裡面對於顏色的指定可以使用顏色本來的英文名稱、RGB值、HLS值以及HEX值。所有的定義都可以去 SAS TS-688 的網站去看:

http://support.sas.com/techsup/technote/ts688/ts688.html

下面用一個簡單範例介紹:
%let color1=lightgreen ; 
%let color2=cornsilk ;
proc template;
   define style styles.grade;
   parent=styles.listing;
   style twocolorramp / startcolor=&color1 endcolor=&color2;
   end; 
run;
DISTANCE CALCULATIONS


計算點和點之間的距離是空間分析裡面的基本功,當知道兩個地點的經緯度(不知道的可以去用google map查詢)後,利用下面這個巨集程式便可輕鬆計算出來:
%macro geodist(lat1,long1,lat2,long2);
%let pi180=0.0174532925199433;
    7921.6623*arsin(sqrt((sin((&pi180*&lat2-&pi180*&lat1)/2))**2+
    cos(&pi180*&lat1)*cos(&pi180*&lat2)*(sin((&pi180*&long2-    
    &pi180*&long1)/2))**2));
%mend
如果是做美國的地理研究,則只需要知道每個 county 的 ZIP code 即可。在SAS V9.2版裡面有一個 ZIPCITYDIST() 函式就是專門用來計算兩個 ZIP code 之間的距離。範例如下:
%let zip1=27513;
%let zip2=21202;
/* latitude and longitude coordinates */
proc sql;
create table zip1 as select * from sashelp.zipcode where zip=&zip1;
create table zip2 as select * from sashelp.zipcode where zip=&zip2;
select x into :long1 from zip1;
select y into :lat1 from zip1;
select x into :long2 from zip2;
select y into :lat2 from zip2;
quit; run;
%let city1=%sysfunc(zipcity(&zip1));
%let city2=%sysfunc(zipcity(&zip2));
/* two new 9.2 functions */
%let zipdist=%sysfunc(zipcitydistance(&zip1,&zip2));
%let geodist=%sysfunc(geodist(&lat1, &long1, &lat2, &long2, 'M'),comma10.5);

GRAPHIC OVERLAYS


上面這張圖是利用 PROC GCHART 畫好柱狀圖,然後用 PROC GMAP 畫好地圖,之後把柱狀圖疊在地圖上面。首先把地圖畫好:
pattern v=me c=gray98;
proc gmap data=maps.us (obs=1) map=maps.us all anno=region_outline;
      id state;
      choro state / discrete nolegend coutline=gray98;
run;
quit;
再來把柱狀圖下面的文字設定好:
axis1 label=('NEW ENGLAND'        j=c '%CHANGE 48.7') ;
axis2 label=('MIDDLE ATLANTIC'    j=c '%CHANGE 21.1') ;
axis3 label=('SOUTH ATLANTIC'     j=c '%CHANGE 18.7') ;
axis4 label=('WEST NORTH CENTRAL' j=c '%CHANGE 63.3') ;
axis5 label=('EAST NORTH CENTRAL' j=c '%CHANGE 33.3') ;
axis6 label=('EAST SOUTH CENTRAL' j=c '%CHANGE 48.4') ;
axis7 label=('WEST SOUTH CENTRAL' j=c '%CHANGE 22.5') ;
axis8 label=('MOUNTAIN'           j=c '%CHANGE 50.9') ;
axis9 label=('PACIFIC'            j=c '%CHANGE 23.9') ;
axis10 order=0 to 20 by 5 label=none value=none style=0 major=none minor=none;
其中,axis1~axis9 自然就是九個柱狀圖 X 軸的文字和統計數據(這自己另外算好先),而 axis10 則表示用來控制 axis1~axis9 的 Y 軸共同設定。然後依序畫出每張柱狀圖,以右上角的 New England 為例:
proc gchart data=tpay;
goptions horigin=8.75 in vorigin=6.75 in;
where region eq 1;
vbar time / type=mean discrete sumvar=x mean raxis=axis10 maxis=axis1;
run;
quit;
資料集 tpay 裡面因為有所有 region 的資料,所以要用一個 where statement 限制只使用 New England (region = 1) 的數據。而 raxis 和 maxis 就是設定 Y 軸和 X 軸的選項,所以把 axis10 和 axis1 分別填入。當要畫別區的柱狀圖時,只要把 maxis 後面的設定換掉即可。至於 goptions 則是用來控制柱狀圖出現的位置。這就必須自己一直更換不同的數據來看位置對不對。你也可以用一個巨集程式來簡化步驟:
%macro bars(reg,h,v);
goptions horigin=&h in vorigin=&v in;
where region eq ®
vbar time / type=mean discrete sumvar=x raxis=axis10 maxis=axis&reg mean;
run;
%mend;
因此九張柱狀圖就可以用下面短短的程式碼一口氣完成:
proc gchart data=tpay;
%bars(1,8.75,6.45) %bars(2,7.75,4.90) %bars(3,7.50,3.00) 
%bars(4,4.00,5.00) %bars(5,6.00,4.25) %bars(6,6.10,2.50) 
%bars(7,4.00,1.85) %bars(8,1.95,4.20) %bars(9,0.00,2.00)
quit; 
最後一個步驟就是要把這些柱狀圖疊到地圖上。我們先用 ODS 開一個空白的 pdf 檔:
ods pdf file='z:\map_bars.pdf' notoc startpage=never;
後面的 notoc 和 startpage 主要是抑制 pdf 檔的書籤出現以及讓所有圖形都出現在同一頁。然後進行 pdf 檔裡面的環境設定:
goptions reset=all ftext='Helvetica/bo' htext=3.5 gunit=pct ctext=blue lfactor=3;

這些都是來設定字形、字色還有字體大小。這樣就大功告成了。整個流程可以用下列四步驟來簡單描述:

步驟一:用 GOPTIONS 進行圖形整體環境設定。
步驟二:用 ODS 開啟一個空白的 pdf 檔。
步驟三:用 PROC GMAP 畫地圖。
步驟四:用 PROC GCHART 畫柱狀圖

最後用 ODS CLOSE 把整個圖形打包起來。

REFERENCES & RECOMMENDED READING

以下是一些不錯的線上參考資料:

http://support.sas.com/documentation/onlinedoc/index.html
http://support.sas.com/rnd/datavisualization/mapsonline/html/home.html
http://support.sas.com/rnd/papers
http://support.sas.com/samples
http://ftp.sas.com/techsup/download/sample/graph/other-color.html
http://arcdata.esri.com/data/tiger2000/tiger_download.cfm
http://www.census.gov/geo/www/cob/index.html
http://www.cdc.gov/epiinfo/shape.htm
http://research.cip.cgiar.org/gis/index.php
http://support.sas.com/techsup/technote/ts688/ts688.html
http://www.colorbrewer.org
http://robslink.com/SAS/Home.htm

CONTACT INFORMATION
Your comments and questions are valued and encouraged.  Contact the authors at:
Robert Allison
robert.allison@sas.com
Louise Hadden
louise_hadden@abtassoc.com
Mike Zdeb
msz03@albany.edu
CODE { display: block; /* fixes a strange ie margin bug */ font-family: Courier New; font-size: 8pt; overflow:auto; background: #f0f0f0 url(http://klcintw.images.googlepages.com/Code_BG.gif) left top repeat-y; border: 1px solid #ccc; padding: 10px 10px 10px 21px; max-height:200px; height:200px; // for IE6 line-height: 1.2em; }