AIM : Compare global estimates of benthic light (bottompar, Gattuso et al. (2020)) vs NIWA Ebed. Climatology 1998-2018 for benthic PAR and 07_2002-06_2020 for Ebed. Look at points of difference.

Table of content

  1. Global BottomPar vs Ebed

  2. KdPAR Global vs KdPAR NIWA

  1. Depth Global vs NIWA Depth

  2. PAR Global vs NIWA PAR

  3. Bibliography

R Packages

#knitr::opts_chunk$set(echo = TRUE)
library(raster)
## Loading required package: sp
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

Global BottomPar vs Ebed

The Global BottomPAR climatology (1998-2018) has been downloaded using the R package CoastalLight:: (Gentili and Gattuso 2021). It is projecting on the same grid and projection reference system than Ebed (NZ UTM) using bilinear interpolation (raster::projectRaster, (Hijmans 2020)). Values of BottomPAR and Ebed are joined and compared for the same pixels (using the tidyverse ecosystem of packages, (Wickham et al. 2019)).

## NIWA Ebed
ebed <- raster('./CoastalLight.data/A20021822020182_YC_EBED_exp_coastal_v01.nc')
## Loading required namespace: ncdf4
## PAR BOTTOM - Gattuso2020 - 1998 to 2018
parbottom <- raster('./CoastalLight.data/1998_2018_parbottom_YC_NZ_Gattuso2020.nc')

system.time(
  parbottom_nztm <- projectRaster(parbottom,ebed,crs(ebed), method="bilinear")
)#13s
##    user  system elapsed 
##   13.66    1.92   16.05
max_pb <- max(values(parbottom_nztm),na.rm=T)
max_ebed <- max(values(ebed),na.rm=T)

einstein_values <- c(0.0001,0.01,0.1,1,10,max_pb)
einstein_values_log <- log10(einstein_values)

einstein_values_ebed <- c(0.0001,0.01,0.1,1,10,max_ebed)
einstein_values_log_ebed <- log10(einstein_values)

par(mfrow=c(1,2))
plot(log10(parbottom_nztm),col=viridis(100),legend=FALSE,zlim=c(einstein_values_log[1],einstein_values_log[length(einstein_values_log)]))
plot(log10(parbottom_nztm),zlim=c(einstein_values_log[1],einstein_values_log[length(einstein_values_log)]),col=viridis(100), legend.only=TRUE,legend.width=1, legend.shrink=0.75,
     axis.args=list(at=einstein_values_log,
                    labels=einstein_values, 
                    cex.axis=0.6),
     legend.args=list(text='Bottom PAR E/m2/d', side=4, font=2, line=2.5, cex=0.8))

plot(log10(ebed),col=viridis(100),legend=FALSE,zlim=c(einstein_values_log_ebed[1],einstein_values_log_ebed[length(einstein_values_log_ebed)]))
plot(log10(ebed),zlim=c(einstein_values_log_ebed[1],einstein_values_log_ebed[length(einstein_values_log_ebed)]),col=viridis(100), legend.only=TRUE,legend.width=1, legend.shrink=0.75,
     axis.args=list(at=einstein_values_log_ebed,
                    labels=einstein_values_ebed, 
                    cex.axis=0.6),
     legend.args=list(text='Ebed E/m2/d', side=4, font=2, line=2.5, cex=0.8))
Figure 1 - Map of Global Benthic PAR estimates (PAR Bottom) vs NIWA Ebed.

Figure 1 - Map of Global Benthic PAR estimates (PAR Bottom) vs NIWA Ebed. 

parbottom_tb <- tabularaster::as_tibble(parbottom_nztm,xy=T,cell=F) %>% 
  drop_na() %>% 
  rename(ParBottom=cellvalue)

ebed_tb <- tabularaster::as_tibble(ebed,xy=T,cell=F) %>% 
  drop_na() %>% 
  filter(!cellvalue==0) %>% 
  rename(Ebed=cellvalue)

ebed_parbottom_join <- inner_join(ebed_tb,parbottom_tb,by=c('x','y'))

## Histograms of PARb vs EBED
ebed_parbottom_join_long <- ebed_parbottom_join %>% 
  pivot_longer(-c(x,y),values_to ='SeabedLight',names_to='Source') %>% 
  group_by(Source)

p <- ggplot(ebed_parbottom_join_long,aes(x=SeabedLight,fill=Source)) +
  geom_histogram(position = 'identity', bins=100, alpha=0.6) +
  scale_x_continuous(trans='log10',limits = c(1e-8,100)) +
  theme_bw()
ggplotly(p)

Figure 2 - Distribution of Global Benthic PAR estimates (PAR Bottom) vs NIWA Ebed. 

## ScatterPlots
p <- ggplot(ebed_parbottom_join,aes(x=Ebed,y=ParBottom)) +
  geom_hex(bins = 70) +
  scale_fill_continuous(type = "viridis") +
  scale_x_continuous(trans='log10',limits = c(1e-8,100)) +
  scale_y_continuous(trans='log10',limits = c(1e-8,100)) +
  geom_abline(slope=1,intercept = 0,linetype = "dashed",col='darkred',size=2) +
  theme_bw()
ggplotly(p)

Figure 3 - Scatterplot of Global Benthic PAR estimates (PAR Bottom) vs NIWA Ebed.

Global estimates of benthic PAR seem to “overestimate” NIWA Ebed, and in greater proportion at low values (<.01). Why is that?

Let’s investigate the first 3 points in the following.

KdPAR Global vs KdPAR NIWA

Histogram and scatterplot

The Global KdPAR climatology (1998-2018) has been downloaded using the R package CoastalLight:: (Gentili and Gattuso 2021). It is projecting on the same grid and projection reference system than NIWA KdPAR (NZ UTM, 07-2002 to 06-2020 climatology) using bilinear interpolation (raster::projectRaster, (Hijmans 2020)). Values of Global KdPAR and NIWA KdPAR are joined and compared for the same pixels (using the tidyverse ecosystem of packages, (Wickham et al. 2019)).

kpar_niwa <- raster('./CoastalLight.data/A20021822020182_YC_KPAR_exp_coastal_v01.nc')

kpar_global <- raster('./CoastalLight.data/1998_2018_kpar_YC_NZ_Gattuso2020.nc')

system.time(
  kpar_global_nztm <- projectRaster(kpar_global,kpar_niwa,crs(kpar_niwa), method="bilinear")
)#15s
##    user  system elapsed 
##   11.89    2.19   14.75
kpar_global_tb <- tabularaster::as_tibble(kpar_global_nztm,xy=T,cell=F) %>% 
  drop_na() %>% 
  rename(kpar_global=cellvalue)

kpar_niwa_tb <- tabularaster::as_tibble(kpar_niwa,xy=T,cell=F) %>% 
  drop_na() %>% 
  rename(kpar_niwa=cellvalue)

kpar_join <- inner_join(kpar_global_tb,kpar_niwa_tb,by=c('x','y'))

## Histograms of Global vs NIWA
kpar_join_long <- kpar_join %>% 
  pivot_longer(-c(x,y),values_to ='KdPAR',names_to='Source') %>% 
  group_by(Source)

p <- ggplot(kpar_join_long,aes(x=KdPAR,fill=Source)) +
  geom_histogram(position = 'identity', bins=100, alpha=0.6) +
  #scale_x_continuous(trans='log10',limits = c(1e-8,100)) +
  xlim(c(0,.5)) +
  theme_bw()
ggplotly(p)

Figure 4 - Distribution of Global KdPAR estimates vs NIWA KdPAR 

## ScatterPlots
p <- ggplot(kpar_join,aes(x=kpar_niwa,y=kpar_global)) +
  geom_hex(bins = 70) +
  scale_fill_continuous(type = "viridis") +
  xlim(c(0.05,0.5)) + ylim(c(0.05,0.5)) +
  #scale_x_continuous(trans='log10',limits = c(1e-8,100)) +
  #scale_y_continuous(trans='log10',limits = c(1e-8,100)) +
  geom_abline(slope=1,intercept = 0,linetype = "dashed",col='grey',size=2) +
  theme_bw()
ggplotly(p)

Figure 5 - Scatterplot of Global KdPAR estimates vs NIWA Ebed.

NIWA KdPAR >> Global KdPAR. So global spatial resolution which “misses” very coastal and higher KdPAR values could potentially explain the underestimation of KdPAR leading to the overestimation of BottomPAR.

KdPAR vs Distance to Shore

Let’s have a look at the relationship between KdPAR and Distance to Shore

dshore <- raster('NZ_Dist2Shore.tif')

dshore_tb <- tabularaster::as_tibble(dshore,xy=T,cell=F) %>% 
  drop_na() %>% 
  rename(DShore=cellvalue)

kpar_join_long_dshore <- inner_join(kpar_join_long,dshore_tb,by=c('x','y'))

kpar_join_long_dshore_median <- kpar_join_long_dshore %>% 
  group_by(DShore,Source) %>% 
  summarise(KdPAR_median = median(KdPAR,na.rm=T))
## `summarise()` has grouped output by 'DShore'. You can override using the `.groups` argument.
p <- ggplot(kpar_join_long_dshore_median,aes(x=DShore,y=KdPAR_median,col=Source)) +
  geom_point() +
  theme_bw()
ggplotly(p)

Figure 6 - Global KdPAR estimates and NIWA KdPAR vs Distance to Shore (m)

Similar relationship of KdPAR vs Dshore for both global and NIWA products, i.e. gradient with Dshore. BUT consistent difference, why? Due to algorithms?

–> After meeting with Matt and Mark: there could be issues in the log-scaling of case 1- case 2 waters blending algorithm.

kpar_join_dshore <- inner_join(kpar_join,dshore_tb,by=c('x','y'))

kpar_join_dshore_median <- kpar_join_dshore %>% 
  group_by(DShore) %>% 
  summarise(KdPAR_global_median = median(kpar_global,na.rm=T),
            KdPAR_niwa_median = median(kpar_niwa,na.rm=T))

p <- ggplot(kpar_join_dshore_median,aes(x=KdPAR_niwa_median,y=KdPAR_global_median,col=DShore)) +
  geom_point() +
  scale_colour_continuous(type = "viridis",direction =-1) +
  geom_smooth(method = lm) +
  geom_abline(slope=1,intercept = 0,linetype = "dashed",col='grey',size=2) +
  xlim(0.1,0.35)  + ylim(0.1,0.35) +
  theme_bw()
ggplotly(p)
## `geom_smooth()` using formula 'y ~ x'

Figure 6b - Scatterplot of the median Kd_niwa vs Kd_global. 

--> linear regression is parallel to diagonal
    ○ Correction possible?
    ○ Offset due to log-scaling issue (case 1 not working?)
    ○ Good to see points very close to shore stepping away from linear relationship ==> Underestimate of coastal Kd in Gattuso (so over-estimate Ebed)

Depth Global vs NIWA Depth

The Global Depth has been downloaded using the R package CoastalLight:: (Gentili and Gattuso 2021). It is projecting on the same grid and projection reference system than NIWA Depth (NZ UTM) using bilinear interpolation (raster::projectRaster, (Hijmans 2020)). Values of Global Depth and NIWA Depth are joined and compared for the same pixels (using the tidyverse ecosystem of packages, (Wickham et al. 2019)).

depth_niwa <- raster('./CoastalLight.data/nzbathymetry_2016_transform.tif')
depth_niwa[values(depth_niwa)>0] <- NA

depth_global <- raster('./CoastalLight.data/Depth_NZ_Gattuso2020.nc')

system.time(
  depth_global_nztm <- projectRaster(depth_global,depth_niwa,crs(depth_niwa), method="bilinear")
)#15s
##    user  system elapsed 
##   12.83    1.65   14.64
depth_global_tb <- tabularaster::as_tibble(depth_global_nztm,xy=T,cell=F) %>% 
  drop_na() %>% 
  rename(depth_global=cellvalue)

depth_niwa_tb <- tabularaster::as_tibble(depth_niwa,xy=T,cell=F) %>% 
  drop_na() %>% 
  rename(depth_niwa=cellvalue)

depth_join <- inner_join(depth_global_tb,depth_niwa_tb,by=c('x','y'))

## Histograms of Global vs NIWA
depth_join_long <- depth_join %>% 
  pivot_longer(-c(x,y),values_to ='Depth',names_to='Source') %>% 
  group_by(Source)

p <- ggplot(depth_join_long,aes(x=Depth,fill=Source)) +
  geom_histogram(position = 'identity', bins=100, alpha=0.6) +
  #scale_x_continuous(trans='log10',limits = c(1e-8,100)) +
  theme_bw()
ggplotly(p)

Figure 6 - Distribution of Global Bathymetry vs NIWA Bathymetry. 

## ScatterPlots
p <- ggplot(depth_join,aes(x=depth_niwa,y=depth_global)) +
  geom_hex(bins = 70) +
  scale_fill_continuous(type = "viridis") +
  xlim(c(-200,0)) + ylim(c(-200,0)) +
  #scale_x_continuous(trans='log10',limits = c(1e-8,100)) +
  #scale_y_continuous(trans='log10',limits = c(1e-8,100)) +
  geom_abline(slope=1,intercept = 0,linetype = "dashed",col='grey',size=2) +
  theme_bw()
ggplotly(p)

Figure 7 - Scatterplot of Global Bathymetry vs NIWA Bathymetry.

The bathymetry used for the global estimates of Benthic PAR is very similar than the bathymetry used for NIWA Ebed. So difference in BottomPAR and Ebed not due to difference in Bathymetry.

PAR Global vs NIWA PAR

The Global PAR climatology (1998-2018) has been downloaded using the R package CoastalLight:: (Gentili and Gattuso 2021). The “NIWA” PAR has been downloaded from NASA Ocean Colour, Monthly composite, level 3, 4km before being averaged over the 08-2002 to 07-2020 period (Frouin et al. (2012)). It is projected on the same grid and projection reference system than Global KdPAR (lon lat, 0.004166667 degree resolution) using bilinear interpolation (raster::projectRaster, (Hijmans 2020)). Values of Global PAR and NIWA PAR are joined and compared for the same pixels (using the tidyverse ecosystem of packages, (Wickham et al. 2019)).

par_niwa <- raster('./CoastalLight.data/082002_07_2020_PAR_4km_Mean.tif')

par_global <- raster('./CoastalLight.data/1998_2018_par_YC_NZ_Gattuso2020.nc')

system.time(
  par_niwa_int <- projectRaster(par_niwa,par_global,crs(par_global), method="bilinear")
)#15s
##    user  system elapsed 
##   23.66   12.82   40.45
par_global_tb <- tabularaster::as_tibble(par_global,xy=T,cell=F) %>% 
  drop_na() %>% 
  rename(par_global=cellvalue)

par_niwa_tb <- tabularaster::as_tibble(par_niwa_int,xy=T,cell=F) %>% 
  drop_na() %>% 
  rename(par_niwa=cellvalue)

par_join <- inner_join(par_global_tb,par_niwa_tb,by=c('x','y'))

## Histograms of Global vs NIWA
par_join_long <- par_join %>% 
  pivot_longer(-c(x,y),values_to ='PAR',names_to='Source') %>% 
  group_by(Source)

p <- ggplot(par_join_long,aes(x=PAR,fill=Source)) +
  geom_histogram(position = 'identity', bins=100, alpha=0.6) +
  theme_bw()
ggplotly(p)

Figure 8 - Distribution of Global PAR estimates vs NIWA PAR 

## ScatterPlots
p <- ggplot(par_join,aes(x=par_niwa,y=par_global)) +
  geom_hex(bins = 70) +
  scale_fill_continuous(type = "viridis") +
  geom_abline(slope=1,intercept = 0,linetype = "dashed",col='grey',size=2) +
  theme_bw()
ggplotly(p)

Figure 9 - Scatterplot of Global PAR estimates vs NIWA PAR

Global and NIWA PAR quite similar. Slightly higher estimate for NIWA PAR than Global PAR at higher values, and lower estimates for lower values. It does not seem like that PAR is affecting much the difference in benthic light estimates.

Bibliography

Frouin, Robert, John McPherson, Kyozo Ueyoshi, and Bryan A Franz. 2012. “A Time Series of Photosynthetically Available Radiation at the Ocean Surface from Seawifs and Modis Data.” In Remote Sensing of the Marine Environment Ii, 8525:852519. International Society for Optics; Photonics.

Gattuso, Jean-Pierre, Bernard Gentili, David Antoine, and David Doxaran. 2020. “Global Distribution of Photosynthetically Available Radiation on the Seafloor.” Earth System Science Data 12 (3): 1697–1709.

Gentili, Bernard, and Jean-Pierre Gattuso. 2021. CoastalLight: Available Light on the Coastal Ocean [0-200m] Floor.

Hijmans, Robert J. 2020. Raster: Geographic Data Analysis and Modeling. https://CRAN.R-project.org/package=raster.

Wickham, Hadley, Mara Averick, Jennifer Bryan, Winston Chang, Lucy D’Agostino McGowan, Romain François, Garrett Grolemund, et al. 2019. “Welcome to the tidyverse.” Journal of Open Source Software 4 (43): 1686. https://doi.org/10.21105/joss.01686.