身為一個專門從事空間時間序列分析的學者我來講，每天跟不同的地圖奮鬥是很稀鬆平常的是。目前市面上最熱門的地圖分析軟體莫過於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;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;計算點和點之間的距離是空間分析裡面的基本功，當知道兩個地點的經緯度（不知道的可以去用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%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;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;%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® 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 file='z:\map_bars.pdf' notoc startpage=never;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
沒有留言:
張貼留言
要問問題的人請在文章下方的intensedebate欄位留言，請勿使用blogger預設的意見表單。今後用blogger意見表單留言的人我就不回應了。