标题: Wow! You Did That Map with SAS/GRAPH®? [打印本页] 作者: shiyiming 时间: 2011-8-17 06:27 标题: Wow! You Did That Map with SAS/GRAPH®? From LCChien's blog on blogspot
原文載點：<a href="http://www.nesug.info/Proceedings/nesug07/np/np01.pdf"><!-- m --><a class="postlink" href="http://www.nesug.info/Proceedings/nesug07/np/np01.pdf">http://www.nesug.info/Proceedings/nesug07/np/np01.pdf</a><!-- m --></a><br /><br />身為一個專門從事空間時間序列分析的學者我來講，每天跟不同的地圖奮鬥是很稀鬆平常的是。目前市面上最熱門的地圖分析軟體莫過於ArcGIS，但這套軟體僅提供有限的分析方法，而其主要的功能僅僅就是繪製地圖而已。如果你用了別的軟體進行分析，就必須要將結果另存新檔，然後把結果跟地圖檔一起輸入到ArcGIS裡面才能開始繪圖。想當然爾，這裡面也沒有太多強大的資料整理功能，所以一切的一切都還是要經過另外的軟體將資料結果整理好才行。如果分析和繪製地圖的動作可以一次在SAS裡面完成的話將會省下許多時間和精力。NESUG 2007年有一篇教學文件，提供一些在SAS繪製地圖的資訊。<br /><br /><a name='more'></a><b><span class="Apple-style-span" style="font-size: large;">PROC MAPIMPORT</span></b><br /><br />這是SAS專門用來輸入地圖資料的語法，不過必須先知道去哪裡可以找到地圖資料。以下是作者列出的幾個不錯的網頁資源：<br /><br />(1) ESRI (a GIS software company)<br /><a href="http://arcdata.esri.com/data/tiger2000/tiger_download.cfm"><!-- m --><a class="postlink" href="http://arcdata.esri.com/data/tiger2000/tiger_download.cfm">http://arcdata.esri.com/data/tiger2000/ ... wnload.cfm</a><!-- m --></a><br /><br />(2) US Census Bureau<br /><a href="http://www.census.gov/geo/www/cob/index.html"><!-- m --><a class="postlink" href="http://www.census.gov/geo/www/cob/index.html">http://www.census.gov/geo/www/cob/index.html</a><!-- m --></a><br /><br />(3) CDC (Centers for Disease Control)<br /><a href="http://www.cdc.gov/epiinfo/shape.htm"><!-- m --><a class="postlink" href="http://www.cdc.gov/epiinfo/shape.htm">http://www.cdc.gov/epiinfo/shape.htm</a><!-- m --></a><br /><br />(4) CIP (International Potato Center, part of the Consultative Group on International Agricultural Research)<br /><a href="http://research.cip.cgiar.org/gis/index.php"><!-- m --><a class="postlink" href="http://research.cip.cgiar.org/gis/index.php">http://research.cip.cgiar.org/gis/index.php</a><!-- m --></a><br /><br />其中，ESRI 就是 ArcGIS 的母公司，而美國地圖大都可以在 US Census Bureau 的網站上免費下載。CDC 的網站則是可以下載一個名叫Epi Info的軟體，裡面除了有地圖外還有其他的流病資料。CIP 的網站則是除了地圖以外還提供氣候資料。個人的經驗是 US Census Bureau 網站就足夠提供所有需要使用美國資料來進行研究的地圖。如果想進行台灣的研究，則可以到交通部的運研所下載：<br /><br /><a href="http://www.iot.gov.tw/ct.asp?xItem=154948&ctNode=1091"><!-- m --><a class="postlink" href="http://www.iot.gov.tw/ct.asp?xItem=154948&ctNode=1091">http://www.iot.gov.tw/ct.asp?xItem=154948&ctNode=1091</a><!-- m --></a><br /><br />如果要別的國家的地圖資料，有兩個網站可以下載：<br /><br /><a href="http://www.diva-gis.org/gData"><!-- m --><a class="postlink" href="http://www.diva-gis.org/gData">http://www.diva-gis.org/gData</a><!-- m --></a><br /><a href="http://www.vdstech.com/map_data.htm"><!-- m --><a class="postlink" href="http://www.vdstech.com/map_data.htm">http://www.vdstech.com/map_data.htm</a><!-- m --></a><br /><br />接下來的範例是用北卡（North Carolina）的資料畫地圖。先去 US Census Bureau 找北卡的 shape 檔：<br /><br /><span class="Apple-style-span" style="-webkit-border-horizontal-spacing: 15px; -webkit-border-vertical-spacing: 15px; font-family: arial, helvetica; font-size: 13px;">North Carolina -</span><span class="Apple-style-span" style="-webkit-border-horizontal-spacing: 15px; -webkit-border-vertical-spacing: 15px; font-family: arial, helvetica; font-size: 13px;"> </span><span class="Apple-style-span" style="-webkit-border-horizontal-spacing: 15px; -webkit-border-vertical-spacing: 15px; font-family: arial, helvetica; font-size: 13px;"><a href="http://www.census.gov/geo/cob/bdy/zt/z500shp/zt37_d00_shp.zip">zt37_d00_shp.zip</a></span><span class="Apple-style-span" style="-webkit-border-horizontal-spacing: 15px; -webkit-border-vertical-spacing: 15px; font-family: arial, helvetica; font-size: 13px;"> </span><span class="Apple-style-span" style="-webkit-border-horizontal-spacing: 15px; -webkit-border-vertical-spacing: 15px; font-family: arial, helvetica; font-size: 13px;">(1,672,705 bytes)</span><br /><br />解壓縮後存在自定目錄，然後執行下列程式：<br /><pre><code>/*Version 9.1 */<br />proc mapimport <br /> datafile="<span class="Apple-style-span" style="color: red;">path</span>/zt37_d00.shp" <br /> out=mymap91;<br />run;<br />/* Version 9.2 */<br />proc mapimport;<br /> datafile="<span class="Apple-style-span" style="color: red;">path</span>/zt37_d00.shp" <br /> out=mymap92;<br /> id zcta;<br />run;</code></pre>記得把上面的 path 改為自己的路徑。如此一來 SAS 就可以把 shapefile 裡面的資料都輸入進去。<br /><br /><b><span class="Apple-style-span" style="font-size: large;">MAPSONLINE </span></b><br /><br />這是 SAS 官網提供的線上工具，裡面也有地圖資料和一些其他資訊。網址如下：<br /><br /><a href="http://support.sas.com/rnd/datavisualization/mapsonline/html/home.html"><!-- m --><a class="postlink" href="http://support.sas.com/rnd/datavisualization/mapsonline/html/home.html">http://support.sas.com/rnd/datavisualiz ... /home.html</a><!-- m --></a><br /><br />重點是如果你不小心把 sashelp 資料庫裡面一個叫做 ZIPCODE 的檔案刪掉的話，可以去上面這個網址重新下載，因為這個資料裡面包含美國所有的地理資訊，要畫美國地圖一定得連結到這個檔案。裡面還有一個很威的資料名為 WORLDCTS，包含全世界各大城市的基本資料，可以到這個連結直接下載：<br /><br /><a href="http://support.sas.com/rnd/datavisualization/mapsonline/html/sampledata.html"><!-- m --><a class="postlink" href="http://support.sas.com/rnd/datavisualization/mapsonline/html/sampledata.html">http://support.sas.com/rnd/datavisualiz ... edata.html</a><!-- m --></a><br /><br /><b><span class="Apple-style-span" style="font-size: large;">CHOOSING MAP COLORS</span></b><br /><br />在繪製地圖時，定義顏色是相當重要的一件事情。如何讓地圖上的顏色不會眼花撩亂並且能夠使讀者一目了然，這已經無關統計分析能力，而是一門藝術了。SAS 裡面對於顏色的指定可以使用顏色本來的英文名稱、RGB值、HLS值以及HEX值。所有的定義都可以去 SAS TS-688 的網站去看：<br /><br /><a href="http://support.sas.com/techsup/technote/ts688/ts688.html"><!-- m --><a class="postlink" href="http://support.sas.com/techsup/technote/ts688/ts688.html">http://support.sas.com/techsup/technote ... ts688.html</a><!-- m --></a><br /><br />下面用一個簡單範例介紹：<br /><pre><code>%let color1=lightgreen ; <br />%let color2=cornsilk ;<br />proc template;<br /> define style styles.grade;<br /> parent=styles.listing;<br /> style twocolorramp / startcolor=&color1 endcolor=&color2;<br /> end; <br />run;</code></pre><b><span class="Apple-style-span" style="font-size: large;">DISTANCE CALCULATIONS</span></b><br /><br /><br />計算點和點之間的距離是空間分析裡面的基本功，當知道兩個地點的經緯度（不知道的可以去用google map查詢）後，利用下面這個巨集程式便可輕鬆計算出來：<br /><pre><code>%macro geodist(lat1,long1,lat2,long2);<br />%let pi180=0.0174532925199433;<br /> 7921.6623*arsin(sqrt((sin((&pi180*&lat2-&pi180*&lat1)/2))**2+<br /> cos(&pi180*&lat1)*cos(&pi180*&lat2)*(sin((&pi180*&long2- <br /> &pi180*&long1)/2))**2));<br />%mend</code></pre>如果是做美國的地理研究，則只需要知道每個 county 的 ZIP code 即可。在SAS V9.2版裡面有一個 ZIPCITYDIST() 函式就是專門用來計算兩個 ZIP code 之間的距離。範例如下：<br /><pre><code>%let zip1=27513;<br />%let zip2=21202;<br />/* latitude and longitude coordinates */<br />proc sql;<br />create table zip1 as select * from sashelp.zipcode where zip=&zip1;<br />create table zip2 as select * from sashelp.zipcode where zip=&zip2;<br />select x into :long1 from zip1;<br />select y into :lat1 from zip1;<br />select x into :long2 from zip2;<br />select y into :lat2 from zip2;<br />quit; run;<br />%let city1=%sysfunc(zipcity(&zip1));<br />%let city2=%sysfunc(zipcity(&zip2));<br />/* two new 9.2 functions */<br />%let zipdist=%sysfunc(zipcitydistance(&zip1,&zip2));<br />%let geodist=%sysfunc(geodist(&lat1, &long1, &lat2, &long2, 'M'),comma10.5);</code></pre><br /><b><span class="Apple-style-span" style="font-size: large;">GRAPHIC OVERLAYS</span></b><br /><br /><div class="separator" style="clear: both; text-align: center;"><a href="http://2.bp.blogspot.com/-NpYERIPYae0/TkrZhYnSiuI/AAAAAAAAGR4/eMRd_vqXYXg/s1600/us_map_area.png" imageanchor="1" style="margin-left: 1em; margin-right: 1em;"><img border="0" height="425" src="http://2.bp.blogspot.com/-NpYERIPYae0/TkrZhYnSiuI/AAAAAAAAGR4/eMRd_vqXYXg/s640/us_map_area.png" width="640" /></a></div><br />上面這張圖是利用 PROC GCHART 畫好柱狀圖，然後用 PROC GMAP 畫好地圖，之後把柱狀圖疊在地圖上面。首先把地圖畫好：<br /><pre><code>pattern v=me c=gray98;<br />proc gmap data=maps.us (obs=1) map=maps.us all anno=region_outline;<br /> id state;<br /> choro state / discrete nolegend coutline=gray98;<br />run;<br />quit;</code></pre>再來把柱狀圖下面的文字設定好：<br /><pre><code>axis1 label=('NEW ENGLAND' j=c '%CHANGE 48.7') ;<br />axis2 label=('MIDDLE ATLANTIC' j=c '%CHANGE 21.1') ;<br />axis3 label=('SOUTH ATLANTIC' j=c '%CHANGE 18.7') ;<br />axis4 label=('WEST NORTH CENTRAL' j=c '%CHANGE 63.3') ;<br />axis5 label=('EAST NORTH CENTRAL' j=c '%CHANGE 33.3') ;<br />axis6 label=('EAST SOUTH CENTRAL' j=c '%CHANGE 48.4') ;<br />axis7 label=('WEST SOUTH CENTRAL' j=c '%CHANGE 22.5') ;<br />axis8 label=('MOUNTAIN' j=c '%CHANGE 50.9') ;<br />axis9 label=('PACIFIC' j=c '%CHANGE 23.9') ;<br />axis10 order=0 to 20 by 5 label=none value=none style=0 major=none minor=none;</code></pre>其中，axis1~axis9 自然就是九個柱狀圖 X 軸的文字和統計數據（這自己另外算好先），而 axis10 則表示用來控制 axis1~axis9 的 Y 軸共同設定。然後依序畫出每張柱狀圖，以右上角的 New England 為例：<br /><pre><code>proc gchart data=tpay;<br />goptions <span class="Apple-style-span" style="color: blue;">horigin=8.75 in vorigin=6.75 in</span>;<br /><span class="Apple-style-span" style="color: red;">where region eq 1;</span><br />vbar time / type=mean discrete sumvar=x mean<span class="Apple-style-span" style="color: red;"> raxis=axis10 maxis=axis1</span>;<br />run;<br />quit;</code></pre>資料集 tpay 裡面因為有所有 region 的資料，所以要用一個 where statement 限制只使用 New England (region = 1) 的數據。而 raxis 和 maxis 就是設定 Y 軸和 X 軸的選項，所以把 axis10 和 axis1 分別填入。當要畫別區的柱狀圖時，只要把 maxis 後面的設定換掉即可。至於 goptions 則是用來控制柱狀圖出現的位置。這就必須自己一直更換不同的數據來看位置對不對。你也可以用一個巨集程式來簡化步驟：<br /><pre><code>%macro bars(reg,h,v);<br />goptions horigin=&h in vorigin=&v in;<br />where region eq &reg;<br />vbar time / type=mean discrete sumvar=x raxis=axis10 maxis=axis&reg mean;<br />run;<br />%mend;</code></pre>因此九張柱狀圖就可以用下面短短的程式碼一口氣完成：<br /><pre><code>proc gchart data=tpay;<br />%bars(1,8.75,6.45) %bars(2,7.75,4.90) %bars(3,7.50,3.00) <br />%bars(4,4.00,5.00) %bars(5,6.00,4.25) %bars(6,6.10,2.50) <br />%bars(7,4.00,1.85) %bars(8,1.95,4.20) %bars(9,0.00,2.00)<br />quit; </code></pre>最後一個步驟就是要把這些柱狀圖疊到地圖上。我們先用 ODS 開一個空白的 pdf 檔：<br /><pre><code>ods pdf file='z:\map_bars.pdf' notoc startpage=never;</code></pre>後面的 notoc 和 startpage 主要是抑制 pdf 檔的書籤出現以及讓所有圖形都出現在同一頁。然後進行 pdf 檔裡面的環境設定：<br /><pre><code>goptions reset=all ftext='Helvetica/bo' htext=3.5 gunit=pct ctext=blue lfactor=3;</code></pre><br />這些都是來設定字形、字色還有字體大小。這樣就大功告成了。整個流程可以用下列四步驟來簡單描述：<br /><br />步驟一：用 GOPTIONS 進行圖形整體環境設定。<br />步驟二：用 ODS 開啟一個空白的 pdf 檔。<br />步驟三：用 PROC GMAP 畫地圖。<br />步驟四：用 PROC GCHART 畫柱狀圖<br /><br />最後用 ODS CLOSE 把整個圖形打包起來。<br /><br /><b><span class="Apple-style-span" style="font-size: large;">REFERENCES & RECOMMENDED READING</span></b><br /><br />以下是一些不錯的線上參考資料：<br /><br /><a href="http://support.sas.com/documentation/onlinedoc/index.html"><!-- m --><a class="postlink" href="http://support.sas.com/documentation/onlinedoc/index.html">http://support.sas.com/documentation/on ... index.html</a><!-- m --></a><br /><a href="http://support.sas.com/rnd/datavisualization/mapsonline/html/home.html"><!-- m --><a class="postlink" href="http://support.sas.com/rnd/datavisualization/mapsonline/html/home.html">http://support.sas.com/rnd/datavisualiz ... /home.html</a><!-- m --></a><br /><a href="http://support.sas.com/rnd/papers"><!-- m --><a class="postlink" href="http://support.sas.com/rnd/papers">http://support.sas.com/rnd/papers</a><!-- m --></a><br /><a href="http://support.sas.com/samples"><!-- m --><a class="postlink" href="http://support.sas.com/samples">http://support.sas.com/samples</a><!-- m --></a><br /><a href="http://ftp.sas.com/techsup/download/sample/graph/other-color.html"><!-- m --><a class="postlink" href="http://ftp.sas.com/techsup/download/sample/graph/other-color.html">http://ftp.sas.com/techsup/download/sam ... color.html</a><!-- m --></a><br /><a href="http://arcdata.esri.com/data/tiger2000/tiger_download.cfm"><!-- m --><a class="postlink" href="http://arcdata.esri.com/data/tiger2000/tiger_download.cfm">http://arcdata.esri.com/data/tiger2000/ ... wnload.cfm</a><!-- m --></a><br /><a href="http://www.census.gov/geo/www/cob/index.html"><!-- m --><a class="postlink" href="http://www.census.gov/geo/www/cob/index.html">http://www.census.gov/geo/www/cob/index.html</a><!-- m --></a><br /><a href="http://www.cdc.gov/epiinfo/shape.htm"><!-- m --><a class="postlink" href="http://www.cdc.gov/epiinfo/shape.htm">http://www.cdc.gov/epiinfo/shape.htm</a><!-- m --></a><br /><a href="http://research.cip.cgiar.org/gis/index.php"><!-- m --><a class="postlink" href="http://research.cip.cgiar.org/gis/index.php">http://research.cip.cgiar.org/gis/index.php</a><!-- m --></a><br /><a href="http://support.sas.com/techsup/technote/ts688/ts688.html"><!-- m --><a class="postlink" href="http://support.sas.com/techsup/technote/ts688/ts688.html">http://support.sas.com/techsup/technote ... ts688.html</a><!-- m --></a><br /><a href="http://www.colorbrewer.org/"><!-- m --><a class="postlink" href="http://www.colorbrewer.org">http://www.colorbrewer.org</a><!-- m --></a><br /><a href="http://robslink.com/SAS/Home.htm"><!-- m --><a class="postlink" href="http://robslink.com/SAS/Home.htm">http://robslink.com/SAS/Home.htm</a><!-- m --></a><br /><br /><b><span class="Apple-style-span" style="font-size: large;">CONTACT INFORMATION</span></b><br />Your comments and questions are valued and encouraged. Contact the authors at:<br />Robert Allison<br /><a href="mailto:robert.allison@sas.com"><!-- e --><a href="mailto:robert.allison@sas.com">robert.allison@sas.com</a><!-- e --></a><br />Louise Hadden<br /><a href="mailto:louise_hadden@abtassoc.com"><!-- e --><a href="mailto:louise_hadden@abtassoc.com">louise_hadden@abtassoc.com</a><!-- e --></a><br />Mike Zdeb<br /><a href="mailto:msz03@albany.edu"><!-- e --><a href="mailto:msz03@albany.edu">msz03@albany.edu</a><!-- e --></a><div class="blogger-post-footer"><img width='1' height='1' src='https://blogger.googleusercontent.com/tracker/6268919072942670865-2375885139926829096?l=sugiclub.blogspot.com' alt='' /></div>