空間分析期中考參考答案

前置作業

rm(list=ls())
setwd("D:/SA_Mid-term")
library(rgdal)#read shp
library(spatstat) #ppp
library(sp)#座標系統設定
library(maptools)#window
library(splancs)#as.points
crime <- readOGR(dsn = "Data", layer = "crime1", encoding="big5")
TP <- readOGR(dsn = "Data", layer = "Taipei_district", encoding="big5")
climate<-read.csv("Data/Climate_data.csv")

1. 台北市各行政區「犯罪率」的定義：該區犯罪個案總數 / 該區人口數。請繪製以下兩小題的地圖： (1). 繪製台北行政區犯罪率的面量圖Choropleth Map （需以Equal Interval方式，將屬性分3層，並提供圖例說明）。(10%)

首先進行資料處理

#檢查CRS是否相同
proj4string(TP)

## [1] "+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_defs"

proj4string(crime)

## [1] "+proj=tmerc +lat_0=0 +lon_0=121 +k=0.9999 +x_0=250000 +y_0=0 +ellps=GRS67 +units=m +no_defs"

#轉換CRS
crime2 <- spTransform(crime, CRS("+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_defs"))

#spatial join
crime2<-crime2[TP, ] # selecting pts over polygon
crime_agg<-aggregate(x = crime2["ID"], by = TP, FUN = length) # calculating the counts
TP$crime_Count <- crime_agg$ID
TP$crime_rate<-TP$crime_Count/TP$CENSUS*100

接著產生equal interval

x1<-TP$crime_rate
interval<-max(x1)-min(x1)
breaks=c(0,interval/3,interval*2/3,interval)+min(x1)

面量圖色碼與分組

greys <- paste0("grey", c(seq(from = 100, to = 0, by = round(-100/3)))[1:3])
q1 <- cut(x1, breaks , include.lowest=T,right=F)
clr <- as.character(factor(q1, labels = greys))

設定小數位數，以供legend使用

breaks<-sprintf("%.4f",breaks)
plot(TP,col=clr,main="台北市犯罪率分布")
LenTEXT1<-c(paste0("<",breaks[2],"%"),paste0(breaks[3],"%","-",breaks[4],"%"),paste0(">",breaks[3],"%"))
legend(legend = LenTEXT1, fill = greys, "bottomright",cex=0.6,title="犯罪率")

(2). 定義若某區的犯罪率高於台北市各區犯罪率分佈的Q3 (第三個四分位數，upper quantile)，即為高犯罪區；若低於台北市各區犯罪率分佈的Q1 (第一個四分位數，lower quantile)，即為低犯罪區。將符合高犯罪區標示成紅色，符合低犯罪區標示成綠色，其餘標示為灰色，繪製「台北市犯罪地圖」(含圖例說明)。(15%)

###Q1-2  quantile
breaks1<-c(quantile(x1,probs=c(0,0.25,0.75,1)))
q1<<-cut(x1,breaks1,include.lowest=T,right=F)
#設定legend顏色
color <- c("green","gray","red")
#圖例色碼
clr <- as.character(factor(q1, labels = color))
plot(TP,col=clr,main="台北市犯罪地圖")
#設計legend
LenTEXT1<-c("低犯罪區","一般","高犯罪區")
legend(legend = LenTEXT1, fill = color, "bottomright",cex=0.8)

針對台北市的犯罪資料，計算犯罪事件兩兩之間的直線距離，找出兩事件之間最短距離為何？這兩起事件位於哪一個行政區？(10%)

讀取時空矩陣函數

stmatrix<-function(pxy,pt){
  xymatrix<-matrix(0,length(pxy[,1]),length(pxy[,1]))
  tmatrix<-matrix(0,length(pt),length(pt))
  for(i in 1:length(pxy[,1])){
    for(j in i:length(pxy[,1])){
      x1<-pxy[i,1]-pxy[j,1]
      y1<-pxy[i,2]-pxy[j,2]
      xymatrix[i,j]<-(x1^2+y1^2)^(1/2)
      t1<-abs(pt[i]-pt[j])
      tmatrix[i,j]<-t1
      for (j in 1:i){
        tmatrix[i,j]<-NA
        xymatrix[i,j]<-NA
      }
    }
  }
  return(list(xymatrix,tmatrix))
}

產生crime的時空矩陣

cmatrix<-stmatrix(crime2@coords,crime2$Time)
#最短距離
min(cmatrix[[1]],na.rm=T)

## [1] 43.5779

#找最近點
minp<-c()
index<-max(cmatrix[[1]],na.rm=T)
for(i in 1:116){
  for(j in 1:116){
    if(!is.na(cmatrix[[1]][i,j])){
      if(cmatrix[[1]][i,j] < index){
        minp[1:2]<-c(i,j)
        index<-cmatrix[[1]][i,j]
      }}   
  }}
#顯示最近的點ID
minp

## [1] 85 86

plot(TP)
plot(crime[85:86,],add=T)

結果顯示最短距離為43.5779，最近的兩點(ID=85,86)在松山區。

針對台北市的犯罪資料，以發生時間在第11週(含)以前的個案，作為參考組(reference group)，利用k function來評估在12週之後新發生的個案，再調整參考組的個案分佈之後，是否呈現顯著的群聚型態？(15%)

先將所有資料依週次分組

sel<-crime$Time<=11
sel2<-crime$Time>11
crime2.r<-crime2[sel,] #參照組
crime2.c<-crime2[sel2,]

資料處理(ppp格式)

crime.o<-as(TP,"owin")
crime.cppp<-as(crime2.c,"ppp")
crime.rppp<-as(crime2.r,"ppp")
cpattern<-ppp(crime.cppp$x,crime.cppp$y, window=crime.o)
rpattern<-ppp(crime.rppp$x,crime.rppp$y, window=crime.o)
#繪製一般的k-function
cKF<-Kest(cpattern)
cKFE<-envelope(cpattern, Kest, nsim=99)
plot(cKFE,main="crimes after week 12(unadjusted)")

運用參照組評估非均值過程，並得到調整後的結果

mydense<-density(rpattern)
Kenv2<-envelope(cpattern, Kinhom, simulate = expression(rpoispp(mydense)), nrank=2, nsim=99)# nrank=2  95 C.I

## Generating 99 simulations by evaluating expression  ...
## 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15,
## 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30,
## 31, 32, 33, 34, 35, 36, 37, 38, 39, 40, 41, 42, 43, 44, 45,
## 46, 47, 48, 49, 50, 51, 52, 53, 54, 55, 56, 57, 58, 59, 60,
## 61, 62, 63, 64, 65, 66, 67, 68, 69, 70, 71, 72, 73, 74, 75,
## 76, 77, 78, 79, 80, 81, 82, 83, 84, 85, 86, 87, 88, 89, 90,
## 91, 92, 93, 94, 95, 96, 97, 98,  99.
## 
## Done.

plot(Kenv2,main = "Inhomogeneous K function")

##       lty col   key                            label
## obs     1   1   obs italic({hat(K)[inhom]^{obs}}(r))
## mmean   2   2 mmean       italic({bar(K)[inhom]}(r))
## hi      1   8    hi  italic({hat(K)[inhom]^{hi}}(r))
## lo      1   8    lo  italic({hat(K)[inhom]^{lo}}(r))
##                                                        meaning
## obs             observed value of K[inhom](r) for data pattern
## mmean              sample mean of K[inhom](r) from simulations
## hi    upper pointwise envelope of K[inhom](r) from simulations
## lo    lower pointwise envelope of K[inhom](r) from simulations

** 4. 某警察分局以台北市南區（包括：文山、萬華、中正、大安、信義與南港區）為犯罪分析的研究區，以該研究區內的這些犯罪事件，進行以下的時空分析： (1). 計算這些犯罪事件之間的space-distance and time-distance，請繪製Space-time Scatter Plot。 (10%) **

#轉投影
crime = spTransform(crime,CRS(proj4string(TP)))

先建立計算空間距離和時間距離的function

dd = function(x1,y1,x2,y2)
{
  distance = sqrt((x1-x2)^2+(y1-y2)^2)
  return (distance)
}
tt = function(t1,t2)
{
  td = abs(t1-t2)
  return (td)
}

進行資料選取

Ct = crime$TOWN
Dt = TP$TOWN
selC = Ct == "中正區"| Ct == "大安區" | Ct ==  "信義區" | Ct ==  "萬華區" |  Ct == "文山區" |  Ct == "南港區"
selD = Dt == "中正區"| Dt == "大安區" | Dt ==  "信義區" | Dt ==  "萬華區" |  Dt == "文山區" |  Dt == "南港區"
crime4 = crime[selC,]
TP4 = TP[selD,]

產生矩陣

smatrix4 = matrix(0,nrow=nrow(crime4),ncol=nrow(crime4))
tmatrix4 = matrix(0,nrow=nrow(crime4),ncol=nrow(crime4))
for(i in 1:nrow(crime4)) 
{
  for(j in i:nrow(crime4))
  {
    smatrix4[i,j] = dd(crime4@data[i,4],crime4@data[i,5],crime4@data[j,4],crime4@data[j,5])
    tmatrix4[i,j] = tt(crime4@data[i,6],crime4@data[j,6])
  }
}

畫散佈圖

plot(smatrix4,tmatrix4,xlim=range(smatrix4),ylim=range(tmatrix4),type="p",main="Space-Time association plot",xlab="Space_Distance",ylab="Time_Distance",pch=1,col="red")

(2). 利用space-time K function，評估這些犯罪事件是否呈現顯著的時空群聚？（需呈現D0(d,t)與Standardized residuals的計算結果）(15%)

擷取選取行政區的最大邊界

TaipeiS = aggregate(TP4, by=list(TP4@data$COUNTY), FUN=length) 
boundary4 = TaipeiS@polygons[[1]]@Polygons[[1]]@coords 
boundary4 = as.points(boundary4[,1], boundary4[,2])

計算空間半徑 1000 到 10000(1000) 時間半徑 1 ~ 20(2)

spk = stkhat(as.points(crime4$X, crime4$Y), crime4@data$Time, boundary4, range(crime4@data$Time), seq(1000,10000,1000),seq(0, 20, 2))
s = matrix(spk$ks)
t = matrix(spk$kt)
st = s %*% t(t)
D = spk$kst - st
persp(spk$s, spk$t, D, theta=-30, phi = 15, expand = 0.5, xlab="spatial distance", ylab="temporal distance", zlab="D", ticktype ="detailed" )

D0 = spk$kst/st-1.0
persp(spk$s, spk$t, D0, theta=-30,phi = 15, expand = 0.5, xlab="spatial distance", ylab="temporal distance", zlab="D0", ticktype ="detailed" )

散佈圖

se= stsecal(as.points(crime4$X, crime4$Y), crime4@data$Time, boundary4, range(crime4@data$Time), seq(1000,10000,1000),seq(0, 20, 2))
Res= D / se
plot( st, Res ,ylim=c(-3, 3),xlab="K1(s)K2(t)", ylab=" standardized residuals R(s,t)")
abline(h=c(-2,2), lty=2)

根據試題提供的氣象資料，進行以下統計圖表繪製： (1). 計算2013年1月到12月的月平均氣溫及其標準差，並繪製合適的統計圖表來呈現這兩個統計數據的逐月趨勢。(10%)

Tmm = c() #月均溫
Tmv = c() #月溫標準差
sel4<-list() #欄位查詢
cli.M<-list() #每月資料

排除有誤的資料並按照月份分開同時算月均溫和標準差

  for(j in 1:length(climate$Year))
  {
  ifelse(climate$Temperature[j]<(-1000),climate$Temperature[j]<-NA,climate$Temperature[j]<-climate$Temperature[j])
  }

  for(i in 1:12){
      sel4[[i]]<-climate$Month==i
      cli.M[[i]]<-climate[sel4[[i]],]
      Tmm[i] = mean(cli.M[[i]]$Temperature,na.rm=T)
      Tmv[i] = sd(cli.M[[i]]$Temperature,na.rm=T)
  }
Tmm

##  [1] 15.58524 17.26587 18.55216 20.65334 25.38448 28.51632 29.10094
##  [8] 29.10176 27.37149 23.59704 20.96457 16.17439

Tmv

##  [1] 2.701512 2.901808 3.689734 3.469699 3.300154 2.444232 2.597699
##  [8] 2.832235 3.051021 2.599147 4.157043 3.566506

plot(Tmm, xlim=range(1,12), ylim=c(0,40),xlab="月份", ylab="溫度", main="2013年逐月氣溫趨勢", type="l")
lines(Tmm+Tmv, col="red")
lines(Tmm-Tmv, col="blue")
legend(1,40,c("均溫+標準差","平均氣溫","均溫-標準差"),lty=c(1,1),col=c("red","black","blue"))

(2). 繪製2013年夏季溫度的直方圖 (histogram)，並套疊一常態分佈曲線（該常態分佈的平均值：夏季溫度平均值；標準差：夏季溫度標準差） (15%) 說明：「夏季」的定義是從夏至到秋分，亦即 6/21-9/23

head(climate) #先觀察資料型態，找可用的欄位

##   Site_No       Time Year Month Day hour Temperature Relative_Humidity
## 1  466900 2013010101 2013     1   1    1        11.3                87
## 2  466900 2013010102 2013     1   1    2        11.3                89
## 3  466900 2013010103 2013     1   1    3        11.6                90
## 4  466900 2013010104 2013     1   1    4        11.8                89
## 5  466900 2013010105 2013     1   1    5        12.2                88
## 6  466900 2013010106 2013     1   1    6        12.1                90

sel5<-climate$Time>=2013062100 & climate$Time<2013092400
Summerdata<-climate[sel5,]$Temperature
hist(Summerdata, col="grey", breaks=seq(20,40,by=1), freq=FALSE, xlab="溫度")
curve(dnorm(x, mean=mean(Summerdata,na.rm=T), sd=sd(Summerdata,na.rm=T)), add=TRUE, col="red", lwd=2)

空間分析期中考 參考答案

空間分析期中考參考答案