5. Extract ocean color data in optically-shallow waters

Code written by Thomas Oliver, NOAA PIFSC.

Remotely sensed ocean color algorithms are calibrated for optically-deep waters, where the signal received by the satellite sensor originates from the water column without any bottom contribution.

Optically shallow waters are those in which light reflected off the seafloor contributes significantly to the water-leaving signal, such as coral reefs, atolls, lagoons. This is known to affect geophysical variables derived by ocean-color algorithms, often leading to biased values in chlorophyll-a concentration for example.

In the tropical Pacific, optically-deep waters are typically deeper than 15 – 30 m. It is recommended to remove shallow-pixels (<30m depth) from the study area before computing ocean color metrics.

In this tutorial, we will extract chl-a concentration data at survey locations around Wake island in the Pacific Ocean. For survey points located in waters shallower than 30m, we will find chl-a pixels in deeper water and extract those values instead.

The survey locations can be downloaded here: https://oceanwatch.pifsc.noaa.gov/files/wake.csv

Load packages

library(raster) library(sp) library(rerddap) library(lubridate)

Load survey data

Wake=read.csv('wake.csv')

Wake

	Deploy_Longitude	Deploy_Latitude	Year	Island	Site
1	166.6073	19.29178	2014	WAK	WAK06
2	166.5983	19.31627	2014	WAK	WAK08
3	166.6516	19.27068	2014	WAK	WAK09
4	166.5983	19.31621	2017	WAK	WAK08
5	166.6272	19.31605	2017	WAK	WAK23
6	166.6278	19.28066	2017	WAK	WAK01
7	166.6511	19.30614	2017	WAK	WAK24

Let's transform our survey data into a SpatialPointsDataFrame, which makes it easier to work with geospatial rasters:

coordinates(Wake)<- ~Deploy_Longitude+Deploy_Latitude

Wake class : SpatialPointsDataFrame features : 7 extent : 166.5983, 166.6516, 19.27068, 19.31627 (xmin, xmax, ymin, ymax) crs : NA variables : 3 names : Year, Island, Site min values : 2014, WAK, WAK01 max values : 2017, WAK, WAK24

Download bathymetry and chl-a concentration data for the survey area

scale=.05 CW_u='https://coastwatch.pfeg.noaa.gov/erddap/' ETOPO1_id='etopo180' ETOPO1_info=info(datasetid = ETOPO1_id,url = CW_u) WakeBathy=griddap(ETOPO1_info,url=CW_u, latitude=(range(Wake$Deploy_Latitude)+c(-1,1)*scale), longitude=(range(Wake$Deploy_Longitude)+c(-1,1)*scale), store=disk(),fmt = "nc")

OW_u='https://oceanwatch.pifsc.noaa.gov/erddap/' VIIRS_id='noaa_snpp_chla_monthly' VIIRS_info=info(datasetid = VIIRS_id,url = OW_u) var=VIIRS_info$variable$variable_name

WakeVIIRS=griddap(url=OW_u, VIIRS_id, time = c('2016-04-01', '2017-04-01'), latitude = range(Wake$Deploy_Latitude)+c(-1,1)*scale, longitude = range(Wake$Deploy_Longitude)+c(-1,1)*scale, fields = var[1], store=disk(),fmt = "nc" )

Convert bathymetry and chl-a data to rasters

rWB=raster(WakeBathy$summary$filename)

WakeVIIRS contains multiple timesteps, so we need to use the stack function:

rVI=stack(WakeVIIRS$summary$filename,varname="chlor_a") rVI.mean=mean(rVI,na.rm=TRUE)

Let's look at our rasters:

blue.col <- colorRampPalette(c("darkblue", "lightblue")) chl.col=colorRampPalette(c("#00ffff","#00e600")) plot(rWB,main="ETOPO1 Bathymetry Raster",col=blue.col(255)) contour(rWB,levels=c(-30,-1000,-2000),add=TRUE) plot(Wake,add=TRUE,pch=16) plot(log(rVI.mean),main="VIIRS CHLA Raster (log scale)",col=chl.col(255)) plot(Wake,add=TRUE,pch=16)

Depth at survey sites

Wake$Depth=raster::extract(x=rWB,y=Wake) Wake$Depth [1] 3 53 2 53 2 2 5

The bathymetry data has ~2-km resolution, whereas the chl-a data has a ~4-km resolution. We are going to build a function to look at the chl-a data in our survey area and determine whether to call each chl-a pixel "shallow" based on how many shallow (<30m) bathymetry pixels it overlaps. 1. Convert the depth raster to a SpatialPointsDataFrame : extent(rWB)=extent(rVI.mean) spWB=data.frame(rasterToPoints(rWB)) coordinates(spWB)=~x+y plot(rWB,main="ETOPO1 Bathymetry Raster",col=blue.col(255)) contour(rWB,levels=c(-30),add=TRUE) plot(spWB,add=TRUE,pch=20)

plot(log(rVI.mean),main="VIIRS Chl-a Raster",col=chl.col(255)) contour(rWB,levels=c(-30),add=TRUE) plot(spWB,add=TRUE,pch=20)

Then we plot those SpatialPoints on top of the chl-a raster:

Each chl-a pixel overlaps 4 bathymetry pixels

2. Define a function to consider a pixel necessary to mask:

count_shallow_pixels=function(depths,threshold=-30,na.rm=T){ return(length(which(depths>threshold))) } count_all=function(x,na.rm=T){ return(length(x)) }

3. Build a raster of the chl-a grid, using the function to count how many (smaller) depth pixels in each (larger) Chla pixel are "too shallow" (out of 4):

rVI.N_SHALLOW=rasterize(x = spWB,y=rVI.mean,field="Altitude",fun=count_shallow_pixels) plot(rVI.N_SHALLOW,main="N Shallow Pixels",asp=1) contour(rWB,levels=c(-30),add=TRUE) plot(spWB,add=T,pch=20)

Masking

Decide what threshold means a pixel is 'bad', then generate the masked chl-a layer. For example, let's mask all chl-a pixels that overlap 2 or more shallow bathymetry pixels:

rVI.mean.masked=rVI.mean rVI.mean.masked[rVI.N_SHALLOW>=2]=NA plot(log(rVI.mean.masked),main="Masked VIIRS Chl-A",col=chl.col(255)) plot(Wake,add=T,pch=20)

Values for survey locations in masked pixels

For the points located in shallow pixels, we now need to decide which value of chl-a to extract using chl-a pixels that are further off-shore.

The following function takes a spatial raster, and a spatial data frame (SpDF) of in situ points. Then it will fill any NA value in the SpDF with the first-discovered non-NA values from the raster:

ExpandingExtract=function(r,SpDF,Dists=c(500,1000,2000,4000,8000)){ OutDF=data.frame(values=rep(NA,nrow(SpDF)), Dist=rep(NA,nrow(SpDF)), N=rep(NA,nrow(SpDF))) nDists=length(Dists) cnt=1 NAi=which(is.na(OutDF$values)) NAsLeft=length(NAi)>0 while(cnt<=nDists&NAsLeft){ NAi=which(is.na(OutDF$values)) pull=raster::extract(x=r,y=SpDF[NAi,], buffer=Dists[cnt], small=TRUE, na.rm=TRUE) OutDF$values[NAi]=unlist(lapply(pull,mean,na.rm=TRUE)) OutDF$Dist[NAi]=Dists[cnt] NAi=which(is.na(OutDF$values)) NAsLeft=length(NAi)>0 cnt=cnt+1 } return(OutDF) }

Let's call our ExpandingExtract function: EEdf=ExpandingExtract(rVI.mean.masked,Wake,Dists=c(500,1000,2000,4000,6000)) EEdf

	values	Dist	N
1	0.107738	500	NA
2	0.107738	4000	NA
3	0.049559	500	NA
4	0.107738	4000	NA
5	0.064098	6000	NA
6	0.077915	4000	NA
7	0.050901	500	NA

Note: the distances are in km.

Let's save those chl-a values in our Wake data frame.

Wake$VIIRS_CHLA=EEdf$values

Plot

Plot color scale using the scale.R file:

library(maps) library(mapdata) library(maptools) source('scale.R')

x11(width=4.8,height=3.55)

xlim=c(166.58,166.7) ylim=c(19.25,19.35)

layout(matrix(c(1,2,3,0,4,0), nrow=1, ncol=2), widths=c(4,1), heights=4) layout.show(2)

par(mar=c(3,3,3,1))

land <- maps::map('worldHires', fill=TRUE, xlim=xlim, ylim=ylim, plot=FALSE)

ids <- sapply(strsplit(land$names, ":"), function(x) x[1]) bPols <- map2SpatialPolygons(land, IDs=ids, proj4string=CRS("+proj=longlat +datum=WGS84"))

plot(bPols, col="grey", axes=FALSE,xlim=xlim,ylim=ylim,cex.axis=3,xaxs='i',yaxs='i',asp=1)

x=seq(166.5,166.7,0.05) axis(1,x,x) y=seq(19.2,19.4,0.05) axis(2,y,y) box()

breaks=seq(-3.01,-2.22,0.01) n=length(breaks)-1 c=chl.col(n) for (i in 1:length(Wake$VIIRS_CHLA)) { I=which(breaks>log(Wake$VIIRS_CHLA[i])) points(coordinates(Wake)[i,1],coordinates(Wake)[i,2], pch=20,col=c[I[1]-1],cex=3) }

par(mar=c(3,1,3,3)) image.scale(Wake$VIIRS_CHLA, col=c, breaks=breaks, horiz=FALSE, yaxt="n",xlab='',ylab='',main='log(chl)') axis(4,las=1) box()