functions.R

# Giacomo Roversi 2019

# functions to manipulate time strings
s2p <- function(str, dtformat="%Y%m%d%H%M"){return(as.POSIXct(str, format=dtformat, tz='GMT'))}
p2s <- function(datetime, dtformat="%Y%m%d%H%M"){return(as.character(datetime, format=dtformat))}
n2s <- function(unixtime, dtformat="%Y%m%d%H%M"){return(format(as.POSIXct(unixtime, origin = "1970-01-01", tz = 'GMT'), dtformat))}

# hourly accumulation on CMLs
accu1hr <- function(CmlRainfall, TIMESTEP){
  frac <- 60/TIMESTEP
  ids  <- unique(CmlRainfall$ID)
  dts  <- sort(unique(CmlRainfall$DateTime))
  
  totdts <- length(dts)
  totids <- length(ids)
  
  CmlHourlyData <- data.frame()
  
  
  startTime <- proc.time()
  
  for(selid in ids){                       # Ciclo sugli ID
    bool1 <- CmlRainfall$ID == selid
    numids <- which(ids == selid)
    
    elapsed_start <- round((proc.time()-startTime)[3],digits=1)
    remaining_start <- elapsed_start/numids*(totids-numids)
    
    lapTime <- proc.time()
    hourlydepth   <- 0
    count         <- 0
    
    for(seldt in dts){                         # Ciclo sui TIMESTEP
      bool2 <- CmlRainfall$DateTime == seldt
      numdts <- which(dts == seldt)
      # print(selid, quote = F)             # [debug]
      
      j <- which(bool1 & bool2)
      stopifnot(length(j) <= 1)              # Verifico univocità
      
      depth <- CmlRainfall$RainfallDepthPath[j]
      if(length(j) == 0 ){
        depth <- NA
      }
      
      hourlydepth <- hourlydepth + depth
      count       <- count + 1
      
      if(substr(seldt,11,12)=='00'          # - Cerco le ore intere
         & numdts > frac - 1                # - Attendo che sia passata almeno la prima ora
         & count == frac                    # - Ho tutte le 4 misure dell'ora
         & length(j) != 0){                 # - L'ora esiste per il link selezionato
        # print(seldt, quote = F)             # [debug]
        CmlHourlyData_row <- cbind(CmlRainfall[j,
                                               c("Frequency", "PathLength",  "XStart", 
                                                 "YStart", "XEnd", "YEnd", "Label", 
                                                 "Polarization", "Direction",  "ID")],
                                   DateTime = seldt,
                                   HourlyRainfallDepth = hourlydepth)
        
        CmlHourlyData <- rbind(CmlHourlyData, CmlHourlyData_row)
        
        hourlydepth <- 0
        count       <- 0
      }
      
      if(count > frac){        # in caso di errore (ad esempio per ora piena mancante)
        hourlydepth <- depth   # riparto con la somma delle accumulate
        count       <- 1       # come se fosse il primo intervallo
      } 
      # Remaining time estimation
      elapsed_lap <- round((proc.time()-lapTime)[3],digits=1)
      remaining_lap <- elapsed_lap/numdts*(totdts-numdts)
      cat(sprintf("\r%d:%d\t[%.1f]\t%d:%d\t[%.1f]", 
                  numdts, 
                  totdts, 
                  remaining_lap,
                  numids, 
                  totids, 
                  remaining_start))
    }
  }
  cat(sprintf("\nElapsed %.1f s", elapsed_start))
  return(CmlHourlyData)
}# VERY SLOW



# hourly accumulation on 15 min CML data
fast50x_accu1hr <- function(CmlRainfall){
  require(zoo)
  startTime <- proc.time()
  
  CmlRainfall <- CmlRainfall[order(CmlRainfall$ID, CmlRainfall$DateTime),]
  
  log1 <- substr(CmlRainfall$DateTime,11,12)=='00'
  sign <- c(T,F,F,F,T)
  
  wholehour <- rollapply(log1, width = 5, by = 1,
                         FUN = function(x) identical(x,y=sign),
                         align="right", 
                         fill = NA)
  
  samelink <- rollapply(CmlRainfall$ID, width = 4, by = 1,
                        FUN = function(x) length(unique(x)) == 1,
                        align="right", 
                        fill = NA)

  sel <- wholehour & samelink
  
  
  fouravg <- rollapplyr(CmlRainfall$RainfallMeanInt, width=4, by=1,
                        FUN = mean,
                        align="right", 
                        fill = NA)
  
  CmlHourlyData <- cbind(CmlRainfall[sel,
                                     c("DateTime","Frequency", "PathLength",  "XStart", 
                                       "YStart", "XEnd", "YEnd", "Label", 
                                       "Polarization", "Direction",  "ID")],
                         HourlyRainfallDepth = fouravg[sel])
  
  elapsed_start <- round((proc.time()-startTime)[3],digits=1)
  cat(sprintf("\nElapsed %.1f s", elapsed_start))
  
  
  filter <- which(is.na(CmlHourlyData$Frequency) 
                  & is.na(CmlHourlyData$DateTime)
                  & is.na(CmlHourlyData$PathLength)
                  & is.na(CmlHourlyData$XStart)
                  & is.na(CmlHourlyData$YStart)
                  & is.na(CmlHourlyData$XEnd)
                  & is.na(CmlHourlyData$YEnd)
                  & is.na(CmlHourlyData$Label)
                  & is.na(CmlHourlyData$Direction)
                  & is.na(CmlHourlyData$ID))
  
  require(txtplot)
  txtboxplot(CmlHourlyData$HourlyRainfallDepth[CmlHourlyData$HourlyRainfallDepth > 0.1])
  
  return(CmlHourlyData[-filter,])
}

# daily accumulation on 15 min CML data
fast50x_accu24hr <- function(CmlRainfall){
  require(zoo)
  startTime <- proc.time()
  
  CmlRainfall <- CmlRainfall[order(CmlRainfall$ID, CmlRainfall$DateTime),]
  
  log1 <- substr(CmlRainfall$DateTime,9,12)=='0000'
  sign <- c(T, rep(F, 95) ,T) # 24*4+1
  
  wholeday <- rollapply(log1, width = 97, by = 1,
                         FUN = function(x) identical(x,y=sign),
                         align="right", 
                         fill = NA)
  
  samelink <- rollapply(CmlRainfall$ID, width = 96, by = 1,
                        FUN = function(x) length(unique(x)) == 1,
                        align="right", 
                        fill = NA)
  
  sel <- wholeday & samelink
  
  
  fouravg <- rollapplyr(CmlRainfall$RainfallMeanInt, width=96, by=1,
                        FUN = mean,
                        align="right", 
                        fill = NA)
  
  CmlDailyData <- cbind(CmlRainfall[sel,
                                     c("DateTime","Frequency", "PathLength",  "XStart", 
                                       "YStart", "XEnd", "YEnd", "Label", 
                                       "Polarization", "Direction",  "ID")],
                         DailyRainfallDepth = fouravg[sel])
  
  elapsed_start <- round((proc.time()-startTime)[3],digits=1)
  cat(sprintf("\nElapsed %.1f s", elapsed_start))
  
  
  filter <- which(is.na(CmlDailyData$Frequency) 
                  & is.na(CmlDailyData$DateTime)
                  & is.na(CmlDailyData$PathLength)
                  & is.na(CmlDailyData$XStart)
                  & is.na(CmlDailyData$YStart)
                  & is.na(CmlDailyData$XEnd)
                  & is.na(CmlDailyData$YEnd)
                  & is.na(CmlDailyData$Label)
                  & is.na(CmlDailyData$Direction)
                  & is.na(CmlDailyData$ID))
  
  require(txtplot)
  txtboxplot(CmlDailyData$DailyRainfallDepth[CmlDailyData$DailyRainfallDepth > 0.1])
  
  return(CmlDailyData[-filter,])
}

# daily accumulation on hourly rasters
rast50x_accu24hr <- function(hrast){
  require(zoo)
  startTime <- proc.time()
  
  # CmlRainfall <- CmlRainfall[order(CmlRainfall$ID, CmlRainfall$DateTime),]
  
  log1 <- substr(names(hrast),10,13)=='0000'
  sign <- c(T, rep(F, 23) ,T) # 24+1
  
  wholeday <- rollapply(log1, width = 25, by = 1,
                        FUN = function(x) identical(x,y=sign),
                        align="right", 
                        fill = NA)
  
  samelink <- rollapply(CmlRainfall$ID, width =24, by = 1,
                        FUN = function(x) length(unique(x)) == 1,
                        align="right", 
                        fill = NA)
  
  sel <- wholeday & samelink
  
  
  fouravg <- rollapplyr(CmlRainfall$RainfallMeanInt, width=96, by=1,
                        FUN = mean,
                        align="right", 
                        fill = NA)
  
  CmlDailyData <- cbind(CmlRainfall[sel,
                                    c("DateTime","Frequency", "PathLength",  "XStart", 
                                      "YStart", "XEnd", "YEnd", "Label", 
                                      "Polarization", "Direction",  "ID")],
                        DailyRainfallDepth = fouravg[sel])
  
  elapsed_start <- round((proc.time()-startTime)[3],digits=1)
  cat(sprintf("\nElapsed %.1f s", elapsed_start))
  
  
  filter <- which(is.na(CmlDailyData$Frequency) 
                  & is.na(CmlDailyData$DateTime)
                  & is.na(CmlDailyData$PathLength)
                  & is.na(CmlDailyData$XStart)
                  & is.na(CmlDailyData$YStart)
                  & is.na(CmlDailyData$XEnd)
                  & is.na(CmlDailyData$YEnd)
                  & is.na(CmlDailyData$Label)
                  & is.na(CmlDailyData$Direction)
                  & is.na(CmlDailyData$ID))
  
  require(txtplot)
  txtboxplot(CmlDailyData$DailyRainfallDepth[CmlDailyData$DailyRainfallDepth > 0.1])
  
  return(CmlDailyData[-filter,])
}

################################################################################################

## POLYGONS GRIDS GENERATOR FUNCTION
# Creates polygons for printing function from the interpolation grid.
# Giacomo Roversi, 13 sept 2016 
# Corrected 29 oct 2017 
# Revised 31 jul 2019

PolyGridGen <- function(IntpGrid, SaveToFile, FileName){
  cat(sprintf('Loading...\n'))
  require(sp)
  require(rgdal)
  
  ok   <- 0
  jump <- 0
  
  while (ok != 1){
    width <- 1
    print(jump)
    
    # look for X-width of the grid
    while(IntpGrid[(jump+width+1),1] != IntpGrid[(jump+1),1]){
      width <- width+1
    }
    
    # search for a width larger than 10, otherwise skip to the next line
    if (width < 5){
      jump <- jump + width
      next()
    }else{
      ok <- 1
    }
    
    cat(sprintf('First acceptable line width: %d points\n', width)) 
  }
  step <- abs(IntpGrid[(jump+2),] - IntpGrid[(jump+width+1),]) # vector of the grid-side increments
  halfstep <- step/2						                               # compute the semi-diagonal
  polygons <- data.frame(X=numeric(), Y=numeric())
  totalpoints <- length(IntpGrid[,1])
  #cat(sprintf('IntpGrid dimensions: \n lon %d ° \n lat %d °\n', step[1], step[2]))
  print(step)
  
  for(i in 1:totalpoints){
    center <- IntpGrid[i,]
    r      <- 6*(i-1)
    
    polygons[r+1,] <- center + halfstep*c(-1,1)
    polygons[r+2,] <- center + halfstep*c(1,1)
    polygons[r+3,] <- center + halfstep*c(1,-1)
    polygons[r+4,] <- center + halfstep*c(-1,-1)
    polygons[r+5,] <- center + halfstep*c(-1,1)
    polygons[r+6,] <- NA
    
    cat(sprintf('Generating %d of %d polygons\r', i, totalpoints))
  }
  cat(sprintf('Generating %d of %d polygons\n', i, totalpoints))
  stopifnot(length(polygons[,1])==(6*totalpoints))
  
  if(SaveToFile){
    write.table(polygons, file=FileName, row.names=F)
  }
  
  return(polygons)
}



################################################################################################

# NLA Alternative function
# Giacomo Roversi, August 24, 2021
# Thakns to Elia Covi (NLA suspected wrong definition)

revWetDryNearbyLinkApMinMaxRSL <- function(Data,CoorSystemInputData=NULL,MinHoursPmin=6,PeriodHoursPmin=24,
                                        Radius=15,Step8=TRUE,ThresholdMedian=-1.4,ThresholdMedianL=-0.7,ThresholdNumberLinks=3,ThresholdWetDry=2)
{
  # Determine the middle of the area over which there are data 
  # (for reprojection onto a Cartesian coordinate system)
  if (!is.null(CoorSystemInputData))
  {
    Coor <- data.frame(x = c(min(Data$XStart, Data$XEnd), max(Data$XStart, Data$XEnd)), 
                       y = c(min(Data$YStart, Data$YEnd), max(Data$YStart, Data$YEnd)))
    coordinates(Coor) <- c("x", "y")
    proj4string(Coor) <- CRS(CoorSystemInputData) 
    CRS.latlon <- CRS("+proj=longlat +ellps=WGS84")
    Coor.latlon <- spTransform(Coor, CRS.lotlon)
    XMiddle <- (Coor.latlon$x[1] + Coor.latlon$x[2]) / 2
    YMiddle <- (Coor.latlon$y[1] + Coor.latlon$y[2]) / 2
  } else {
    XMiddle <- (min(Data$XStart, Data$XEnd) + max(Data$XStart, Data$XEnd)) / 2
    YMiddle <- (min(Data$YStart, Data$YEnd) + max(Data$YStart, Data$YEnd)) / 2
    CoorSystemInputData <- "+proj=longlat +ellps=WGS84"
  }
  
  # Set projection string
  projstring <- paste("+proj=aeqd +a=6378.137 +b=6356.752 +R_A +lat_0=",YMiddle,
                      " +lon_0=",XMiddle," +x_0=0 +y_0=0",sep="")
  
  # Set link IDs and time intervals
  Data$ID <- as.character(Data$ID)
  IDLink <- unique(Data$ID)
  N_links <- length(IDLink)
  t <- sort(unique(Data$DateTime))
  N_t <- length(t)
  
  # Make numeric representation of time in seconds from an arbitrary origin
  t_sec <- as.numeric(as.POSIXct(as.character(t), format = "%Y%m%d%H%M"))
  
  # Determine time interval length (in seconds)
  dt <- min(diff(t_sec))
  
  #Determine time indices for each entry
  t_ind <- rep(NA, length(Data$DateTime))
  for (i in 1 : N_t)
  {
    ind <- which(Data$DateTime == t[i])
    t_ind[ind] <- i
  }
  
  # Initialize arrays and vectors
  PminLink <- array(NA, c(N_t, N_links))
  array_ind <- array(NA, c(N_t, N_links))
  XStartLink <- rep(NA, N_links)
  YStartLink <- rep(NA, N_links)
  XEndLink <- rep(NA, N_links)
  YEndLink <- rep(NA, N_links)
  LengthLink <- rep(NA, N_links)
  
  # Loop over all links for coordinate transformation and putting data in an array
  for (p in 1 : N_links)
  {
    # Find indices corresppnding to this link
    Cond <- which(Data$ID == IDLink[p])
    
    #Convert coordinates to a system in km, centered on the area covered by the links
    Coor <- data.frame(x = c(Data$XStart[Cond[1]], Data$XEnd[Cond[1]]), 
                       y = c(Data$YStart[Cond[1]], Data$YEnd[Cond[1]]))
    coordinates(Coor) <- c("x", "y")
    proj4string(Coor) <- CRS(CoorSystemInputData) 
    CRS.cart <- CRS(projstring)
    Coor.cart <- spTransform(Coor, CRS.cart)
    XStartLink[p] <- Coor.cart$x[1]  # Easting (in km)
    YStartLink[p] <- Coor.cart$y[1]  # Northing (in km)
    XEndLink[p] <- Coor.cart$x[2]  # Easting (in km)
    YEndLink[p] <- Coor.cart$y[2]  # Northing (in km)
    
    LengthLink[p] <- Data$PathLength[Cond[1]] 
    
    # Store data from the considered link in an array
    PminLink[t_ind[Cond],p] <- Data$Pmin[Cond]
    array_ind[t_ind[Cond], p] <- Cond
  }
  
  # Initialize arrays
  PminLink_max <- array(NA, c(N_t, N_links))
  DeltaP <- array(NA, c(N_t, N_links))
  DeltaPL <- array(NA, c(N_t, N_links))
  ind_PrevPeriod <- rep(1, N_t)
  for (i in 2 : N_t)
  {
    # Determine index of time at most PeriodHoursPmin before current time interval
    int.ind = which(t_sec[ind_PrevPeriod[i - 1] : (i - 1)] > (t_sec[i] - PeriodHoursPmin * 3600))
    if (length(int.ind) > 0) {
      ind_PrevPeriod[i] <- min(int.ind) + ind_PrevPeriod[i - 1] - 1
      
      # Compute the time for which valid data are available, and check if this is sufficient
      t_valid <- colSums(!is.na(PminLink[ind_PrevPeriod[i] : i, ])) * dt
      links_valid <- which(t_valid >= (MinHoursPmin * 3600))
      if (length(links_valid) > 0)
      {
        for (j in links_valid)
        {
          # Compute maximum of Pmin over previous PeriodHoursPmin
          PminLink_max[i, j] <- max(PminLink[ind_PrevPeriod[i] : i, j], na.rm = TRUE)
        }
        
        # Compute Delta P and Delta P_L
        DeltaP[i, links_valid] <- PminLink[i, links_valid] - PminLink_max[i, links_valid]
        DeltaPL[i, links_valid] <- DeltaP[i, links_valid] / LengthLink[links_valid]
      }
    }
  }
  
  # Initialize dry and F vectors
  dry_vec <- rep(NA, length(Data$DateTime))
  F_vec <- rep(NA, length(Data$DateTime))
  for (i in 1 : N_links)
  {
    # Compute distances
    Distance1 <- sqrt( (XStartLink[i]-XStartLink)^2 + (YStartLink[i]-YStartLink)^2 )
    Distance2 <- sqrt( (XEndLink[i]-XStartLink)^2 + (YEndLink[i]-YStartLink)^2 ) 	
    Distance3 <- sqrt( (XStartLink[i]-XEndLink)^2 + (YStartLink[i]-YEndLink)^2 )
    Distance4 <- sqrt( (XEndLink[i]-XEndLink)^2 + (YEndLink[i]-YEndLink)^2 ) 
    # SelectDist <- which(Distance1 < Radius & Distance2 < Radius & Distance3 < Radius & 
    #                       Distance4 < Radius )
    
    # REVISED SELECTDIST
    SelectDist <- which( (Distance1 < Radius | Distance3 < Radius) & 
                           (Distance2 < Radius | Distance4 < Radius) )
    
    
    # Loop over all time intervals to compute medians and F values
    medianDeltaP <- rep(NA, N_t)
    medianDeltaPL <- rep(NA, N_t)
    for (j in 1 : N_t)
    {
      # Check if enough links are available for median and F computation
      if (sum(!is.na(DeltaP[j, SelectDist])) >= ThresholdNumberLinks)
      {
        medianDeltaP[j] = median(DeltaP[j, SelectDist], na.rm = TRUE)
        medianDeltaPL[j] = median(DeltaPL[j, SelectDist], na.rm = TRUE)
        F[j] <- sum(DeltaPL[ind_PrevPeriod[j] : j, i] - 
                      medianDeltaPL[ind_PrevPeriod[j] : j], na.rm = TRUE) * dt / 3600
      }
    }
    
    # Set dry indicator variable
    dry  <- rep(0, N_t)
    dry[medianDeltaP >= ThresholdMedian | medianDeltaPL >= ThresholdMedianL] <- 1
    dry[is.na(medianDeltaP) | is.na(medianDeltaPL)] <- NA
    
    # Perform step 8 if desired
    if (Step8)
    {
      ind_wet <- which(dry == 0 & DeltaP[, i] < (-1 * ThresholdWetDry))
      int_dry <- dry
      dry[ind_wet[ind_wet > 1] - 1] <- 0
      dry[ind_wet[ind_wet > 2] - 2] <- 0
      dry[ind_wet[ind_wet < length(dry)] + 1] <- 0
      dry[is.na(int_dry)] <- NA
    }
    
    # Map arrays of dry and F to vectors corresponding to input data frame
    dry_vec[array_ind[!is.na(array_ind[, i]), i]] <- dry[!is.na(array_ind[, i])]
    F_vec[array_ind[!is.na(array_ind[, i]), i]] <- F[!is.na(array_ind[, i])]
  }
  # Set return data frame
  return_value <- data.frame(Dry = dry_vec, F = F_vec)
  return(return_value)
  
}