R語言與GIS實踐:Shapefile導入與地圖顯示
1. R項目簡介
R是一個開源的統計計算和圖形編程語言和軟件環境，廣泛應用于統計分析和數據挖掘[1]。R項目遵循GNU協議，它的軟件環境源碼由C,Fortran和R寫就。R可以運行在多個平臺，包括各種Unix，Linux發行版本，Window和MacOS上, 國內有中科院、廈門大學、北京交通大學等多個鏡像[2]。R是對S語言的一種實現，由Ross Ihaka 和 Robert Gentleman 所創建，R取之他們名字的首字母。R具有由用戶貢獻的大量的類庫，能夠處理各種科學計算的問題。
2.R與GIS
與地學問題相關的R類庫也十分繁多，包括Spatial data,Maptools,RGdal等。 不過似乎StackOverflow的主站關于R在GIS的問題比GIS分站更為活躍。原因未知……之所以要使用R是因為R能夠以Hadoop集成[3],希望利用R更高效地處理空間聚類問題。雖然Mashout也是一個很好的解決方案，但是感覺用戶更加需要腳本式的問題處理工具。R的GIS的書籍比較少，基本都是以文檔的形式存在，其中Applied spatial data analysis with R是一本挺好的教材。多使用Example命令和Google也是學習R的一個好途徑。
3.Shapefile導入與地圖顯示
R語言中導入Shapefile是非常容易的，可以使用的包，包括RGdal,maptools,PBSmapping等。在這里我們采用maptools和sp兩個類庫來實現地圖數據的導入和顯示。在R中心選擇鏡像并下載。
加載這兩個類庫，輸入如下代碼：
- library(maptools)
- library(sp)
讀取shapefile文件，可以直接輸入文件路徑或者通過file.choose函數來選擇文件，這里采用國家基礎地理信息系統1比400萬數據中的省界國界。
- vent.map <-readShapeSpatial(file.choose())
或者選擇點、線、面類型輸入
- point.mp <- readShapePoints(file.choose())
- line.mp <- readShapeLines(file.choose())
- poly.mp <- readShapePoly(file.choose())
注意，這個命令并不導入Shapefiel的Proj文件，所以導入的是無坐標信息的數據。
- print(proj4string(vent.map))
輸出是：
[1] NA
查看shapefile屬性，添加坐標信息
- summary(vent.map)
輸出是：
- Object of class SpatialLinesDataFrame
- Coordinates:
- min max
- x 73.446960 135.08583
- y 3.408477 53.55793
- Is projected: NA
- proj4string : [NA]
- Data attributes:
- FNODE_ TNODE_ LPOLY_ RPOLY_ LENGTH
- Min. : 1.0 Min. : 1.0 Min. : 1.0 Min. : 1.0 Min. : 0.0040
- st Qu.: 433.0 1st Qu.: 432.0 1st Qu.: 1.0 1st Qu.: 1.0 1st Qu.: 0.0300
- Median : 876.0 Median : 880.0 Median : 1.0 Median : 65.0 Median : 0.0700
- Mean : 927.6 Mean : 929.9 Mean :128.8 Mean :210.8 Mean : 0.5421
- rd Qu.:1344.0 3rd Qu.:1344.0 3rd Qu.:204.0 3rd Qu.:327.0 3rd Qu.: 0.2350
- Max. :2131.0 Max. :2131.0 Max. :925.0 Max. :926.0 Max. :16.7130
- BOU2_4M_ BOU2_4M_ID GBCODE
- Min. : 1 Min. : 1 Min. :26010
- st Qu.: 447 1st Qu.: 219 1st Qu.:26010
- Median : 893 Median :26010 Median :26010
- Mean : 893 Mean :17260 Mean :35583
- rd Qu.:1339 3rd Qu.:26010 3rd Qu.:61010
- Max. :1785 Max. :61156 Max. :99001
添加坐標信息
proj4string(vent.map) <- "+proj=longlat +datum=WGS84"
顯示地圖
plot(vent.map, axes=TRUE, border="gray")
R就會彈出新的窗口顯示:如圖-1所示：
原文鏈接：http://www.cnblogs.com/ventlam/archive/2012/11/15/2771619.html
4.參考文獻
[1]R (programming language) 2012 http://en.wikipedia.org/wiki/R_(programming_language)
[2]The R Project for Statistical Computing 2012
http://www.r-project.org/
[3]RevolutionAnalytics / RHadoop 2012 https://github.com/RevolutionAnalytics/RHadoop
