# Week 3 Using R as a GIS 0323
> 1. Spatial intersection
> 2. Buffering & Merging spatial features
> 3. Data join
> 4. Point-in-polygon and area calculations
> 5. Distance analysis
>
## 1. Spatial Intersection
### Step 1. 利用屬性查詢取得指定資料(Area of Interest)
```r=
# index
index <- us_states$STATE_NAME == "Texas" | us_states$STATE_NAME == "New Mexico" |
us_states$STATE_NAME == "Oklahoma" | us_states$STATE_NAME == "Arkansas"
AoI <- us_states[index,] # 指定的資料圖層
```
### Step 2. Intersect using gIntersection() 圖層取交集
```r=
AoI.torn <- gIntersection(AoI, torn, byid = TRUE) # 點圖層，在指定範圍內的龍捲風點位
```
### Step 3. attach attributes 抓資料出來對
處理 index:
```r=
get.X = lapply(I.split , function(x)[1])
get.Y = lapply(I.split , function(x)[2])
X.id = unlist(get.X); Y.id = (unlist get.Y)
```
利用 x,y的id去對原本x,y的資料表抓資料
製作取交集後圖層的資料表：
```r=
XY.Inter = gIntersection(X,Y,byid=T) # 圖層
# 取id
XY.name = strsplit(names(XY.Inter), " ")
X.id = unlist(lapply(XY.name, function(x) x[1]))
Y.id = unlist(lapply(XY.name, function(x) x[2]))
# 建立 XY.Inter 的資料表
XY.data = data.frame(X.id, Y.id, row.names = names(XY.Inter))
XY.Inter = SpatialPolygonsDataFrame(XY.Inter, XY.data, match.ID=F) # row.names = names(XY.Inter) 跟 match.ID=F 選一個寫
```
1. 透過新圖形特徵做運算：
i.e.各截切小區的面積：poly.areas(XY,Inter)
2. 透過原本表格做運算：
i.e.抓取行政區的人口：
* id是數字的話，且從1開始 X$pop\[as.numeric(X.id)]
* id是數字且從0開始 X$pop\[as.numeric(X.id)+1]
* id不是數字格式或資料不連續--->保持字串的格式 X[X.id,]$pop
## 2. Buffering
```r=
AoI.buf = gBuffer(AoI, width = 25000) # 對整個 AOI做 BUFFER
buf.t = gBuffer(georgia2, width = 5000, byid = T, id = gerogia2$Name) # 對喬治亞州的每個郡做 buffer
```
## 3. Merging (dissolve)
```r=
AoI.merge = gUnaryUnion(us_states)
```
```r=
TW.County <- gUnaryUnion(Popn.TWN, as.character(Popn.TWN$COUNTY_ID)) # 屬於同樣id的才會 merge
```
## 4. Point-in-Polygon, data join and area calculataions
### Step 1. Point-in-Polygon using poly.counts()
```r=
torn.count <- poly.counts(torn, us_states)
```
### Step 2-1. create a new table
```r=
n<-49
new_stateid<-c(1:n); new_statename<-c(1:n)
for (i in 1:n) {
new_stateid[i]<-stateid[i]
new_statename[i]<-as.character(us_states$STATE_NAME[as.numeric(stateid[i])])
}
# create new table
dfnew=data.frame(new_stateid,STATE_NAME=new_statename, torn.count)
```
### Step 2-2 Create us_state attribute using left_join()
```r=
us_states@data<- left_join(us_states@data, dfnew)
new_us.attr<-us_states@data
for (i in 1:n) {
if( is.na(us_states$torn.count[i]) ) {us_states$torn.count[i]=0}
}
```
left_join(x, y, by = c("name.x" = "name.y"))
x 的值會被保留，若只在x，則會顯示na；若只在y，該欄會消失
### Step 3. calculating area and density using poly.areas()
```r=
us_states$AREA.KM2<-poly.areas(us_states2) / (1000 * 1000)
```
## 5. Distance Analysis using gDistance() and gWithinDistance()
### 利用gDistance
到圖形中心點的距離：
```r=
centroids = gCentroid(blocks, byid=T, id=rownames(blocks))
distances = ft2miles(gDistance(places, centroids, byid=T))
```
輸出的每一列(row)是每一個中心點
每一行(col)是每一個指定地點(places)
找出到醫院平均距離最短的點：
每個里中心到醫院的平均距離：
```r=
nearest = apply(distances, 1, mean) # 計算每一列的平均值
```
到醫院平均距離最短的里的index：
```r=
nb = which.min(nearest)
nearest[nb] # 該里的平均距離
blocks@data[nb,] # 該里的屬性資料
```
### 利用gWithinDistance()
```r=
distances_2 = gWithinDistance(places, blocks, byid=T, dist = miles2ft(1.2)) # 回傳 TRUE 跟 FALSE 的陣列
```
#### Mapping serviced areas of hospital
```r=
Hosp1 = distances_2[,1]
plot(blocks)
plot(blocks[Hops1,], col="yellow", add = TRUE)
points(places[1,], col="red", pch = 16, cex = 1.2)
```
分析：
1. 研究區內村里總數：
length(blocks)
2. 研究區內hosp1服務村里的總數
length(blocks\[hosp1,])
3. 研究區內hosp1服務村里的總人數
sum(blocks$POP1990\[Hops1])
4. 研究區內hosp1服務村里的白人總數
sum(blocks\$POP1990[Hosp1] * block$P_WHITE\[Hosp]/100)
實習：
Q1
1. 挑出麥當勞的分布
2. 1 KM buffer 的麥當勞連鎖密度，標示連鎖密度最高的麥當勞店家名稱
Q2
1. 各里中心是否涵蓋麥當勞的服務範圍
2. 找出被最多麥當勞服務的里，標示位置和店家
土地利用 & 淹水程度 spatial
Sign in with Wallet
Connect another wallet