AIM : Compare interpolation outputs from 3 methods: (1) Moving window (MW) following the coastline and using the relationships of Kd with Distance to shore (2) Nearest-Neighbor (NN) (3) Inverse-Distance Weighting (IDW).

Table of content

  1. Moving Window (MW)

  2. Nearest-Neighbor(NN)

  1. Inverse-Distance Weighting (IDW)

  2. Comparison

  1. Semi-Variogram exploration

  2. Bibliography

R Packages

#knitr::opts_chunk$set(echo = TRUE)
library(sp)
library(raster)
library(tidyverse)
## -- Attaching packages --------------------------------------- tidyverse 1.3.1 --
## v ggplot2 3.3.5     v purrr   0.3.4
## v tibble  3.1.2     v dplyr   1.0.7
## v tidyr   1.1.3     v stringr 1.4.0
## v readr   1.4.0     v forcats 0.5.1
## -- Conflicts ------------------------------------------ tidyverse_conflicts() --
## x tidyr::extract() masks raster::extract()
## x dplyr::filter()  masks stats::filter()
## x dplyr::lag()     masks stats::lag()
## x dplyr::select()  masks raster::select()
library(sf)
## Linking to GEOS 3.9.0, GDAL 3.2.1, PROJ 7.2.1
library(plotly)
## 
## Attaching package: 'plotly'
## The following object is masked from 'package:ggplot2':
## 
##     last_plot
## The following object is masked from 'package:raster':
## 
##     select
## The following object is masked from 'package:stats':
## 
##     filter
## The following object is masked from 'package:graphics':
## 
##     layout
library(viridis)
## Loading required package: viridisLite
library(tabularaster)
## Registered S3 method overwritten by 'spatstat.geom':
##   method     from
##   print.boxx cli
library(RColorBrewer)
library(gstat)
library(broom)
library(knitr)
library(stargazer)
## 
## Please cite as:
##  Hlavac, Marek (2018). stargazer: Well-Formatted Regression and Summary Statistics Tables.
##  R package version 5.2.2. https://CRAN.R-project.org/package=stargazer

Moving Window (MW) 

## Part 4: Moving the window, storing estimated values, producing full estimated map - Size of window consistent vs DShore (keep 15km of pixels)
#####
## Loading data - KPAR raster, crop in KK zone - Coastline shp file
#####
## Zone
xmin <- 1600000
xmax <- 1730000
ymin <- 5200000
ymax <- 5390000
##

## Load raster of KPAR + crop KK
kpar <- raster('A20021822002212_MO_KPAR_exp_coastal_v01.nc')

kpar_KK <- crop(kpar,extent(xmin, xmax, ymin, ymax))

kpar_KK_sp <- rasterToPoints(kpar_KK,spatial=T)
kpar_KK_sf <- st_as_sf(kpar_KK_sp) %>% rename(KPAR=DownwellingPARattenuation)
##

## Coastline + crop KK
coastline <- st_read('Coastline_Reproj_Line.shp') %>% 
  st_set_crs(.,crs(kpar_KK_sf))

coastline_KK <- st_crop(coastline,xmin=xmin,ymin=ymin,xmax=xmax,ymax=ymax)
##

## Convert Coastline Multinestring to Raster
coastline_KK_ras <- rasterize(coastline_KK,kpar_KK)
##

## Load DShore from previous Bathymetry grid
dShore_bathy <- raster('NZ_Dist2Shore.tif')

dShore_bathy_KK <- crop(dShore_bathy,extent(xmin, xmax, ymin, ymax))

dShore_KK_sp <- rasterToPoints(dShore_bathy_KK,spatial=T)
dShore_KK_sf <- st_as_sf(dShore_KK_sp) %>% rename(DShore=NZ_Dist2Shore) %>% 
  st_set_crs(.,crs(kpar_KK_sf))
##

## Moving the window, storing estimated values, producing full estimated map - Size of window consistent vs DShore (keep ~10km of pixels)
rot <- function(a){
  # Rotation matrix, 
  # Takes angle, a, in radian as input
  # cf https://r-spatial.github.io/sf/articles/sf3.html
  # and https://geocompr.robinlovelace.net/geometric-operations.html
  matrix(c(cos(a), sin(a), -sin(a), cos(a)), 2, 2) 
}
##

## Initialization of Box - Width, Length, Buff_dist
start_pix_ID <- 1
length <- 50#20 #Length factor of normal line
buff_dist <- 2500

dShore_max <- 5000#10000

kpar_KK_sp <- rasterToPoints(kpar_KK,spatial=T)
kpar_KK_sf <- st_as_sf(kpar_KK_sp) %>% 
  rename(KPAR=DownwellingPARattenuation) %>% 
  st_join(dShore_KK_sf,.,by=c('geometry'),left=T) %>% #Merge with DShore straight away here, keep NA in KPAR
  dplyr::filter(DShore<=dShore_max)    #+ filter by DShore<dShore_max
##

## Step 1: Get index (not x,y as before) of coastline pixels
#idx_coastline <- which(!values(coastline_KK_ras)==0)
idx_coastline <- which(!is.na(values(coastline_KK_ras)))
#end_domain <- nrow(idx_coastline)
end_domain <- length(idx_coastline)
#The box will go from start_pix_ID to end_domain
##


## Matrix of cells (x,y) which need to be filled (which have KPAR==NA)
kpar_coastal_NA <- xyFromCell(coastline_KK_ras,which(is.na(values(kpar_KK))==T)) %>% 
  as_tibble() %>% add_column(Pred = NA)

kpar_coastal_NA_ls <- list()
##

system.time(
  #for (i in start_pix_ID:end_domain){
  for (i in 20:500){
    print(i)
    
    ## Step 2: Select pixel #1 (or i), get adjacent pixels, only keep adjacent with not NA, get x/y
    idx_adj <- adjacent(coastline_KK_ras,idx_coastline[i],directions=8,include=T)#Issue of taking all ngb cells, regardless if it's coastline (non 0) or not
    idx_adj_coast <- idx_adj[which(!is.na(coastline_KK_ras)[idx_adj[,2]]),2]#Here only keep pixels without NA, i.e, coastline
    
    xy_box <- xyFromCell(coastline_KK_ras,idx_adj_coast)
    ##
    
    ## Step 3: Keep First and Last Points
    pixs_coastline_fl_box <- xy_box[c(1,nrow(xy_box)),]
    #--> Can be more thought of as first and last points not always best, cf plots 
    ##
    
    ## Step 4: Draw lines between first and last points in box
    pixs_coastline_line <- st_linestring(pixs_coastline_fl_box) %>% 
      st_sfc() %>% st_set_crs(crs(kpar_KK))
    ##
    
    ## Step 5: Create normal/perpendicular line by rotating by 90deg + extend by factor length +
    # cut normal line in 2 lines around centroid and keep offshore line
    cntrd <- st_centroid(pixs_coastline_line) %>% 
      st_sfc() %>% st_set_crs(crs(kpar_KK))
    
    pixs_coastline_normal <- (pixs_coastline_line - cntrd) * rot(pi/2) *length + cntrd
    pixs_coastline_normal <- pixs_coastline_normal %>% 
      st_sfc() %>% st_set_crs(crs(kpar_KK))
    
    pixs_coastline_normal_split <- st_split(pixs_coastline_normal,cntrd) %>% 
      st_collection_extract(.,"LINESTRING")
    
    # Let's continue by only keeping bit==1, meanig that what follows (cut, pixs_coastline_normal_fixLength) not needed
    pixs_coastline_normal <- pixs_coastline_normal_split[1]
    ##
    
    ## Step 6: Create buffer around normal/perpendicular line
    pixs_coastline_normal_buff <- st_buffer(pixs_coastline_normal,dist=buff_dist)
    ##
    
    ## Step 7: Crop KPAR pixels in buffered polygon
    kpar_dshore_buff_allpix <- st_intersection(kpar_KK_sf,pixs_coastline_normal_buff)
    ##
    if (sum(!is.na(kpar_dshore_buff_allpix$KPAR))>1){#nkpar_min
      ## Step 8: Estimate KPAR at missing pixels (where DShore is), shouldnt be issue of land pixels
      loess_out <- loess(KPAR ~ DShore, kpar_dshore_buff_allpix,control = loess.control(surface = "direct"))
      
      ## Step 9: Create list of predictions at each pixels 
      # --> Will be use to have several estimates per pixels (keep median at the end)
      kpar_dshore_buff_NA <- kpar_dshore_buff_allpix %>% 
        dplyr::filter(is.na(KPAR)==T) %>% 
        mutate(KPAR_pred=predict(loess_out,DShore)) %>% 
        mutate(x = sf::st_coordinates(.)[,1],y = sf::st_coordinates(.)[,2])
      
      kpar_coastal_NA_ls[[i]] <- kpar_dshore_buff_NA %>% dplyr::select(-KPAR)
      # THIS LIST KEEPS IN MEMORY THE X, Y and KPAR Prediction AT EACH PIXELS THAT HAVE KPAR==NA
      ##
      
    }else{
      kpar_coastal_NA_ls[[i]] <- NA
    }
  }
)#156s for 602i

## Unlist and unnest kpar_coastal_NA_ls, calculate median of KPAR_Pred to each pixels (x,y) 
# in new object containing pred at NA pixels, make sf and merge with kpar_dShore_KK_sf 
  
kpar_KK_pred_sf <- tibble(
  x = map(kpar_coastal_NA_ls, "x"),
  y = map(kpar_coastal_NA_ls, "y"),
  DShore = map(kpar_coastal_NA_ls, "DShore"),
  KPAR_pred = map(kpar_coastal_NA_ls, "KPAR_pred")) %>% #Good for switching from list to tibble which contains lists
  unnest(cols = c(x, y, DShore, KPAR_pred)) %>% 
  group_by(x, y,DShore) %>% 
  #summarise(KPAR_pred_list = list(KPAR_pred)) %>% #To check number of predictions per pixels
  summarise(KPAR = median(KPAR_pred)) %>% #Name of KPAR column should match name column in kpar_KK_sf for easy rbind
  st_as_sf(.,coords = c('x','y'), crs=crs(kpar_KK))
# --> Yeeeew


kpar_KK_sf_noNA <- kpar_KK_sf %>% 
  dplyr::filter(!is.na(KPAR))

kpar_KK_merge_sf <- rbind(kpar_KK_sf_noNA,kpar_KK_pred_sf)

st_write(kpar_KK_sf,'KPARKK_Before_Large_Win4_v2.shp')
st_write(kpar_KK_merge_sf,'KPARKK_After_Large_Win4_v2.shp') 
st_write(kpar_KK_pred_sf,'KPARKK_Estimated_Large_Win4_v2.shp') 

## Map of KPAR after/before smoothing

# --> To load kpar_KK_sf and kpar_KK_merge_sf here
kpar_KK_sf <- st_read('KPARKK_Before_Large_Win4_v2.shp')
## Reading layer `KPARKK_Before_Large_Win4_v2' from data source 
##   `C:\Users\thoralf\OneDrive - NIWA\Documents\Sat_Data\NETCDF_ORIGINALS\KPAR\MO\Gap_Fill_Trials\KPARKK_Before_Large_Win4_v2.shp' 
##   using driver `ESRI Shapefile'
## Simple feature collection with 4505 features and 2 fields
## Geometry type: POINT
## Dimension:     XY
## Bounding box:  xmin: 1600500 ymin: 5226500 xmax: 1711000 ymax: 5390000
## CRS:           unknown
kpar_KK_merge_sf <- st_read('KPARKK_After_Large_Win4_v2.shp') #%>% 
## Reading layer `KPARKK_After_Large_Win4_v2' from data source 
##   `C:\Users\thoralf\OneDrive - NIWA\Documents\Sat_Data\NETCDF_ORIGINALS\KPAR\MO\Gap_Fill_Trials\KPARKK_After_Large_Win4_v2.shp' 
##   using driver `ESRI Shapefile'
## Simple feature collection with 4338 features and 2 fields
## Geometry type: POINT
## Dimension:     XY
## Bounding box:  xmin: 1600500 ymin: 5226500 xmax: 1711000 ymax: 5390000
## CRS:           unknown
  #distinct()#Get rid of duplicated points, why?
kpar_KK_estimates_sf <- st_read('KPARKK_Estimated_Large_Win4_v2.shp')
## Reading layer `KPARKK_Estimated_Large_Win4_v2' from data source 
##   `C:\Users\thoralf\OneDrive - NIWA\Documents\Sat_Data\NETCDF_ORIGINALS\KPAR\MO\Gap_Fill_Trials\KPARKK_Estimated_Large_Win4_v2.shp' 
##   using driver `ESRI Shapefile'
## Simple feature collection with 1025 features and 2 fields
## Geometry type: POINT
## Dimension:     XY
## Bounding box:  xmin: 1625500 ymin: 5252000 xmax: 1708500 ymax: 5390000
## CRS:           unknown
## Zone
xmin <- 1600000
xmax <- 1730000
ymin <- 5200000
ymax <- 5390000
##

## Coastline + crop KK
coastline <- st_read('Coastline_Reproj_Line.shp') %>% 
  st_set_crs(.,crs(kpar_KK_sf))
## Reading layer `Coastline_Reproj_Line' from data source 
##   `C:\Users\thoralf\OneDrive - NIWA\Documents\Sat_Data\NETCDF_ORIGINALS\KPAR\MO\Gap_Fill_Trials\Coastline_Reproj_Line.shp' 
##   using driver `ESRI Shapefile'
## Simple feature collection with 322 features and 7 fields
## Geometry type: MULTILINESTRING
## Dimension:     XY
## Bounding box:  xmin: 1089971 ymin: 4748121 xmax: 2089558 ymax: 6221557
## CRS:           unknown
coastline_KK <- st_crop(coastline,xmin=xmin,ymin=ymin,xmax=xmax,ymax=ymax)
##

## Load DShore from previous Bathymetry grid
dShore_bathy <- raster('NZ_Dist2Shore.tif')

dShore_bathy_KK <- crop(dShore_bathy,extent(xmin, xmax, ymin, ymax))

# DSHore as tibble
dShore_KK_tb <- dShore_bathy_KK %>% 
  as_tibble(.,xy=T) %>% 
  drop_na() %>% 
  rename(DShore = cellvalue)
##

# Before
p1 <- ggplot() + 
  geom_sf(data=kpar_KK_sf,aes(color=KPAR)) +
  scale_color_viridis(direction=1, limits = c(0.04, .35),na.value='transparent') +
  geom_sf(data=coastline_KK)
#ggplotly(p)

# After
p2 <- ggplot() + 
  geom_sf(data=kpar_KK_merge_sf,aes(color=KPAR)) +
  scale_color_viridis(direction=1, limits = c(0.04, .35),na.value='transparent') +
  geom_sf(data=coastline_KK)
#ggplotly(p)
##

subplot(p1, p2, nrows = 1, margin = 0.04)

Figure 1 - KdPAR before/after MW interpolation.

Nearest-Neighbor (NN)

NN 5 

## Part 2: Perform NN interpolation
# cf https://rspatial.org/raster/analysis/4-interpolation.html#calfornia-air-pollution-data
#####
kpar_KK_ras <- raster('KPARKK_Large.tif')

# Convert KPAR raster to sp points
kpar_KK_sp <- rasterToPoints(kpar_KK_ras, fun=NULL, spatial = T)

#create grid from kpar raster
grid_r <- raster(kpar_KK_ras)

gs <- gstat(formula=KPARKK_Large ~1, locations=kpar_KK_sp, nmax=5, set=list(idp = 0)) 
#using 5 neighboors/points + equal weight idp=0 
nn_5 <- interpolate(grid_r, gs)

nnmsk_5 <- mask(nn_5, dShore_bathy_KK)# Use Dshore to mask nn

writeRaster(nnmsk_5,'KPARKK_NN5_mask_Large.tif','GTiff')
nnmsk_5 <- raster('KPARKK_NN5_mask_Large.tif')
kpar_KK_ras <- raster('KPARKK_Large.tif')

par(mfrow=c(1,2))
plot(kpar_KK_ras,col=viridis(100),zlim=c(0.08,.35))
plot(coastline_KK,col='black',add=T,legend = F)

plot(nnmsk_5,col=viridis(100),zlim=c(0.08,.35))
plot(coastline_KK,col='black',add=T,legend = F)
Figure 2 - KdPAR after NN interpolation (ngb 5).

Figure 2 - KdPAR after NN interpolation (ngb 5).

NN 9 

gs <- gstat(formula=DownwellingPARattenuation~1, locations=kpar_KK_sp, nmax=9, set=list(idp = 0)) 
#using 9 neighboors/points + equal weight idp=0 
nn <- interpolate(grid_r, gs)

nnmsk <- mask(nn, dShore_bathy_KK)# Use Dshore to mask nn

writeRaster(nnmsk,'KPARKK_NN9_mask_Large.tif','GTiff')
nnmsk_9 <- raster('KPARKK_NN9_mask_Large.tif')

par(mfrow=c(1,2))
plot(kpar_KK_ras,col=viridis(100),zlim=c(0.08,.35))
plot(coastline_KK,col='black',add=T,legend = F)

plot(nnmsk_9,col=viridis(100),zlim=c(0.08,.35))
plot(coastline_KK,col='black',add=T,legend = F)
Figure 3 - KdPAR after NN interpolation (ngb 9).

Figure 3 - KdPAR after NN interpolation (ngb 9).

NN 25 

gs <- gstat(formula=DownwellingPARattenuation~1, locations=kpar_KK_sp, nmax=25, set=list(idp = 0)) 
#using 9 neighboors/points + equal weight idp=0 
nn <- interpolate(grid_r, gs)

nnmsk <- mask(nn, dShore_bathy_KK)# Use Dshore to mask nn

writeRaster(nnmsk,'KPARKK_NN25_mask_Large.tif','GTiff')
nnmsk_25 <- raster('KPARKK_NN25_mask_Large.tif')

par(mfrow=c(1,2))
plot(kpar_KK_ras,col=viridis(100),zlim=c(0.08,.35))
plot(coastline_KK,col='black',add=T,legend = F)

plot(nnmsk_25,col=viridis(100),zlim=c(0.08,.35))
plot(coastline_KK,col='black',add=T,legend = F)
Figure 4 - KdPAR after NN interpolation (ngb 25).

Figure 4 - KdPAR after NN interpolation (ngb 25).

Inverse-Distance Weighting (IDW) 

gs <- gstat(formula=DownwellingPARattenuation~1, locations=kpar_KK_sp) #idp=2 by default
idw <- interpolate(grid_r, gs)

idwmsk <- mask(idw, dShore_bathy_KK)# Use Dshore to mask nn

writeRaster(idwmsk,'KPARKK_IDW_mask_Large.tif','GTiff')
idwmsk <- raster('KPARKK_IDW_mask_Large.tif')

par(mfrow=c(1,2))
plot(kpar_KK_ras,col=viridis(100),zlim=c(0.08,.35))
plot(coastline_KK,col='black',add=T,legend = F)

plot(idwmsk,col=viridis(100),zlim=c(0.08,0.35))
plot(coastline_KK,col='black',add=T,legend=F)
Figure 3 - KdPAR after IDW interpolation.

Figure 3 - KdPAR after IDW interpolation.

Comparison

Kd vs DShore 

kpar_KK <- kpar_KK_sf %>% as_tibble() %>% rename(KPAR_before = KPAR)
kpar_KK_merge <- kpar_KK_merge_sf %>% as_tibble() %>% rename(KPAR_after = KPAR)

kpar_kk_median_before <- kpar_KK %>% 
  dplyr::filter(DShore<15000) %>% 
  group_by(DShore) %>% 
  summarise(Kd_median_Before = median(KPAR_before,na.rm=T))

kpar_kk_median_mw <- kpar_KK_merge %>% 
  dplyr::filter(DShore<15000) %>% 
  group_by(DShore) %>% 
  summarise(Kd_median_MW = median(KPAR_after,na.rm=T))

## Interpolated Kd as tibble - NN5
nnmsk5_tb <- nnmsk_5 %>% 
  as_tibble(.,xy=T) %>%
  drop_na() %>% 
  rename(Kd_NN = cellvalue)

# Join Interpolated Kd and Dshore
kd_Dshore_nn5 <- inner_join(nnmsk5_tb,dShore_KK_tb)
## Joining, by = c("cellindex", "x", "y")
# Plot median Kd vs DShore
kd_Dshore_median_15km_nn5 <- kd_Dshore_nn5 %>% 
  dplyr::filter(DShore<15000) %>% 
  group_by(DShore) %>% 
  summarise(Kd_median_NN_5 = median(Kd_NN,na.rm=T))
##

## Interpolated Kd as tibble - NN9
nnmsk9_tb <- nnmsk_9 %>% 
  as_tibble(.,xy=T) %>%
  drop_na() %>% 
  rename(Kd_NN = cellvalue)

# Join Interpolated Kd and Dshore
kd_Dshore_nn9 <- inner_join(nnmsk9_tb,dShore_KK_tb)
## Joining, by = c("cellindex", "x", "y")
# Plot median Kd vs DShore
kd_Dshore_median_15km_nn9 <- kd_Dshore_nn9 %>% 
  dplyr::filter(DShore<15000) %>% 
  group_by(DShore) %>% 
  summarise(Kd_median_NN_9 = median(Kd_NN,na.rm=T))
##

## Interpolated Kd as tibble - NN25
nnmsk25_tb <- nnmsk_25 %>% 
  as_tibble(.,xy=T) %>%
  drop_na() %>% 
  rename(Kd_NN = cellvalue)

# Join Interpolated Kd and Dshore
kd_Dshore_nn25 <- inner_join(nnmsk25_tb,dShore_KK_tb)
## Joining, by = c("cellindex", "x", "y")
# Plot median Kd vs DShore
kd_Dshore_median_15km_nn25 <- kd_Dshore_nn25 %>% 
  dplyr::filter(DShore<15000) %>% 
  group_by(DShore) %>% 
  summarise(Kd_median_NN_25 = median(Kd_NN,na.rm=T))
##


## Interpolated Kd as tibble - IDW
idwmsk_tb <- idwmsk %>% 
  as_tibble(.,xy=T) %>%
  drop_na() %>% 
  rename(Kd_IDW = cellvalue)

# Join Interpolated Kd and Dshore
kd_Dshore_idw <- inner_join(idwmsk_tb,dShore_KK_tb)
## Joining, by = c("cellindex", "x", "y")
# Plot median Kd vs DShore
kd_Dshore_median_15km_idw <- kd_Dshore_idw %>% 
  dplyr::filter(DShore<15000) %>% 
  group_by(DShore) %>% 
  summarise(Kd_median_IDW = median(Kd_IDW,na.rm=T))
##

kd_Dshore_median_15km_idw_idps <- read_csv('IDW_DShore_Median_FewIDP_KPAR_KK.csv')
## 
## -- Column specification --------------------------------------------------------
## cols(
##   DShore = col_double(),
##   Kd_median_IDW_01 = col_double(),
##   Kd_median_IDW_05 = col_double(),
##   Kd_median_IDW_1 = col_double(),
##   Kd_median_IDW_2 = col_double(),
##   Kd_median_IDW_5 = col_double()
## )
kd_Dshore_median_mw_nn_idw <- left_join(kpar_kk_median_before,kpar_kk_median_mw,by='DShore') %>% 
  left_join(.,kd_Dshore_median_15km_nn5,by='DShore') %>% 
  left_join(.,kd_Dshore_median_15km_nn9,by='DShore') %>% 
  left_join(.,kd_Dshore_median_15km_nn25,by='DShore') %>% 
  left_join(.,kd_Dshore_median_15km_idw,by='DShore') %>% 
  left_join(.,kd_Dshore_median_15km_idw_idps,by='DShore') %>% 
  pivot_longer(-DShore,names_to='Method',values_to='KPAR')

p <- ggplot(kd_Dshore_median_mw_nn_idw,aes(DShore,KPAR,col=Method)) +
  geom_line() +
  theme_bw()
ggplotly(p)

Figure 4 - KdPAR vs DShore. 

#nn_mw_idw_tb_median <- nn_mw_idw_tb_long %>% 
#  group_by(DShore,Method) %>% 
#  summarise(Kd_median = median(Kd,na.rm=T), Kd_90 = quantile(Kd,.9,na.rm=T), Kd_10 = quantile(Kd,.1,na.rm=T))

#p <- ggplot(nn_mw_idw_tb_median,aes(x=DShore,y=Kd_median,col=Method)) +
#    geom_ribbon(aes(ymin=Kd_10, ymax=Kd_90), fill = "grey70", alpha=.1) +
#    geom_line()
#ggplotly(p)
  • Kd values really differ in the first ~2km.
  • IDW lower Kd values than NN and IDW + water clearer close to shore, wtf? Due to weight given to offshore and far pixels
  • NN seem more similar, makes sense
  • MW gives lower Kd values close to shore, not representative! Work to do to constrain loess fit when inshore estimate < offshore estimate?

MW vs NN 9 

#mw_sp <- st_read('KPARKK_After_Large_Win4.shp') %>% as_Spatial()
#mw_tb <- mw_sp %>% as_tibble() %>% 
#  rename(x = coords.x1, y = coords.x2, Kd_MW = KPAR)
mw_tb <- kpar_KK_estimates_sf %>%  #Only keep estimates
  dplyr::mutate(x = sf::st_coordinates(.)[,1],
                y = sf::st_coordinates(.)[,2]) %>% 
  rename(Kd_MW = KPAR) %>% 
  as_tibble() %>% dplyr::select(-geometry)

nn_mw_tb <- inner_join(kd_Dshore_nn9,mw_tb) #%>% 
## Joining, by = c("x", "y", "DShore")
  #mutate(diff = Kd_NN-Kd_MW) #%>% 
  #dplyr::filter(!between(diff,-1e-10,1e-10)) #ONly keep rows of differing estimates?
  #dplyr::distinct(Kd_NN,Kd_MW,.keep_all=T) #Delete rows with same Kd, mostly from true Kd, not estimates

nn_mw_idw_tb <- inner_join(nn_mw_tb,kd_Dshore_idw) %>% 
  dplyr::select(-cellindex) 
## Joining, by = c("cellindex", "x", "y", "DShore")
nn_mw_idw_tb_long <- nn_mw_idw_tb %>% 
  pivot_longer(-c(DShore,x,y),names_to = 'Method',values_to = 'Kd')

p <- ggplot(nn_mw_idw_tb,aes(x=Kd_NN,y=Kd_MW,col=DShore)) +
  geom_point() +
  geom_smooth(method = 'lm') +
  geom_abline(slope=1,intercept = 0,col='darkgrey') +
  scale_color_viridis_c(direction=-1) +
  theme_bw()
ggplotly(p)
## `geom_smooth()` using formula 'y ~ x'

Figure 6 - MW vs NN_9 interpolated values.

R2 of linear fit:

lm_fit <- lm(Kd_MW ~ Kd_NN, data = nn_mw_idw_tb) #%>% tidy() %>% kable()
#summary(lm_fit)#$r.squared

#lm(Kd_MW ~ Kd_NN, data = nn_mw_idw_tb) %>% kable()

stargazer(lm_fit,type = "text", style = 'default')#"html"
## 
## ===============================================
##                         Dependent variable:    
##                     ---------------------------
##                                Kd_MW           
## -----------------------------------------------
## Kd_NN                        0.706***          
##                               (0.038)          
##                                                
## Constant                     0.049***          
##                               (0.008)          
##                                                
## -----------------------------------------------
## Observations                   1,025           
## R2                             0.250           
## Adjusted R2                    0.249           
## Residual Std. Error      0.039 (df = 1023)     
## F Statistic          340.888*** (df = 1; 1023) 
## ===============================================
## Note:               *p<0.1; **p<0.05; ***p<0.01

MW vs IDW 

p <- ggplot(nn_mw_idw_tb,aes(x=Kd_IDW,y=Kd_MW,col=DShore)) +
  geom_point() +
  geom_smooth(method = 'lm') +
  geom_abline(slope=1,intercept = 0,col='darkgrey') +
  scale_color_viridis_c(direction=-1) +
  theme_bw()
ggplotly(p)
## `geom_smooth()` using formula 'y ~ x'

Figure 7 - MW vs IDW interpolated values.

R2 of linear fit:

lm_fit <- lm(Kd_MW ~ Kd_IDW, data = nn_mw_idw_tb) #%>% tidy() %>% kable()
stargazer(lm_fit,type = "text", style = 'default')#"html"
## 
## ===============================================
##                         Dependent variable:    
##                     ---------------------------
##                                Kd_MW           
## -----------------------------------------------
## Kd_IDW                       1.700***          
##                               (0.113)          
##                                                
## Constant                     -0.105***         
##                               (0.020)          
##                                                
## -----------------------------------------------
## Observations                   1,025           
## R2                             0.182           
## Adjusted R2                    0.181           
## Residual Std. Error      0.040 (df = 1023)     
## F Statistic          228.051*** (df = 1; 1023) 
## ===============================================
## Note:               *p<0.1; **p<0.05; ***p<0.01

NN 9 vs IDW 

p <- ggplot(nn_mw_idw_tb,aes(x=Kd_IDW,y=Kd_NN,col=DShore)) +
  geom_point() +
  geom_smooth(method = 'lm') +
  geom_abline(slope=1,intercept = 0,col='darkgrey') +
  scale_color_viridis_c(direction=-1) +
  theme_bw()
ggplotly(p)
## `geom_smooth()` using formula 'y ~ x'

Figure 8 - NN_9 vs IDW interpolated values.

R2 of linear fit:

lm_fit <- lm(Kd_NN ~ Kd_IDW, data = nn_mw_idw_tb) #%>% tidy() %>% kable()
stargazer(lm_fit,type = "text", style = 'default')#"html"
## 
## ===============================================
##                         Dependent variable:    
##                     ---------------------------
##                                Kd_NN           
## -----------------------------------------------
## Kd_IDW                       1.926***          
##                               (0.064)          
##                                                
## Constant                     -0.134***         
##                               (0.011)          
##                                                
## -----------------------------------------------
## Observations                   1,025           
## R2                             0.467           
## Adjusted R2                    0.466           
## Residual Std. Error      0.023 (df = 1023)     
## F Statistic          896.383*** (df = 1; 1023) 
## ===============================================
## Note:               *p<0.1; **p<0.05; ***p<0.01

Map of Relative Difference

IDW - MW/MW and NN - MW/MW

ndiff_sf <- nn_mw_idw_tb %>% 
  mutate(diff=(Kd_IDW-Kd_MW)/Kd_MW, diff2=(Kd_NN-Kd_MW)/Kd_MW) %>% 
  st_as_sf(.,coords=c('x','y'),crs=crs(kpar_KK_ras))

p <- ggplot() +
  geom_sf(data=coastline_KK) +
  geom_sf(data=ndiff_sf,aes(col=diff))+
  scale_colour_gradient2(low = "blue",mid = "white",high = "red",midpoint = 0,na.value = 'transparent') 
#ggplotly(p)

p2 <- ggplot() +
  geom_sf(data=coastline_KK) +
  geom_sf(data=ndiff_sf,aes(col=diff2))+
  scale_colour_gradient2(low = "blue",mid = "white",high = "red",midpoint = 0,na.value = 'transparent') 

subplot(p, p2, nrows = 1, margin = 0.04)

Figure 9 - Map of relative difference between IDW and MW (left), NN and MW (right).

Semi-Variogram exploration 

gs <- gstat(formula=DownwellingPARattenuation~1, locations=kpar_KK_sp) 
#idw <- interpolate(grid_r, gs)
#idwmsk <- mask(idw, dShore_bathy_KK)# Use Dshore to mask nn

system.time(
  gs_var <- variogram(gs, width=50)
)#80s for width=20 
write_csv(gs_var,'Variogram_KparKK_Width50.csv')
gs_var <- read.csv('Variogram_KparKK_Width50.csv')
class(gs_var) <- c("gstatVariogram", "data.frame")
plot(gs_var)
Figure 10 - Semi-Variogram

Figure 10 - Semi-Variogram

The Semi-Variogram does not level off, sign of anisotropy?

Let’s have a look at directional Semi-Variograms 45 degrees would be along-shore, 135 cross-shore direction.

var_dir <- variogram(DownwellingPARattenuation~1, kpar_KK_sp, alpha = c(45, 90, 135, 180), width=50)

write_csv(var_dir,'Directional_Variogram_KparKK_Width50.csv')
var_dir <- read.csv('Directional_Variogram_KparKK_Width50.csv')
class(var_dir) <- c("gstatVariogram", "data.frame")
plot(var_dir)
Figure 11 - Directional Semi-Variogram.

Figure 11 - Directional Semi-Variogram.

Different ranges across different directions = Anisotropy!

Bibliography