R/f_spatial_interpolation_R.R

Defines functions spatial_interpolation

Documented in spatial_interpolation

# spatial interpolation -------------------------------------------------------

#' Run data munging script.
#' @export
#' @param type Sensor type
#' @param type_sub Number of sensors considered 
#' @param interpolation Interpolation type
#' @param adjustment Type of adjustment 
#' @param adjustment_sub Subtype of adjustment
spatial_interpolation <- function(type,type_sub,interpolation,adjustment,adjustment_sub) {
  
  if(type=="") stop("Warning: choose measurement type")
  
  # get meta data
  sp::coordinates(position_rg) <- ~coord.X+coord.Y
  sp::proj4string(position_rg) <- sp::CRS("+proj=somerc +lat_0=46.95240555555556 +lon_0=7.439583333333333 +k_0=1 +x_0=600000 +y_0=200000 +ellps=bessel +units=m +no_defs")
  prj.LatLong <- sp::CRS("+proj=longlat +ellps=WGS84 +datum=WGS84")
  position_rg <- sp::spTransform(position_rg, prj.LatLong)
  position_rg <- as.data.frame(position_rg)
  subcatchments@data <- dplyr::left_join(subcatchments@data,subc,by=c("Subname"="ID"))
  subcatchments <- subcatchments[!is.na(subcatchments@data$X),]
  prj.LatLong <- sp::CRS("+proj=longlat +ellps=WGS84 +datum=WGS84")
  subcatchments <- sp::spTransform(subcatchments, prj.LatLong)
  subcatchments.df <- subcatchments@data
  subcatchments.df <- subcatchments.df[,2]
  par0 <- get.ITU.pars(meta$V7,meta$V6)
  w.noise <- meta$length/(par0[,1]*0.03)
  L <- cbind(meta$V10,meta$V11,meta$V12,meta$V13)

  if(interpolation=="_IDW"){
    for(intr in 1:35){
      assign("rg_data",get(load(file = paste0(baseLoc,"/../data-processed/",type,"/",type,type_sub,adjustment,adjustment_sub,"_event_",intr,".RData"))))
      for(s in 1:5){
        assign(paste0("rg_",s,"_coord"),position_rg[s,2:3])
        rg_data <- cbind(rg_data,do.call("rbind", replicate(length(rg_data$Date), get(paste0("rg_",s,"_coord")), simplify = FALSE)))
      }
      colnames(rg_data) <- c("Date","rain.intensity.1","rain.intensity.2","rain.intensity.3","rain.intensity.4","rain.intensity.5","coord.X.1","coord.Y.1","coord.X.2","coord.Y.2","coord.X.3","coord.Y.3","coord.X.4","coord.Y.4","coord.X.5","coord.Y.5")
      rg.idw.matrix <- matrix(,nrow=102,ncol=length(rg_data$Date)+1)
      for(i in 1:length(rg_data$Date)){
        rg.data <- rg_data[i,]
        rg.data <- data.frame(Date=rep(rg.data$Date,5),rain.intensity=c(rg.data$rain.intensity.1,rg.data$rain.intensity.2,rg.data$rain.intensity.3,rg.data$rain.intensity.4,rg.data$rain.intensity.5),coord.X=c(rg.data$coord.X.1,rg.data$coord.X.2,rg.data$coord.X.3,rg.data$coord.X.4,rg.data$coord.X.5),coord.Y=c(rg.data$coord.Y.1,rg.data$coord.Y.2,rg.data$coord.Y.3,rg.data$coord.Y.4,rg.data$coord.Y.5))
        sp::coordinates(rg.data) <- ~coord.X+coord.Y
        sp::proj4string(rg.data) <- sp::CRS("+proj=longlat +ellps=WGS84 +datum=WGS84")
        rg.idw <- gstat::idw(rain.intensity~1,rg.data,newdata=subcatchments,idp=2.0)
        rg.idw <- ggplot2::fortify(rg.idw@data)
        rg.idw.matrix[,i+1] <- rg.idw$var1.pred
        print(i)
      }
      rg.idw.df <- as.data.frame(rg.idw.matrix)
      rg.idw.df$V1 <- as.character(subcatchments.df)
      colnames(rg.idw.df) <- c("ID",as.numeric(as.POSIXct(rg_data$Date)))
      save(rg.idw.df,file=paste0(baseLoc,"/../data-processed/",type,"/",type,type_sub,interpolation,adjustment,adjustment_sub,"_event_",intr,".RData"))
    }
  }
  
  if(interpolation=="_Thiessen"){
    library(maptools)
    library(spatstat)
    for(intr in 1:35){
      assign("rg_data",get(load(file = paste0(baseLoc,"/../data-processed/",type,"/",type,type_sub,adjustment,adjustment_sub,"_event_",intr,".RData"))))
      for(s in 1:5){
        assign(paste0("rg_",s,"_coord"),position_rg[s,2:3])
        rg_data <- cbind(rg_data,do.call("rbind", replicate(length(rg_data$Date), get(paste0("rg_",s,"_coord")), simplify = FALSE)))
      }
      colnames(rg_data) <- c("Date","rain.intensity.1","rain.intensity.2","rain.intensity.3","rain.intensity.4","rain.intensity.5","coord.X.1","coord.Y.1","coord.X.2","coord.Y.2","coord.X.3","coord.Y.3","coord.X.4","coord.Y.4","coord.X.5","coord.Y.5")
      rg.thiessen.matrix <- matrix(,nrow=102,ncol=length(rg_data$Date)+1)
      for(i in 1:length(rg_data$Date)){
        rg.data <- rg_data[i,]
        rg.data <- data.frame(Date=rep(rg.data$Date,5),rain.intensity=c(rg.data$rain.intensity.1,rg.data$rain.intensity.2,rg.data$rain.intensity.3,rg.data$rain.intensity.4,rg.data$rain.intensity.5),coord.X=c(rg.data$coord.X.1,rg.data$coord.X.2,rg.data$coord.X.3,rg.data$coord.X.4,rg.data$coord.X.5),coord.Y=c(rg.data$coord.Y.1,rg.data$coord.Y.2,rg.data$coord.Y.3,rg.data$coord.Y.4,rg.data$coord.Y.5))
        sp::coordinates(rg.data) <- ~coord.X+coord.Y
        sp::proj4string(rg.data) <- sp::CRS("+proj=longlat +ellps=WGS84 +datum=WGS84")
        # masking problem? spatstat::drichlet, maptools::as.ppp?
        th  <- as(dirichlet(as.ppp(rg.data)),"SpatialPolygons")
        sp::proj4string(th) <- sp::proj4string(rg.data)
        th.z <- sp::over(th,rg.data,fn=mean)
        th.spdf  <-  sp::SpatialPolygonsDataFrame(th, th.z)
        th.clp   <- raster::intersect(subcatchments,th.spdf)
        th.clp <- th.clp[!is.na(subcatchments@data$X),]
        rg.thiessen <- ggplot2::fortify(th.clp@data)
        rg.thiessen.matrix[,i+1] <- rg.thiessen$rain.intensity
      }
      rg.thiessen.df <- as.data.frame(rg.thiessen.matrix)
      rg.thiessen.df$V1 <- as.character(subcatchments.df)
      colnames(rg.thiessen.df) <- c("ID",as.numeric(as.POSIXct(rg_data$Date)))
      save(rg.thiessen.df,file=paste0(baseLoc,"/../data-processed/",type,"/",type,type_sub,interpolation,adjustment,adjustment_sub,"_event_",intr,".RData"))
    }
  }
  
  if(interpolation=="_goldsthein"){
    for(intr in 1:35){
      R.5min.date <- get(load(file=paste0(baseLoc,"/../data-processed/",type,"/",type,type_sub,adjustment,adjustment_sub,"_event_",intr,".RData")))
      R.5min <- R.5min.date[,-1]
      CML_goldsthein.5min <- matrix(,nrow=102,ncol=length(R.5min$Sensor_1))
      for (j in 1:length(R.5min$Sensor_1)){
        CML_goldsthein.5min[,j] <- goldsthein(L, R.5min[j,], w.noise, 0.5, subc)
      }
      CML_goldsthein.5min <- cbind(ID=subc$ID,as.data.frame(CML_goldsthein.5min))
      CML_goldsthein.5min$ID <- as.character(CML_goldsthein.5min$ID)
      colnames(CML_goldsthein.5min) <- c("ID",as.numeric(as.POSIXct(R.5min.date$Date)))
      save(CML_goldsthein.5min,file=paste0(baseLoc,"/../data-processed/",type,"/",type,type_sub,interpolation,adjustment,adjustment_sub,"_event_",intr,".RData"))
    }
  }
  
  if(type == "Radar"){
    
    for (intr in 1:length(events$Event.Nr)){
      assign("tmp",get(load(file=paste0(baseLoc,"/../data-processed/",type,"/",type,type_sub,adjustment,adjustment_sub,"_event_",intr,".RData"))))
      Date <- tmp$Date
      tmp <- tmp[,-1]
      tmp <- t(tmp)
      df <- data.frame(ID=radar_meta$ID,Sensor=radar_meta$Sensor)
      Sensor_1 <- tmp["Sensor_1",]
      Sensor_2 <- tmp["Sensor_2",]
      Sensor_3 <- tmp["Sensor_3",]
      matrix <- matrix(nrow=length(df$ID),ncol=length(Date))
      for(ii in 1:length(df$ID)){
        if(df$Sensor[ii]=="Sensor_1") matrix[ii,] <- Sensor_1
        if(df$Sensor[ii]=="Sensor_2") matrix[ii,] <- Sensor_2
        if(df$Sensor[ii]=="Sensor_3") matrix[ii,] <- Sensor_3
      }
      df <- cbind(df,as.data.frame(matrix))
      df <- df[,-2]
      colnames(df) <- c("ID",as.numeric(as.POSIXct(Date)))
      df[,1] <- as.character(df[,1])
      save(df,file=paste0(baseLoc,"/../data-processed/",type,"/",type,type_sub,interpolation,adjustment,adjustment_sub,"_event_",intr,".RData"))
    }
    
  }
  
}
adisch87/COMCORDEcml documentation built on Dec. 22, 2017, 11:11 a.m.