test_codes/spatial_qc.R

In jieyu97/QCwind: Quality control procedures for citizen science wind observations

#' @name spatial_qc
#' @title Spatial quality control of input observations
#' @description Implementing the spatial quality control, and receive the information
#' of which observations faild the quality control.
#' @param data a data.frame of observations by different stations, each
#' column represents a station, and the rows represent standard recording datetimes.
#' @param ref.data a data.frame that has the same settings with `data`; the stations
#' in `ref.data` are official stations used for reference.
#' @param data.nbhd a list, the selected neighboring stations of each station in `data`
#' (obtained by the function `neighbors_selection`).
#' @param data.order a vector that indicating the optimal order of stations in `data`
#' to do spatial quality control (obtained by the function `neighbors_selection`).
#' @param fail.flag a character/string that represents the name of flag where an
#' observation fails the spatial quality control.
#' @param method a character indicating the spatial quality control method to be used;
#' `idw`: the inverse distance weighting method, this is also the default;
#' `srt`: the spatial weighted regression test (! needed much longer processing time).
#' @return a list of four data.frames that have the same size with `data`:
#' `estimation`: the estimation of each target observation generated by neighboring reference
#' observations;
#' `nbhd_sd`: the standard deviations of neighboring reference observations related to
#' each target observation;
#' `flag`: the flagging information that indicates whether a target observation passes
#' spatial quality control of not, `P` is a pass, `fail.flag` is a fail, `isolated` represent
#' isolated observations;
#' `new_data`: the new data of the input `data` that removes the failed observations.
#' @examples
#' wow_spatial = spatial_qc(wow_ws_test, knmi_ws_test, wow_nbhd, wow_order, fail.flag = 'F')

spatial_qc = function(data, ref.data, data.nbhd, data.order, fail.flag, method = 'idw')
{
  # set empty matrix
  estimation_idw_emd_mat = matrix(NA, nrow = 157825, ncol = 39) # estimation 1
  estimation_idw_rmse_mat = estimation_idw_emd_mat # estimation 2
  sd_idw_emd_mat = estimation_idw_emd_mat # nbhd sd 1
  sd_idw_rmse_mat = estimation_idw_emd_mat # nbhd sd 2

  spatial_emd_ws = spatial_ws
  spatial_emd_flag = spatial_ws
  spatial_rmse_ws = spatial_ws
  spatial_rmse_flag = spatial_ws


  estimation_idw_emd_mat_bad = estimation_idw_emd_mat
  sd_idw_emd_mat_bad = sd_idw_emd_mat
  spatial_emd_ws_bad = spatial_emd_ws
  spatial_emd_flag_bad = spatial_emd_flag
  ########################3 huge loop over WOW stations
  # EMD
  for (k in 1:38 ) {
    i = wow_station_order[k]
    # observation WOW
    observation = spatial_emd_ws[,i] # spatial_emd_ws spatial_ws

    estimation_idw_emd = as.numeric(matrix(NA, nrow = 1, ncol = 157825)) # estimation 1
    sd_idw_emd = estimation_idw_emd # nbhd sd 1

    nbhd_label = wow_spatial_nbhd[[k]]
    weight_emd = ( emd_radius^2 - ( emd_matrix_ws[i,nbhd_label] )^2 ) / # dist_rmse_ws emd_matrix_ws
      ( emd_radius^2 + ( emd_matrix_ws[i,nbhd_label] )^2 ) # rmse_radius emd_radius

    for (t in which(!is.na(observation)) ) {
      # get nbhd wind speed & weight
      nbhd_infomation = cbind(nbhd_label, weight_emd, 1)
      nbhd_infomation[,3] = as.numeric( spatial_emd_ws[t,nbhd_label] ) # spatial_emd_ws spatial_ws

      if ( sum(!is.na(nbhd_infomation[,3])) < 3 ) {
        estimation_idw_emd[t] = NA
        sd_idw_emd[t] = NA
      } else {
        estimation_idw_emd[t] = sum( nbhd_infomation[,3] * nbhd_infomation[,2], na.rm = TRUE ) /
          sum( nbhd_infomation[!is.na(nbhd_infomation[,3]),2] )
        sd_idw_emd[t] = sd( nbhd_infomation[,3], na.rm = TRUE )
      }
    }

    flag_after_emd = ifelse( observation > (estimation_idw_emd + 2*sd_idw_emd) |
                               observation < (estimation_idw_emd - 2*sd_idw_emd),
                             1, 0 )
    obs_after_emd = ifelse( observation > (estimation_idw_emd + 2*sd_idw_emd) |
                              observation < (estimation_idw_emd - 2*sd_idw_emd),
                            NA, observation )

    estimation_idw_emd_mat[,i] = estimation_idw_emd
    sd_idw_emd_mat[,i] = sd_idw_emd
    spatial_emd_ws[,i] = obs_after_emd
    spatial_emd_flag[,i] = flag_after_emd
    print(k)
  }
}