6  Fronts (or any variable) by equal-size grid

6.1 Data import

Load the required packages.

library(terra)
library(sf)
library(exactextractr)
library(dplyr)
library(nngeo)
library(stringr)

6.2 Generic Function to replace NAs with nearest neighbor.

This is just a generic function that uses nngeo R package. Feel free to used and adapt it to your needs.

  fCheckNAs <- function(df, vari) {
    if (sum(is.na(pull(df, !!sym(vari))))>0){ # Check if there are NAs
        
        gp <- df %>%
          mutate(isna = is.finite(!!sym(vari))) %>%
          group_by(isna) %>%
          group_split()
        
        out_na <- gp[[1]] # DF with NAs
        out_finite <- gp[[2]] # DF without NAs
        
        d <- st_nn(out_na, out_finite) %>% # Get nearest neighbour
          unlist()
        
        out_na <- out_na %>%
          mutate(!!sym(vari) := pull(out_finite, !!sym(vari))[d])
        
        df <- rbind(out_finite, out_na)
        
      }
      return(df)
    }

6.3 Read, Extract and weighted mean interpolation

# Create the projection
  proj.geo = "+proj=robin +lon_0=0 +datum=WGS84 +units=m +no_defs"
# Reading equal grid file
  shp_file <- st_read("data/PUs_MZ_100km2.shp") %>% 
    st_transform(crs = terra::crs(proj.geo))
# Read raster object
  rs_file <- rast("data/BOAonMUR_SWIO_Y2003-M1-D1.nc")
  crs(rs_file) <- terra::crs("+proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0")
  weight_rs <- terra::cellSize(rs_file)
  rs_file <- terra::project(rs_file, y = terra::crs(proj.geo), method = "near")
  weight_rs <- terra::project(weight_rs, y = terra::crs(proj.geo), method = "near")
# Getting value by polygon
  rs_bypu <- exact_extract(rs_file, 
                           shp_file, 
                           "weighted_mean", 
                           weights = weight_rs, 
                           append_cols = TRUE, 
                           full_colnames = TRUE)
  rs_shp <- dplyr::right_join(shp_file, rs_bypu, "FID")
  colnames(rs_shp) <- c(stringr::str_remove_all(string = names(rs_shp), pattern = "weighted_mean."))

6.4 RUN: replace NAs with nearest neighbor.

nms <- names(rs_shp)
nms <- nms[nms != "geometry" & nms != "FID"]

single <- rs_shp %>% 
  dplyr::select(FID, nms[1])
rs_sfInt <- fCheckNAs(df = single, vari = names(single)[2]) %>% 
  as_tibble() %>%
  dplyr::arrange(FID) %>%
  dplyr::select(-FID, -geometry, -isna)

saveRDS(rs_sfInt, "data/BOAonMUR_SWIO_Y2003-M1-D1.rds")

6.5 Plot the output

pus <- st_read("data/PUs_MZ_100km2.shp") %>%
  st_transform(crs = robin)
sf1 <- readRDS("data/BOAonMUR_SWIO_Y2003-M1-D1.rds")
df1 <- cbind(pus, sf1) %>% 
  st_transform(crs = robin)

p1 <- ggplot() +
  geom_sf(data = df1, aes(fill = temp_gradient.area), colour = NA) +
  geom_sf(data = world_sfRob, size = 0.05, fill = "grey20") +
  theme_bw() +
  theme(legend.title = element_text(angle = 0, size = rel(0.7)),
        plot.title = element_text(face = "plain", size = 22, hjust = 0.5),
        axis.text.x = element_text(size = rel(1), angle = 0),
        axis.text.y = element_text(size = rel(1), angle = 0),
        axis.title = element_blank()) +
  scale_fill_distiller(palette = "Spectral",
                       limits = c(min(df1$temp_gradient.area), max(df1$temp_gradient.area)),
                       direction = -1,
                       oob = scales::squish,
                       guide = guide_colourbar(title.position = "top", title = "thermal gradient fronts"))