基本地圖繪製
B06208030 何承諭
2019年7月15日
以下範例使用台南市學校點資料、台北市速食店點位(假資料)與全台行政區資料作圖 here
還有一些推薦的套件們在這裡可以找到 here
- 環境建置 ── 套件們
GISTools (intersect, buffer, gdistance等相關函數)
rgdal (讀shp)
sp (投影轉換)library(GISTools) library(rgdal) library(sp) library(dplyr) # 這個非必要
- 讀檔
以下讀檔範例 Popn_TWN2.shp, Schools.shp
setwd("C:/Users/smili/Downloads/Spatial_Analysis") # 讀取.shp檔 shapefile是由很多個檔案(.cpg .shx等等)所組成的 # 讀檔範例 Popn_TWN2.shp, Schools.shp TW <- readOGR(dsn = ".", layer = "Popn_TWN2", encoding="unicode") # 台灣村里界 (面資料) Schools <- readOGR(dsn = ".", layer = "Schools", encoding="unicode") # 台南市學校 (點資料) TN <- subset(TW, COUNTY == "臺南市") # 選出臺南市
以下讀檔範例 station.csv (經緯度)
setwd("C:/Users/smili/Downloads/Spatial_Analysis") file <- read.csv("station.csv") head(file)## stno name Lon Lat ## 1 466850 五分山雷達站 121.7812 25.0712 ## 2 466880 板橋 121.4420 24.9976 ## 3 466900 淡水 121.4489 25.1649 ## 4 466910 鞍部 121.5297 25.1826 ## 5 466920 臺北 121.5149 25.0377 ## 6 466930 竹子湖 121.5445 25.1621# Case 1: Datum = TWD97, Coor = TM2(m) 如果資料是TM2度XY座標，就直接代入以下公式速解(CRS都幫你設定好了) # crs <- CRS("+proj=tmerc +lat_0=0 +lon_0=121 +k=0.9999 +x_0=250000 +y_0=0 +ellps=GRS80 +units=m +no_defs") # point = SpatialPointsDataFrame(coords = data.frame(file$X, points$Y), # data = points@data, # proj4string = crs) # Case 2: Datum = TWD97/WGS84, Coor = Lon/Lat 如果你的座標是經緯度，那也是代公式XD crs <- CRS("+init=epsg:4326") # proj4string of coords 經緯度的 epsg(座標代碼) point <- SpatialPointsDataFrame(coords = data.frame(file$Lon, file$Lat), data = file, proj4string = crs) point <- spTransform(point, TN@proj4string) # 座標轉換 # 測試讀檔是否成功 par(mar = c(0,0,0,0)) # 設置畫布位置 plot(TW) plot(point, add=T, col = 'blue', cex = 1, pch=19)
不管哪一種方式
讀進來後都會是 SpatialPolygonsDataframe (面) 或 SpatialPointsDataframe (點)
這兩種資料格式內部大致包含了 投影方式(Projection) 幾何資料(Geometry) 屬性資料框(Dataframe)
投影方式在兩個以上的資料堆疊時要一致
幾何資料拿來呈現畫圖(行政區外框)
屬性資料拿來操作(把你要的資料Merge進他的欄位裡用來面量圖著色)# 屬性資料框 呼叫方式: @data data <- TN@data head(data)## TOWN_ID TOWN COUNTY_ID COUNTY A0A14_CNT A0A14_M A0A14_F A15A64_CNT ## 327 67000010 新營區 67000 臺南市 10171 5289 4882 55555 ## 328 67000020 鹽水區 67000 臺南市 2528 1326 1202 18051 ## 329 67000030 白河區 67000 臺南市 2475 1268 1207 19016 ## 330 67000040 柳營區 67000 臺南市 1818 960 858 15110 ## 331 67000050 後壁區 67000 臺南市 1905 988 917 15853 ## 332 67000060 東山區 67000 臺南市 1779 926 853 14296 ## A15A64_M A15A64_F A65UP_CNT A65UP_M A65UP_F INFO_TIME ## 327 27681 27874 12132 5565 6567 107Y06M ## 328 9716 8335 4908 2258 2650 107Y06M ## 329 10322 8694 6803 3201 3602 107Y06M ## 330 8110 7000 4343 2034 2309 107Y06M ## 331 8530 7323 5774 2712 3062 107Y06M ## 332 7921 6375 4867 2269 2598 107Y06M# 投影方法 呼叫方式: @proj4string TN@proj4string## CRS arguments: ## +proj=tmerc +lat_0=0 +lon_0=121 +k=0.9999 +x_0=250000 +y_0=0 ## +ellps=GRS80 +units=m +no_defs
- 確認投影方式是否相同
# 確認兩者投影方式是否相同 Schools@proj4string## CRS arguments: ## +proj=tmerc +lat_0=0 +lon_0=121 +k=0.9999 +x_0=250000 +y_0=0 ## +a=6378157.5 +b=6356772.2 +units=m +no_defsTN@proj4string## CRS arguments: ## +proj=tmerc +lat_0=0 +lon_0=121 +k=0.9999 +x_0=250000 +y_0=0 ## +ellps=GRS80 +units=m +no_defs# 發現疑似投影不相同的情形(橢球體不同) Schools <- spTransform(Schools, TN@proj4string) # 把Schools的投影替換為TN的使一致
- 畫出來觀察一下 (疊圖操作)
par(mar = c(0,0,0,0)) # 設置畫布位置 plot(TN, col="grey", border="aliceblue") plot(Schools, add=T, pch=19, col="blue", cex=0.3)
圖資校正
這裡我們可以發現有些點在邊界外，在確認投影方式相同的情形下一般視為誤差
因此我會選擇進行:1. 邊緣校正 (此步驟視情形使用) 2. 忽略位於縣市界線外的學校 (本次我使用)# 1. 邊緣校正 # 不捨棄點資料，但可能包含到其他行政區的部分點資料 # 在原本臺南市邊界上做環域擴大原本行政區邊界 單位為投影方式的單位，此處為公尺 TN.buffer <- gBuffer(TN, width=300) par(mar = c(0,0,0,0)) # 設置畫布位置 plot(TN.buffer, col="grey", border="aliceblue") plot(Schools, add=T, pch=19, col="blue", cex=0.3)
# 2. 忽略界線外的學校 直接取交集 # 優點就是直覺與直接 乾淨俐落 # byid 一般取交集都使用TRUE (用在配對兩個資料集之間的id關聯) Schools.cor <- gIntersection(TN, Schools, byid = T) par(mar = c(0,0,0,0)) # 設置畫布位置 plot(TN, col="grey", border="aliceblue") plot(Schools.cor, add=T, pch=19, col="blue", cex=0.3)
面量圖 Choropleth
把要畫的資料整併進來
# 叫出屬性資料表 df <- TN@data # 新增欄位 AREA(km2) poly.areas()求土地面積 df$AREA <- poly.areas(TN) / 10^6 # 新增欄位 Schools poly.counts()求Schools在每個行政區內的點數總和 df$Schools <- poly.counts(Schools, TN) # 計算學校密度 df$DENSITY <- df$Schools / df$AREA head(df)## TOWN_ID TOWN COUNTY_ID COUNTY A0A14_CNT A0A14_M A0A14_F A15A64_CNT ## 327 67000010 新營區 67000 臺南市 10171 5289 4882 55555 ## 328 67000020 鹽水區 67000 臺南市 2528 1326 1202 18051 ## 329 67000030 白河區 67000 臺南市 2475 1268 1207 19016 ## 330 67000040 柳營區 67000 臺南市 1818 960 858 15110 ## 331 67000050 後壁區 67000 臺南市 1905 988 917 15853 ## 332 67000060 東山區 67000 臺南市 1779 926 853 14296 ## A15A64_M A15A64_F A65UP_CNT A65UP_M A65UP_F INFO_TIME AREA ## 327 27681 27874 12132 5565 6567 107Y06M 38.60263 ## 328 9716 8335 4908 2258 2650 107Y06M 51.58097 ## 329 10322 8694 6803 3201 3602 107Y06M 128.23901 ## 330 8110 7000 4343 2034 2309 107Y06M 61.54153 ## 331 8530 7323 5774 2712 3062 107Y06M 72.41042 ## 332 7921 6375 4867 2269 2598 107Y06M 124.11432 ## Schools DENSITY ## 327 17 0.44038451 ## 328 13 0.25203094 ## 329 14 0.10917115 ## 330 10 0.16249189 ## 331 11 0.15191184 ## 332 10 0.08057088
面量圖上色
# plot par(mar = c(0,0,3,1)) # 調整下左上右邊界 shades = auto.shading(df$DENSITY,n=6, cols = brewer.pal(6, "Greens")) # 把顏色自動分為6層 choropleth(TN, df$DENSITY, shading=shades) # 面量圖 plot(TN, col=rgb(0,0,0,0), add=T, border="cadetblue1") # 只是設邊框顏色，可忽略 map.scale(152025.1,2537352,10000, "50 km",2,1) # 前面為比例尺的位置設置 可以使用locator()來定位 north.arrow(222425,2580695,1500,col= 'lightblue') # 指北針 title('Schools in Tainan') # 標題 choro.legend(208525.1,2555695, shades, title="Density(per km2)") # 圖例
或者你想要拿算出來的密度去ggplot其他圖表也都可以自行操作
空間分析函數的應用
前置作業─資料匯入整併
# 讀檔 setwd("C:/Users/smili/Downloads/Spatial_Analysis") FF <- readOGR(dsn = ".", layer = "Tpe_Fastfood", encoding="utf8") Taipei_Vill <- readOGR(dsn = ".", layer = "Taipei_Vill", encoding="unicode")# 確認投影 FF@proj4string## CRS arguments: ## +proj=tmerc +lat_0=0 +lon_0=121 +k=0.9999 +x_0=250000 +y_0=0 ## +ellps=GRS80 +units=m +no_defsTaipei_Vill@proj4string## CRS arguments: ## +proj=tmerc +lat_0=0 +lon_0=121 +k=0.9999 +x_0=250000 +y_0=0 ## +ellps=GRS80 +units=m +no_defs# 資料處理 FF_MIC = subset(FF, STORE == "MIC") # 取麥當勞 # 觀察資料 par(mar = c(2,2,2,2)) plot(Taipei_Vill, lwd = 0.5, border = "grey50", bg='aliceblue') plot(FF_MIC, add=T, pch=19, col="blue", cex=0.6)
條件篩選
找出各村里以中心點圓心，半徑1000m內涵蓋最多麥當勞的村里
# 利用村里邊界找出幾何中心點 gCentroid centroids <- gCentroid(Taipei_Vill, byid = T, id = rownames(Taipei_Vill)) # 利用村里中心點找出各村里半徑1000m內有多少麥當勞 df.tmp <- gWithinDistance(FF_MIC, centroids, byid=T, dist=1000) %>% data.frame() df.tmp$SUM <- rowSums(df.tmp) # 計算各行總和 (true / false 相當於 1 / 0) head(df.tmp)## X1 X2 X3 X7 X10 X18 X27 X28 X29 X30 X31 X32 ## 0 FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE ## 1 FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE ## 2 FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE ## 3 FALSE FALSE FALSE FALSE FALSE FALSE TRUE FALSE FALSE FALSE FALSE FALSE ## 4 FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE ## 5 FALSE FALSE FALSE FALSE FALSE FALSE TRUE FALSE FALSE FALSE FALSE FALSE ## X33 X34 X35 X36 X37 X38 X39 X40 X41 X42 X43 X44 ## 0 FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE TRUE TRUE FALSE FALSE ## 1 FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE TRUE FALSE FALSE ## 2 FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE TRUE TRUE FALSE FALSE ## 3 FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE TRUE TRUE FALSE FALSE ## 4 FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE TRUE FALSE FALSE FALSE ## 5 FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE ## X45 X46 X47 X48 X49 X50 X51 X52 X53 X54 X55 X56 ## 0 FALSE FALSE TRUE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE ## 1 FALSE FALSE TRUE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE ## 2 FALSE FALSE TRUE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE ## 3 FALSE FALSE TRUE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE ## 4 FALSE FALSE TRUE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE ## 5 FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE ## X57 X58 X59 X60 X61 X62 X63 X64 X65 X66 X67 X68 ## 0 FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE ## 1 FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE ## 2 FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE ## 3 FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE ## 4 FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE ## 5 FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE ## X69 X70 X71 X72 X73 X74 X75 X76 X77 X78 X79 X80 ## 0 FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE ## 1 FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE ## 2 FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE ## 3 FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE ## 4 FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE ## 5 FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE ## X81 X82 X83 X84 X85 X86 X87 X88 X89 X90 X91 X92 ## 0 FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE ## 1 FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE ## 2 FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE ## 3 FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE ## 4 FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE ## 5 FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE ## X93 X94 X95 X96 X97 X98 SUM ## 0 FALSE FALSE FALSE FALSE FALSE FALSE 3 ## 1 FALSE FALSE FALSE FALSE FALSE FALSE 2 ## 2 FALSE FALSE FALSE FALSE FALSE FALSE 3 ## 3 FALSE FALSE FALSE FALSE FALSE FALSE 4 ## 4 FALSE FALSE FALSE FALSE FALSE FALSE 2 ## 5 FALSE FALSE FALSE FALSE FALSE FALSE 1# 找出總和最大(附近最多麥當勞)的村里是第幾個row 用來當作搜尋村里的index index <- c(df.tmp$SUM == max(df.tmp$SUM)) %>% which() # 要找的那一個村里就可以被搜出來ㄌ 因為她就在行政區圖資的第index行 toFind_Vill = Taipei_Vill[index,] # 那圍繞著這個村里的又是哪些麥當勞呢 答案就是在距離 <1000 內的點資料們 # 這邊一樣可以使用gWithinDistance函數回傳有沒有在1000m內(true/false) # 但也可以使用gDistance函數計算距離回傳距離多遠(numeric) 再篩選出距離 <1000m 的row index2 <- gDistance(FF_MIC ,gCentroid(toFind_Vill), byid = T) %>% as.vector() toFind_MIC <- FF_MIC[index2 <= 1000,]
繪製篩選結果
# 接下來就跟之前作法相似 把它畫出來~ par(mar = c(2,2,2,2)) plot(Taipei_Vill, lwd = 0.5, border = "grey50", bg='aliceblue') plot(toFind_Vill, col="yellow", add=T) plot(toFind_MIC, add=T, col="blue", pch=19, cex = 0.6) map.scale(289702.6,2764163,4000, "4km", 4, 1) # TWD97, 4000m, 4km north.arrow(289702.6,2767730,miles2ft(0.2),col= 'lightblue') title("臺北市被麥當勞服務家數最多之村里") text(302796.3, 2771416, labels=toFind_Vill$VILLAGE, pos=2)
