空間群聚：點資料分析（一）

本文目錄（點擊可跳轉）

• 1. 點資料讀入
  • 1-1. 地理編碼（Geocoding）
• 2. 什麼是空間分佈趨勢（ spatial distribution patterns）
• 3. 空間點型態的熱區分析
  • 3-1. 直接用視覺化呈現
    • 3-1-1. 繪製點資料地圖
    • 3-1-2. 核密度分布地圖（kernel density map）
      • 3-1-2-1. 用 ggmap 套件快速繪製 kernel density map
      • 3-1-2-2. 用 splancs::kernel2d 函式繪製 kernel density map
      • 3-1-2-3. kernel density 的最佳搜尋半徑

參考資料下載

---

1. 點資料讀入

library(magrittr)
library(data.table)

首先讀入今天要使用的點資料們。 tpe_temple.csv 是從台灣寺廟網所抓下來的，臺北市的所有寺廟分佈位置。我們先來看一下檔案內容

data <- read.table("data/tpe_temple.csv", header=TRUE, sep=",") %>% as.data.table
data

##        X             name    deities                            addr
##   1:   1 二二八公園福德宮     土地公         臺北市中正區懷寧街109號
##   2:   2           長慶廟   福德正神          臺北市中正區晉江街34號
##   3:   3         蓮雲禪苑 釋迦牟尼佛       臺北市中正區連雲街74巷5號
##   4:   4           南福宮   天上聖母 臺北市中正區和平西路一段55巷1號
##   5:   5           光華寺   觀音菩薩         臺北市中正區水源路157號
##  ---                                                                
## 264: 265           天恩宮   彌勒古佛   臺北市文山區指南路38巷37之2號
## 265: 266           樟山寺   千手觀音      臺北市文山區老泉街45巷29號
## 266: 267       木柵集應廟   保儀尊王          臺北市文山區保儀路76號
## 267: 268           指南宮   孚佑帝君         臺北市文山區萬壽路115號
## 268: 269       木柵忠順廟   保儀尊王          臺北市文山區中崙街13號

1-1. 地理編碼（Geocoding）

可以看到寺廟的位置是由地址來表示，但是我們必須要將它轉成座標位置才能呈現在地圖上，或者做進一步的分析。

「把地址轉換成經緯度位置」這件事情，叫作「地理編碼（geocoding）」。

雖然 geocoding 聽起來很專業，但其實你我幾乎天天都在用。沒錯，其實就是 google map 的核心技術。

當我們使用 Google map 搜尋等下中午要去吃的餐廳時，我們 key in 的都是餐廳的地址， Google 會先把這個地址在資料庫中比對到相對應的經緯度後，才在地圖上顯示出來，這就是 geocoding 在做的事情。

所以如何將寺廟的地址轉為經緯度位置呢？我們借用 Google 大嬸所提供的 Google Maps Geocoding API 然後我將所有過程包裹成一個 function getLocation 如下：

getLocation = function(addr){
  library(magrittr)
  library(jsonlite)
  library(data.table)
  addr %<>% as.character
  url = 
    paste0("https://maps.googleapis.com/maps/api/geocode/json?address=",
           addr) 
  res = fromJSON(URLencode(url), flatten = TRUE) 
  lat = res$results$geometry.location.lat
  lng = res$results$geometry.location.lng
  result = data.table("addr" = addr,
                      "lat" = lat,
                      "lng" = lng)
  return(result)  
}

寫好之後來應用看看，然後將得到的結果和寺廟分佈位置的檔案存在一起。

table = lapply(data$addr, getLocation) %>% do.call(rbind, .) 

data = merge(data, table, by = "addr", all =F) %>% 
  .[,.("name" = name,
             "deities" = deities,
             "addr" = addr, 
             "lat" = lat, 
             "lng" = lng)] %>% 
  .[!duplicated(., by=c("addr"))]

data

##            name        deities                                 addr
##   1:     圓明宮       瑤池金母  臺北市中山區中山北路二段128巷3之2號
##   2:     圓通巖 南無觀世音菩薩    臺北市中山區中山北路四段1巷1弄3號
##   3: 圓山藥師寺         藥師佛  臺北市中山區中山北路四段1巷1弄6號旁
##   4:     福正宮       福德正神         臺北市中山區中山北路四段67號
##   5:     金剛寺       阿彌陀佛 臺北市中山區中山北路四段71巷2弄臨3號
##  ---                                                               
## 263: 艋舺龍津宮     順正大王公          臺北市萬華區長沙街二段184號
## 264:     廣照宮       飛天大帝               臺北市萬華區長泰街54號
## 265:   妙覺精舍       觀音菩薩           臺北市萬華區雙園街21巷50號
## 266: 台北朝天宮       天上聖母               臺北市萬華區雙園街27號
## 267:     青龍宮       九龍三公           臺北市萬華區雙園街60巷88號
##           lat      lng
##   1: 25.06187 121.5219
##   2: 25.08163 121.5288
##   3: 25.07931 121.5243
##   4: 25.08050 121.5246
##   5: 25.08163 121.5288
##  ---                  
## 263: 25.04130 121.4997
## 264: 25.02332 121.5002
## 265: 25.03057 121.4953
## 266: 25.02952 121.4947
## 267: 25.02880 121.4924

再來讀入台北市邊界資料 tpe_sqr_bnd2.csv

data_bnd = read.table("data/tpe_sqr_bnd.csv", header=TRUE, sep=",") %>% data.table
data_bnd

##       OBJECTID        X       Y
##    1:        1 309234.9 2778760
##    2:        2 309230.9 2778746
##    3:        3 309228.1 2778734
##    4:        4 309224.5 2778717
##    5:        5 309223.3 2778696
##   ---                          
## 4155:     4155 309276.4 2778769
## 4156:     4156 309269.1 2778764
## 4157:     4157 309252.6 2778760
## 4158:     4158 309241.8 2778759
## 4159:     4159 309234.9 2778760

這時候發現，我找到的台北市邊界地圖，看起來是 TWD67 格式？！ 搞資料就是這樣，找到的資料不一定能盡人意，習慣就好， 總之要讓兩種資料的座標格式相同才能繼續玩下去

這邊借用 snexuz 大大所寫好的程式碼再稍做修改，它引用了臺灣生物多樣性資訊入口網提供的 API，可以很方便地把資料作轉換。我選擇把寺廟點位的經緯度座標格式，轉為 TWD67 座標格式，因為 TWD67 座標格式的單位是公尺，在接下來計算距離等等都會比較方便計算。

transl <- function(sourceSystem, targetSystem, x, y){
  library(rvest)
  library(magrittr)
  library(data.table)
  url <- paste('http://taibif.org.tw/BDTools/proj4/convert.php?',
               'source=', sourceSystem,
               '&destination=', targetSystem,
               '&x=', x,'&y=', y, sep='')
  
  result <- read_html(url) %>%
    html_nodes(xpath = '//*[@id="body"]/text()[2]') %>% 
    as.character() %>% strsplit(' ') %>% unlist() 
  
  result <- data.table(lon = result[3] %>% as.numeric, 
                       lat = result[4] %>% as.numeric)
  
  return(result)
}

data.TWD67 = data.table()
for(i in 1:nrow(data)){
  data.TWD67 = rbind(transl(1, 7, data$lng[i], data$lat[i]), data.TWD67)
  # 經測試，如果一次送太多 request ，對方伺服器會不爽噴 error 403
  # 所以這邊要設定，每隔一點時間再送出下一個 request
  # 所以跑這個這邊會花不少時間，可以先去喝杯水。
  Sys.sleep(0.3)
}

如果不想浪費時間跑上面那串，參考資料裡已附上跑完的結果 data2.rds 供大家繼續接下來的任務。

data.TWD67 = readRDS("data/data.TWD67.rds") %>% data.table 

2. 什麼是空間分佈趨勢（ spatial distribution patterns）

在地理研究上，研究目標的空間分佈趨勢是非常非常非常非常非常核心的議題。

隨便舉幾個問題，例如，我們想知道： 台灣人口都聚集在台灣的哪些鄉鎮？ 是不是偏鄉的老年人口比例特別高？ 偏遠地區的 7-11 是不是比全家還要多？ 離島的學校數量是不是比市區還要少？

或者更 advance 一點，加上「時間」的維度： 慣竊是不是有習慣的作案模式？ 新興流行病傳染的模式為何？ 八嘎囧下課之後的移動軌跡和退休老師的移動軌跡有沒有什麼不一樣？

這些隨便想就一大堆的空間問題，最基本的就是，研究者需要釐清， 這些目標的分佈是否呈現「群聚（cluster）」現象， 也就是找出，哪裡是「熱區（hotspot）」就是研究的核心。 然而，如何衡量群聚這件事情， 不僅有多種各有優劣的計算方法、研究尺度問題、 更也得因為資料型態的不同（點、面），而必須選擇不同的方法。

source: 溫在弘等(2010)

3. 空間點型態的熱區分析

對於分佈在空間中的點資料， 怎麼樣是群聚？怎麼樣是分散？又怎麼樣是隨機分布呢？

source: 溫在弘等(2010)

從上面這張圖中，可以明顯看到群聚、分散與隨機分佈所呈現的 pattern 有明顯的差異 順道一提，「均勻分布」是上面的哪一張圖？還是都不是？大家可以自己想一想。

3-1. 直接用視覺化呈現

當我們想要了解資料的分布狀況時，往往會先用視覺化的方法來讓自己有個底， 畫個直方圖、盒鬚圖，大略知道這筆資料是左尾還是右尾？變異多大？有沒有離群值？

在空間分析領域中也是一樣的概念，我們先把資料呈現在地圖上， 先讓我們用眼睛看一看這資料如何分佈，瞭解一下狀況。

3-1-1. 繪製點資料地圖

library(splancs)

TPE_bnd = as.points(data_bnd$X, data_bnd$Y)
Temple_pts = as.points(data.TWD67$lon , data.TWD67$lat)
  
polymap(TPE_bnd, col="gray")
pointmap(Temple_pts, xlab="X", ylab="Y", pch=21, bg="yellow", add=TRUE)

3-1-2. 核密度分布地圖（kernel density map）

好不容易做完基本資料的處理之後， 我們想知道，這些寺廟到底有沒有群聚呢？ 我們試著用「更專業的方法」來呈現台北市寺廟分布的密度差異。

3-1-2-1. 用 ggmap 套件快速繪製 kernel density map

ggmap 是一個專門繪製地圖的套件，它介接了 Google Map、 OpenStreetMap、Stamen Maps、CloudMade Maps 的 APIs ，可以讓使用者透過 R 來下載指定位置的地圖，並且結合 ggplot 套件的應用來完成地圖與資料繪製。

首先繪製底圖，存在 TPEMap 物件裡。

library(ggmap)

TPE_center = as.numeric(c(121.548606, 25.079462))
TPEMap = ggmap(
  get_googlemap(center=TPE_center, 
                scale=2, 
                zoom=12, 
                color = "bw",
                size = c(550, 850),
                maptype = "terrain")
  , extent="normal")

接著透過超強大的 ggplot 繪圖套件，將地圖與資料一層一層的繪製疊合上去。

• geom_density2d ：可以繪製出核密度等高線
• stat_density2d ：進行核密度估計
• scale_fill_gradient ：設定核密度梯度的著色
• scale_alpha ：設定透明度

TPEMap +
  geom_density2d(data = data, aes(x = data$lng, y = lat), size = 0.3) +
  stat_density2d(data = data,
                 aes(x = data$lng, y = data$lat,
                     fill = ..level.., alpha = ..level..), size = 0.01, 
                 bins = 16, geom = "polygon") +
  scale_fill_gradient(low = "blue", high = "red") +
  scale_alpha(range = c(0, 0.3), guide = FALSE) +
  geom_point(aes(x=data$lng, y=lat), data, col="orange", alpha=0.8, size=0.4) 

這個 Kernel Density 的意義是什麼呢？

\[ f(x, y) = {1\over nh^2} \sum_{i=1}^{n}K({d_{i}(x,y)\over h})\]

好，其實我只是想練習一下用 latex 打公式而已。 看不懂沒關係，我也看不懂。

但是簡單講它就是進行以下幾個步驟：

1. Step1: 把研究區域打上均勻網格 讓每一個網格都擁有一個權重 W

2. Step2：對每一個網格設定搜尋半徑（bandwidth） 讓每個網格的權重 \(W = d1w1+d2w2\)
   

   

3. Step3：選擇適合的計算模型 簡單的說，就是要透過哪一種模型來對計算出的 W 們做「平滑化」

模型有很多種選擇，較廣泛使用的是 Quartic 和 Gaussian 以 Quartic function 為例，它的曲線長得像「鐘型」，可以清楚凸顯高峰值與低峰值。 其的公式寫做 \[k = {3\over \pi }(1-{h_{i}^2\over\tau ^2} )\] 而公式裡的 \(\tau\) 就是前面 Step2 中所定義的搜尋半徑（bandwidth）

3-1-2-2. 用 splancs::kernel2d 函式繪製 kernel density map

我們試著捨棄 ggmap 所提供的 kernel density 自動化計算，回到比較傳統的計算方法。雖然畫出來在視覺上比較不那麼美麗，但可以看看，輸入不同參數會造成什麼樣的影響。

由 splancs 套件所提供的 kernel2d() 函式參數如下 kernel2d(pts,poly,h0,nx=20,ny=20,kernel='quartic',quiet=FALSE) 其中 - h0： 輸入 kernel 的 bandwidth - nx、ny：輸入要沿著 x 軸或 y 軸搜尋幾個 points。 - kernel：要使用哪一個 function。在這裡只有 quartic kernel 可以使用

在 R 的實作中要注意的是，我們使用的資料是座標剛才辛苦轉為 TWD67 的寺廟分佈和台北市邊界點資料，而非經緯度資料，因為單位為公尺才能方便做計算。

library(splancs)
par(mfrow=c(1,2))

polymap(TPE_bnd)
Temple.kernel<-kernel2d(Temple_pts, TPE_bnd, h0=400, nx=100, ny=100) 

## Xrange is  295264.9 316374.2 
## Yrange is  2761740 2789392 
## Doing quartic kernel

image(Temple.kernel, add=T)

polymap(TPE_bnd)
Temple.kernel<-kernel2d(Temple_pts, TPE_bnd, h0=3000, nx=100, ny=100) 

## Xrange is  295264.9 316374.2 
## Yrange is  2761740 2789392 
## Doing quartic kernel

image(Temple.kernel, add=T)

從圖中可以明顯看出，當我們輸入不同搜尋半徑 h0 時會有什麼差別。

左圖的搜尋半徑設定為 400 公尺時，看到台北市的寺廟分佈有好幾處出現零星的熱區。

而右圖的搜尋半徑設定為 3000 公尺時，則是出現一塊很大片的熱區在大同、萬華地區，和一比較不明顯的熱區出現在士林一帶

原來搜尋半徑會大大影響到 kernel density 的平滑程度！ 光是 kernel 的搜尋半徑設定不同，在視覺化的判讀上就會出現明顯的不同。 那麼，有沒有「最佳」的搜尋半徑存在呢？

3-1-2-3. kernel density 的最佳搜尋半徑

splancs 套件提供了 mse2d 函式，可以幫助我們尋找最佳 bandwidth。 它的說明文件裡描述，mse2d 可以 Estimate the Mean Square Error for a Kernel Smoothing。

所謂的 Mean Square Error（MSE） 均方差，可以方便地衡量「平均誤差」，以評價數據的變化程度，當 MSE 的值越小，說明預測模型描述實驗數據具有更好的精確度。

Mse2d <- mse2d(Temple_pts, TPE_bnd, nsmse=30, range=5000)
plot(Mse2d$h[1:30], Mse2d$mse[1:30], type="l", 
     xlab="bandwidth", ylab="MSE")

從這張圖中，可以看到 Y 軸的值是 MSE 估計值，根據上面說的，當 MSE 值越小，則此預測模型有較好的精確度，所以我們可以拿到最佳的 bandwidth 為 Mse2d$h[which.min(Mse2d$mse)] 也就是 500 公尺。

將這個結果套用回 kernel2d ，可以得到最佳 kernel density map 如下：

polymap(TPE_bnd)
Temple.kernel<-kernel2d(Temple_pts, TPE_bnd, 
                        h0=Mse2d$h[which.min(Mse2d$mse)], 
                        nx=100, ny=100) 

## Xrange is  295264.9 316374.2 
## Yrange is  2761740 2789392 
## Doing quartic kernel

image(Temple.kernel, add=T)