PSInSAR資料依GS44平移校正，需先求出GS44的平均速度。

pre <- read.csv("/Volumes/QuadX/Evan/pose_seismic/gmp_plot/ps_mean_v_pre_GS44.txt", sep = " ", header = F)
post <- read.csv("/Volumes/QuadX/Evan/pose_seismic/gmp_plot/ps_mean_v_post_GS44.txt", sep = " ", header = F)

GS44_pre_mean <- colMeans(pre[3])
GS44_post_mean <- colMeans(post[3])

使用 awk 高速完成了。程式碼如下：

  awk '{if ($1<=120.4007 && $1>=120.4001 && $2<=23.2225 && $2>=22.2219) print $1, $2, $3}' /Volumes/QuadX/Evan/pose_seismic/gmp_plot/diff.txt > /Volumes/QuadX/Evan/pose_seismic/diff_gs44.txt;

GPS修正

gpsss <- read_delim("/Volumes/QuadX/Evan/pose_seismic/gmp_plot/gpsss", delim = " ", col_names = c("x", "y", "z","site"))
gpsss_pre <- read_delim("/Volumes/QuadX/Evan/pose_seismic/gmp_plot/gpsss_pre", delim = " ", col_names = c("x", "y", "z","site"))
diff_gps <- gpsss
diff_gps[,3] <- gpsss[,3]-gpsss_pre[,3]
#write.table(diff_gps,file="/Volumes/QuadX/Evan/pose_seismic/gmp_plot/diff_gps.txt",sep=" ",row.names = F,col.names = F)
shift <-as.numeric(as.character(gpsss_pre[28,3]))
gpsss_pre_shift <- gpsss_pre
gpsss_pre_shift[,3]<- gpsss_pre_shift[,3]-shift
#write.table(gpsss_pre_shift,file="/Volumes/QuadX/Evan/pose_seismic/gmp_plot/gpsss_pre_shift.txt",sep=" ",row.names = F,col.names = F)
shift <-as.numeric(as.character(gpsss[28,3]))
gpsss_shift <- gpsss
gpsss_shift[,3]<- gpsss_shift[,3]-shift
#write.table(gpsss_shift,file="/Volumes/QuadX/Evan/pose_seismic/gmp_plot/gpsss_pre_shift.txt",sep=" ",row.names = F,col.names = F)

讀檔，將文字檔轉成shp檔，改成TWD97座標

rm(list = ls())
library(readr);library(ggtern);library(rgdal);library(raster);library(sp);library(maptools);library(GISTools);library(gstat);library(magrittr);library(akima);library(dismo);library(fields);library(spdep);library(RANN)

wgs84 <- crs("+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs ")
TWD97 <- 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 ")

pre_eq <- read_delim("/Volumes/QuadX/Evan/pose_seismic/gmp_plot/ps_mean_v_pre_move.txt", delim = " ", col_names = c("x", "y", "z"))
post_eq <- read_delim("/Volumes/QuadX/Evan/pose_seismic/gmp_plot/ps_mean_v_post_move.txt", delim = " ", col_names = c("x", "y", "z"))
pre_eq_84 <- SpatialPointsDataFrame(pre_eq[,1:2],pre_eq, proj4string = wgs84)
post_eq_84 <- SpatialPointsDataFrame(post_eq[,1:2],post_eq, proj4string = wgs84)

pre_eq_TW <- spTransform(pre_eq_84, TWD97)
pre.dataframe <- data.frame(pre_eq_TW@coords,pre_eq_TW@data$z)
post_eq_TW <- spTransform(post_eq_84, TWD97)
post.dataframe <- data.frame(post_eq_TW@coords,post_eq_TW@data$z)

gs44 <- read_delim("/Volumes/QuadX/Evan/pose_seismic/diff_gs44.txt", delim = " ", col_names = c("x", "y", "z"))
v <- colMeans(gs44[,3])

以震前和震後最近的點來計算差值。

stop <- post_eq_TW[1,]
near <- nn2(pre.dataframe,post.dataframe,k=1)
near_value <- pre.dataframe[as.numeric(near$nn.idx[,1]),]
rownames(near_value) <- NULL
near_value[,1:2] <- post.dataframe[,1:2]
diff <- data.frame(post_eq_TW@data[,1:2])
diff[,3]<- post.dataframe[,3]-near_value[,3]-v
#write.table(diff,file="/Volumes/QuadX/Evan/pose_seismic/gmp_plot/diff.txt",sep=" ",row.names = F,col.names = F)

400*400m網格(分析方法2前置作業)

#建立網格

diff <- read_delim("/Volumes/QuadX/Evan/pose_seismic/gmp_plot/diff.txt", delim = " ", col_names = c("x", "y", "z"))
diff_84 <- SpatialPointsDataFrame(diff[,1:2],diff, proj4string = wgs84)
diff_TW <- spTransform(diff_84, TWD97)

raster <- raster()
raster@crs <- TWD97
extent(raster) <- extent(c(153554,215954,2460508,2592108))
ncol(raster) <- 156
nrow(raster) <- 329

diff_rast <- rasterize(diff_TW, raster, diff_TW$z, fun=mean)

#write.table(diff_point,file="/Volumes/LANDWIND/diff_point.csv",sep=",",row.names = F,col.names = F)

分析方法1(6000*6000m網格)

raster1 <- raster()
raster@crs <- TWD97
extent(raster1) <- extent(c(153554,219554,2460508,2592508))
ncol(raster1) <- 11
nrow(raster1) <- 22


pre_rast <- rasterize(pre_eq_TW, raster1, pre_eq_TW$z, fun=mean)
post_rast <- rasterize(post_eq_TW, raster1, post_eq_TW$z, fun=mean)
pre_poly <- rasterToPolygons(pre_rast)
post_poly <- rasterToPolygons(post_rast)

pre_nb <- poly2nb(pre_poly)
pre_nb_w <- nb2listw(pre_nb)
count_pre <- pre_poly$layer
gi <- localG(count_pre, pre_nb_w)
LG <- as.vector(gi)

post_nb <- poly2nb(post_poly)
post_nb_w <- nb2listw(post_nb)
count_post <- post_poly$layer
gi_post <- localG(count_post, post_nb_w)
LG_post <- as.vector(gi_post)

pre_point <- data.frame(rasterToPoints(pre_rast))
pre_point[,3] <- LG
point <- data.frame(202819.713, 2536605.644)
center <- SpatialPointsDataFrame(point,point,proj4string =TWD97)
pre_point_LG <- SpatialPointsDataFrame(pre_point[,1:2],pre_point,proj4string =TWD97)
d0.5 <- gWithinDistance(center,pre_point_LG, 10000, byid=T)
dis0.5 <- sum(subset(pre_point_LG@data$layer,d0.5))
d1 <- gWithinDistance(center,pre_point_LG, 20000, byid=T)
dis1 <- sum(subset(pre_point_LG@data$layer,d1))-dis0.5
d1.5 <- gWithinDistance(center,pre_point_LG, 30000, byid=T)
dis1.5 <- sum(subset(pre_point_LG@data$layer,d1.5))-dis0.5-dis1
d2 <- gWithinDistance(center,pre_point_LG, 40000, byid=T)
dis2 <- sum(subset(pre_point_LG@data$layer,d2))-dis0.5-dis1-dis1.5
d2.5 <- gWithinDistance(center,pre_point_LG, 50000, byid=T)
dis2.5 <- sum(subset(pre_point_LG@data$layer,d2.5))-dis0.5--dis1.5-dis2
d3 <- gWithinDistance(center,pre_point_LG, 60000, byid=T)
dis3 <- sum(subset(pre_point_LG@data$layer,d3))-dis0.5-dis1-dis1.5-dis2-dis2.5

post_point <- data.frame(rasterToPoints(post_rast))
post_point[,3] <- LG_post
post_point_LG <- SpatialPointsDataFrame(post_point[,1:2],post_point,proj4string =TWD97)
d0.5_post <- gWithinDistance(center,post_point_LG, 10000, byid=T)
dis0.5_post <- sum(subset(post_point_LG@data$layer,d0.5_post))
d1_post <- gWithinDistance(center,post_point_LG, 20000, byid=T)
dis1_post <- sum(subset(post_point_LG@data$layer,d1_post))-dis0.5_post
d1.5_post <- gWithinDistance(center,post_point_LG, 30000, byid=T)
dis1.5_post <- sum(subset(post_point_LG@data$layer,d1.5_post))-dis0.5_post-dis1_post
d2_post <- gWithinDistance(center,post_point_LG, 40000, byid=T)
dis2_post <- sum(subset(post_point_LG@data$layer,d2_post))-dis0.5_post-dis1_post-dis1.5_post
d2.5_post <- gWithinDistance(center,post_point_LG, 50000, byid=T)
dis2.5_post <- sum(subset(post_point_LG@data$layer,d2.5_post))-dis0.5_post-dis1_post-dis1.5_post-dis2_post
d3_post <- gWithinDistance(center,post_point_LG, 60000, byid=T)
dis3_post <- sum(subset(post_point_LG@data$layer,d3_post))-dis0.5_post-dis1_post-dis1.5_post-dis2_post-dis2.5_post

movement <- data.frame(distance= c(10000,20000,30000,40000,50000,60000), pre= c(dis0.5,dis1,dis1.5,dis2,dis2.5,dis3), post= c(dis0.5_post,dis1_post,dis1.5_post,dis2_post,dis2.5_post,dis3_post))
ggplot() + geom_line(data=movement,aes(x= distance, y= pre,col ="pre EQ")) + geom_line(data=movement,aes(x= distance, y= post,col="post EQ")) + ggtitle("Gi Value & Distance from Epicenter") +scale_x_continuous(name = "Distance from Epicenter") + scale_y_continuous(name = "Gi statistic excludes sum")

分析方法2

diff_kriging <- read_delim("/Volumes/QuadX/Evan/diff_kriging.txt", delim = " ", col_names = c("x", "y", "z"))
diff_krig <- SpatialPointsDataFrame(diff_kriging[,1:2],diff_kriging, proj4string = TWD97)
diff_krig_raster <- rasterize(diff_krig, raster, diff_krig$z, fun=mean)
diff_krig_raster_poly <- rasterToPolygons(diff_krig_raster)

dkrp_nb <- poly2nb(diff_krig_raster_poly)
dkrp_nb_w <- nb2listw(dkrp_nb, zero.policy = T)
s.count <- diff_krig_raster_poly$layer
LISA <- localmoran(s.count, dkrp_nb_w, zero.policy = T, alternative = "greater")

chk <- s.count- mean(s.count)
zi <- LISA[,4]
quadrant <- c()
quadrant[chk>0 & zi>0] = 1 #H-H
quadrant[chk<0 & zi>0] = 2 #L-L
quadrant[chk>0 & zi<0] = 3 #H-L
quadrant[chk<0 & zi<0] = 4 #L-H

signif = 0.05
quadrant[LISA[,5]> signif & chk>0 & zi>0] =5
colors = c("red", rgb(.95, .95, .95), "lightpink", rgb(.95, .95, .95), "lightpink")
par(family="黑體-繁 中黑")
plot(diff_krig_raster_poly, border="grey", lwd=0.001, col=colors[quadrant], main = "LISA Cluster Map: 震後變形發生區域")
points(202819.713, 2536605.644, pch=22, bg="yellow", col="black", cex=0.8)
legend("bottomright",legend=c("有震後變形","有震後變形但不顯著","無震後變形","震央"), fill=c("red","lightpink", rgb(.95, .95, .95),"yellow" ) ,bty="n",cex=0.7,y.intersp=1,x.intersp=1)