R/ad_hoc.R

In carlosfperuchena/satellite2ground: Local adaptation of modeled solar irradiance series with coincident ground measurements

#'  Local adaptation of solar irradiance series with coincident ground measurements
#'
#'
#' @param latitude Sun elevation angle (degrees) above horizon
#' @param z Site elevation above sea level (m)
#' @param calibrating_period Data frame containing, at least: Time (UTC time stamp, in format "YYYY-MM-DD hh:mm:ss"); GHI_measured (measured GHI values); GHI_model (modeled GHI values); Elevation (Sun elevation above horizon, in degrees); Kc (Clear-sky index (ratio of GHI_model with the concomitant modeled GHI under clear-sky conditions); DNI_measured (measured DNI values, optional parameter); DNI_model (modeled DNI values, optional parameter)
#'
#' @param target_period  Data frame containing, at least: Time (time stamp, in format "YYYY-MM-DD hh:mm:ss"); GHI_model (modeled GHI values); Elevation (Sun elevation above horizon, in degrees); Kc (Clear-sky index (ratio of GHI_model with the concomitant modeled GHI under clear-sky conditions); DNI_model (modeled DNI values, optional parameter)
#'
#' @return Time: UTC Time stamp, in format "YYYY-MM-DD hh:mm:ss"
#' @return GHI_adapted: locally adapted GHI series
#' @return DNI_adapted: locally adapted DNI series (if optional parameters have been provided)
#' @export
#' @import glmulti
#' @import solaR

ad_hoc <- function(latitude,calibrating_period,target_period, z){

  #--------------------------------------------------------------------------------------------------------------------------
  # Optical air mass
  calibrating_period$air_mass_kasten = m_Kasten(calibrating_period$Elevation,z )
  target_period$air_mass_kasten = m_Kasten(target_period$Elevation,z)

  #--------------------------------------------------------------------------------------------------------------------------
  # Top of atmosphere solar irradiance
  calibrating_period$Top_irradiance = TOA_h(latitude,calibrating_period$Elevation,calibrating_period$Time)
  target_period$Top_irradiance = TOA_h(latitude,target_period$Elevation,target_period$Time)

  #--------------------------------------------------------------------------------------------------------------------------
  # Non-dimensional parameters
      # Calibrating period
      calibrating_period$KT_model = calibrating_period$GHI_model / calibrating_period$Top_irradiance
      calibrating_period$KT_measured = calibrating_period$GHI_measured / calibrating_period$Top_irradiance
      calibrating_period$KT_measured_modified = calibrating_period$KT_measured / (1.031*exp( -1.4/(0.9 + 9.4/calibrating_period$air_mass_kasten)      )   + 0.1)
      calibrating_period$kd_measured = (calibrating_period$GHI_measured - calibrating_period$DNI_measured * sin(calibrating_period$Elevation*pi/180))/calibrating_period$GHI_measured
      calibrating_period$kd_model = (calibrating_period$GHI_model - calibrating_period$DNI_model * sin(calibrating_period$Elevation*pi/180))/calibrating_period$GHI_model

      # Target period
      target_period$KT_model = target_period$GHI_model / target_period$Top_irradiance
      target_period$kd_model = (target_period$GHI_model - target_period$DNI_model * sin(target_period$Elevation*pi/180))/target_period$GHI_model

  #--------------------------------------------------------------------------------------------------------------------------
  # The local adaptation only applies to solar elevations above 1 degree:
  target_period_low_elev = target_period[which(target_period$Elevation <= 1),]
  calibrating_period_high_elev = calibrating_period[which(calibrating_period$Elevation > 1),]
  target_period_high_elev = target_period[which(target_period$Elevation > 1),]



  #--------------------------------------------------------------------------------------------------------------------------
  # Local adaptation of GHI
  GHI_adaptation_proc <-  glm(KT_measured ~ KT_model + air_mass_kasten + Kc + Elevation,data=calibrating_period_high_elev)
  Best_GHI_adaptation <- glmulti(GHI_adaptation_proc,level = 2,crit="AIC", family = gaussian(link = "identity"), confsetsize=100 ,plotty = F, report = F)
  model_predicting_kt <- Best_GHI_adaptation@objects[[1]];vars <- names(coef(model_predicting_kt))[-1]
  kt_Locally_adapted = as.numeric(predict(model_predicting_kt, target_period_high_elev, family = gaussian(link = "identity")))
  GHI_Locally_adapted = kt_Locally_adapted*target_period_high_elev$Top_irradiance



  #--------------------------------------------------------------------------------------------------------------------------
  #--------------------------------------------------------------------------------------------------------------------------
  # Required inputs for local adaptation of DNI ARE NOT available
  if(min(length(calibrating_period$DNI_measured), length(target_period$DNI_model)) == 0){

    # Prepare data frames
    df_daytime <- data.frame(
      Time = target_period_high_elev$Time,
      GHI_adapted = as.numeric(GHI_Locally_adapted))

    df_low_elev<- data.frame(
      Time = target_period_low_elev$Time,
      GHI_adapted = target_period_low_elev$GHI_model)

    # Merge of daytime (sun elevation above 1 degree) and low elevation (sun elevation equal or below 1 degree) series
    sal = df_daytime
    if(length(df_low_elev$Time) > 0){sal = rbind(df_daytime,df_low_elev)}
    sal = sal[order(year(sal$Time)*1000 + doy(sal$Time) + (round(hour(sal$Time)/24, 2)) ),]

  }
  #--------------------------------------------------------------------------------------------------------------------------
  #--------------------------------------------------------------------------------------------------------------------------


  #--------------------------------------------------------------------------------------------------------------------------
  #--------------------------------------------------------------------------------------------------------------------------
  # Required inputs for local adaptation of DNI ARE available
  if(min(length(calibrating_period$DNI_measured), length(target_period$DNI_model)) > 0){

  adapted_kt_for_calibrating_period <- predict(model_predicting_kt, calibrating_period_high_elev, family = gaussian(link = "identity"))
  calibrating_period_high_elev$non_dim_kt = as.numeric(adapted_kt_for_calibrating_period)/(1.031*exp( -1.4/(0.9 + 9.4/calibrating_period_high_elev$air_mass_kasten))+ 0.1)
  target_period_high_elev$non_dim_kt = kt_Locally_adapted/(1.031*exp( -1.4/(0.9 + 9.4/target_period_high_elev$air_mass_kasten))+ 0.1)
  DNI_adaptation_proc <-  glm(kd_measured ~ poly(non_dim_kt, 4) +  Elevation + air_mass_kasten + Kc, data=calibrating_period_high_elev)

  Best_DNI_adaptation <- glmulti(DNI_adaptation_proc,level = 2, crit="AIC", family =  "gaussian", confsetsize=100)
  model <- Best_DNI_adaptation@objects[[1]];vars <- names(coef(model))[-1]
  salida_kd <- predict(model, target_period_high_elev)
  DHI_Fit = as.numeric(salida_kd)*GHI_Locally_adapted
  DNI_Fit = (GHI_Locally_adapted - DHI_Fit) / sin(target_period_high_elev$Elevation*pi/180)

  df_daytime <- data.frame(
    Time = target_period_high_elev$Time,
    GHI_adapted = as.numeric(GHI_Locally_adapted),
    DNI_adapted = as.numeric(DNI_Fit))

  # No valid adapted values:
  # - negative DNI values
  df_daytime[which(df_daytime$DNI_adapted < 0), c(2:3)] =
    target_period_high_elev[which(df_daytime$DNI_adapted < 0), c(2:3)]

  # - Extremely high DNI values
  df_daytime[which(df_daytime$DNI_adapted > 1.05*max(calibrating_period$DNI_measured)),c(2:3)] =
    target_period_high_elev[which(df_daytime$DNI_adapted > 1.05*max(calibrating_period$DNI_measured)), c(2:3)]

  df_low_elev<- data.frame(
    Time = target_period_low_elev$Time,
    GHI_adapted = target_period_low_elev$GHI_model,
    DNI_adapted = target_period_low_elev$DNI_model)

  sal = df_daytime
  # Merge of daytime (sun elevation above 10 degrees) and low elevation (sun elevation equal or below 10 degrees) series
  if(length(df_low_elev$Time) > 0){sal = rbind(df_daytime,df_low_elev)}
  sal = sal[order(year(sal$Time)*1000 + doy(sal$Time) + (round(hour(sal$Time)/24, 2)) ),]
  }

  return(sal)

}