footprint_rmse_impact_test.r

#!/usr/bin/env Rscript

library(raster)
library(tidyverse)

# Time integrate footprints to a single layer
r0 <- sum(brick('footprint0.nc'))
r1 <- sum(brick('footprint1.nc'))

control <- values(r0)
interpolated <- values(r1)

r_squared <- cor(control, interpolated) ^ 2
message('r squared: ', signif(r_squared, 4))

# Define error as i - c
e <- interpolated - control

rmse <- function(x) sqrt(mean(x^2))

# Calculate rmse using all cells
base_rmse <- rmse(e)

# Estimate rmse excluding each individual cell
rm_rmse <- array(0, length(e))
for (i in 1:length(e)) {
  rm_rmse[i] <- rmse(e[-i])
}

# Aggregate footprint magnitude and rmse after cell exclusion
df <- data.frame(foot=control, rm_rmse=rm_rmse) %>%
  arrange(desc(foot)) %>%
  mutate(idx = 1:n(),
         drmse = rm_rmse - base_rmse)