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)
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)
Last updated
