5. Extract ocean color data in optically-shallow waters
Code written by Thomas Oliver, NOAA PIFSC.
Last updated
Code written by Thomas Oliver, NOAA PIFSC.
Last updated
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
library(raster)
library(sp)
library(rerddap)
library(lubridate)
Wake=read.csv('wake.csv')
Wake
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
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" )
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)
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:
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)
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)
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
Note: the distances are in km.
Let's save those chl-a values in our Wake data frame.
Wake$VIIRS_CHLA=EEdf$values
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()
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
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
