5. Extract ocean color data in optically-shallow waters - new version

Tom Oliver (NOAA/PIFSC), Melanie Abecassis (NOAA CoastWatch)

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")

For Chl-a, let's first extract a random time step:

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', '2016-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)

rVI=raster(WakeVIIRS$summary$filename,varname="chlor_a")

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),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) 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),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))) }

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,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.masked=rVI rVI.masked[rVI.N_SHALLOW>=2]=NA plot(log(rVI.masked),main="Masked VIIRS Chl-A",col=chl.col(255)) plot(Wake,add=T,pch=20,cex=2,col=2) contour(rWB,levels=c(-30),add=TRUE) text(Wake@coords[,1],Wake@coords[,2],row.names(Wake@coords),pos=4)

Values for survey locations in masked pixels

For the points located in shallow pixels, we now need to pick a new chl-a pixel to substitute for each masked one. To that end, for each survey location, we will select the closest unmasked chl-a pixel. We will then stored those new locations to download Chl-a data.

This will allow you to download Chl-a data for different time ranges for each location if needed.

Let's look at the example of survey location #5:

r=rVI.masked xy=Wake@coords i=5 d=replace(distanceFromPoints(r, xy[i,]), is.na(r), NA) plot(d) points(Wake@coords[i,1],Wake@coords[i,2],pch=20,col=2,cex=2) text(Wake@coords[i,1],Wake@coords[i,2],row.names(Wake@coords)[i],pos=4)

d is the distance of each chl-a pixel to survey location #5.

We need to select the pixel closest (minimum distance) to our location:

new_coords=xyFromCell(r,which.min(d)) plot(d) points(Wake@coords[i,1],Wake@coords[i,2],pch=20,col=2,cex=2) text(Wake@coords[i,1],Wake@coords[i,2],row.names(Wake@coords)[i],pos=4) points(new_coords[1],new_coords[2],pch=20,col=4,cex=2)

The blue dot identifies the closest chl-a pixel.

Now we need to do this for all our locations, and store the new coordinates in the Wake data frame.

new_coords=array(NA,dim=c(dim(xy)[1],2)) for (i in 1:dim(xy)[1]) { d=replace(distanceFromPoints(r, xy[i,]), is.na(r), NA) new_coords[i,]=xyFromCell(r,which.min(d)) } colnames(new_coords)=c("new_lon","new_lat") Wake=cbind(data.frame(Wake),new_coords)

Wake

	Deploy_Longitude	Deploy_Latitude	Year	Island	Site	Depth	optional	new_lon	new_lat
1	166.6073	19.29178	2014	WAK	WAK06	3	TRUE	166.5938	19.29375
2	166.5983	19.31627	2014	WAK	WAK08	53	TRUE	166.5938	19.29375
3	166.6516	19.27068	2014	WAK	WAK09	2	TRUE	166.6688	19.25625
4	166.5983	19.31621	2017	WAK	WAK08	53	TRUE	166.5938	19.29375
5	166.6272	19.31605	2017	WAK	WAK23	2	TRUE	166.5938	19.29375
6	166.6278	19.28066	2017	WAK	WAK01	2	TRUE	166.6313	19.25625
7	166.6511	19.30614	2017	WAK	WAK24	5	TRUE	166.6688	19.29375

plot(log(rVI.masked),main="Masked VIIRS Chl-A",col=chl.col(255)) contour(rWB,levels=c(-30),add=TRUE) points(Wake$new_lon,Wake$new_lat,pch=20,col=2,cex=2) text(Wake$new_lon,Wake$new_lat,row.names(Wake),pos=4)