分析方法1

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

rm(list = ls())
library(readr);library(ggtern);library(rgdal);library(raster);library(sp);library(maptools);library(GISTools);library(magrittr);library(dplyr);library(imputeTS);library(ggplot2)

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 ")

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

list_move <- c()
for(i in 1:79){
  file_PSI <- paste0("/Volumes/QuadX/Evan/Project/TWD97_PSI/twd97/ps_u-dm_twd97.",as.character(i),".xy")
  print(file_PSI)
  GPS_move <- read_delim(file_PSI, delim = " ", col_names = c("x", "y", "z"))
  list_move <- append(list_move, assign(paste0("GPS_move_total-",as.character(i)),SpatialPointsDataFrame(GPS_move[,1:2],GPS_move, proj4string = TWD97)))
}


for(i in 1:26){
  file_PSI_pre <- paste0("/Volumes/QuadX/Evan/project_pre/pre_cut/ps_u-dm.",as.character(i),".xy")
  print(file_PSI_pre)
  GPS_move <- read_delim(file_PSI_pre, delim = " ", col_names = c("x", "y", "z"))
  list_move <- append(list_move, assign(paste0("GPS_move_total_pre-",as.character(i)),spTransform(SpatialPointsDataFrame(GPS_move[,1:2],GPS_move, proj4string = wgs84), TWD97)))
}

天數改日期

days <- read_delim("/Volumes/QuadX/Evan/pose_seismic/days_PSI.txt", delim = " ", col_names = c("day"))
days[,2] <- days[,1]-735894
days[,3] <- days$day.1+as.Date("2014-10-22")
dayy <- data.frame(days[,3])

找出站點範圍，找出時間序列，output以便日後使用

for(i in 2:74){
  stop <- GPS_station_TW2[i,]
  Time_series <- data.frame(x=as.double(),
                            y=as.double())
  for(j in 1:105){
    data_select <- list_move[[j]]
    range_nr <- gWithinDistance(stop, data_select, 500, byid=T)
    mean_move <- mean(subset(data_select@data$z,range_nr))
    Time_series <- rbind(Time_series,mean_move)
  }
  Time_series <- cbind(dayy,Time_series)
  out_file <- paste0("/Volumes/QuadX/Evan/gps_series_PSI/",GPS_station[i,4],".txt")
  write.table(Time_series,file=out_file,sep=" ",row.names = F,col.names = F)
}

合併GPS_LOS

目的是要檢核ＰＳＩ的成果，利用ＧＰＳ作為檢核工具，比較兩者的模式是否一致。

rm(list = ls())
library(readr);library(ggtern);library(rgdal);library(raster);library(sp);library(maptools);library(GISTools);library(magrittr);library(dplyr);library(imputeTS);library(ggplot2)



PSI_LOS <- list.files("/Volumes/QuadX/Evan/pose_seismic/PSI_fitting/PSI_series")
fitting_results <- list.files("/Volumes/QuadX/Evan/pose_seismic/GPS_fitting_subset")

#3-dimension to LOS
a=(-12/180)*pi
b=(32.9/180)*pi
every_GPS_LOS <- read_csv("/Volumes/QuadX/Evan/pose_seismic/PSI_fitting/ZIP.txt", col_names = "days")
for(i in 1:74){
  if(is.character(fitting_results[i]) == F){
    next
  }else{
    in_file <- paste0("/Volumes/QuadX/Evan/pose_seismic/GPS_fitting_subset/",fitting_results[i])
  }
  files <- read_csv(in_file, col_names = c("days", "N", "E", "U"))
  #找出時間段
  files <-subset(files,files$days >= 2014)
  files$LOS <- ((-sin(b)*cos(a))*(files$E)+(sin(b)*sin(a))*(files$N)+(cos(b))*(files$U))
  zip <- data.frame(cbind(files$days,files$LOS))
  names(zip)<- c("days",sub(".csv","",paste0(fitting_results[i])))
  every_GPS_LOS <- left_join(every_GPS_LOS,zip)
}

#GPS justify
GS44  <- data.frame(every_GPS_LOS$GS44)
GS44_fill <- na.interpolation(GS44, option = "linear")
GS44_fill_1 <- GS44_fill
for(i in 1:73){
  GS44_fill_1 <- cbind(GS44_fill_1,GS44_fill)
}
every_GPS_LOS_new <- every_GPS_LOS
every_GPS_LOS_new[,2:75]  <- every_GPS_LOS_new[,2:75] - GS44_fill_1

every_GPS_LOS_new_mean <- data.frame(t(data.frame(colMeans(every_GPS_LOS_new[,2:75], na.rm = T, dims = 1))))
every_GPS_LOS_new_mean_1 <- every_GPS_LOS_new_mean
for(i in 1:1687){
  every_GPS_LOS_new_mean_1 <- rbind(every_GPS_LOS_new_mean_1,every_GPS_LOS_new_mean)
}

every_GPS_LOS_new[,2:75] <- every_GPS_LOS_new[,2:75]-every_GPS_LOS_new_mean_1

every_GPS_LOS_new$dayss <- as.numeric(rownames(every_GPS_LOS_new))
every_GPS_LOS_new$daydate <- every_GPS_LOS_new$dayss+as.Date("2013-12-31")

#作圖
for (i in 1:74) {
  nom <- i+1
  in_file <- paste0("/Volumes/QuadX/Evan/pose_seismic/PSI_fitting/PSI_series/",PSI_LOS[i])
  files <- read_delim(in_file,delim = " " ,col_names = c("days", "LOS"))
  mean_post_GPS<- as.numeric(colMeans(every_GPS_LOS_new[nom][677:1688,],na.rm = T ))
  mean_pre_GPS<- as.numeric(colMeans(every_GPS_LOS_new[nom][1:676,],na.rm = T ))
  mean_post_PSI <- mean(files$LOS[1:79])
  mean_pre_PSI <- mean(files$LOS[80:105])
  mean_post <- mean_post_PSI-mean_post_GPS
  mean_pre <- mean_pre_PSI-mean_pre_GPS
  post <- files$LOS[1:79]-mean_post
  pre <- files$LOS[80:105]-mean_pre
  file_real <- files
  PSI <- append(pre, post)
  file_real$LOS <- PSI
  
  GPS.data <- data.frame(every_GPS_LOS_new$daydate,every_GPS_LOS_new[nom])
  colnames(GPS.data) <- c("days", "GPS")
  GPS_LOS.data <- left_join(GPS.data,file_real, by= "days")
  GPS_LOS.data$days <- as.Date(GPS_LOS.data$days)
  #out_file <- paste0("/Volumes/QuadX/Evan/pose_seismic/PSI_fitting/PSI_GPS_series/",PSI_LOS[i])
  #write.table(Time_series,file=out_file,sep=" ",row.names = F,col.names = F) #輸出
  GPS_LOS.data$pre <- GPS_LOS.data$LOS
  GPS_LOS.data$pre[767:1688] <- NA
  GPS_LOS.data$post <- GPS_LOS.data$LOS
  GPS_LOS.data$post[1:767] <- NA
  name <- substring(PSI_LOS[i], 1, 4)
  title <- paste0("GPS LOS and PSI LOS time series at ",name)
  print(ggplot()+theme_set(theme_bw()) + geom_point(data = GPS_LOS.data, aes(x= days, y=GPS), color = "red") + geom_point(data = GPS_LOS.data, aes(x= days, y=LOS), color = "blue") + ggtitle(title)+xlab("date")+ylab("LOS")+ geom_smooth(data = GPS_LOS.data, aes(x= days, y=pre), method = lm,color="#51A8DD", fill="#51A8DD")+ geom_smooth(data = GPS_LOS.data, aes(x= days, y=post), method = lm,color="#51A8DD", fill="#51A8DD"))
}