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.
library(raster)
library(sp)
library(rerddap)
library(lubridate)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_LatitudeWake
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, WAK24scale=.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_nameWakeVIIRS=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" )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)
Wake$Depth=raster::extract(x=rWB,y=Wake)
Wake$Depth
[1] 3 53 2 53 2 2 5The 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)
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([email protected][,1],[email protected][,2],row.names([email protected]),pos=4)
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
[email protected]
i=5
d=replace(distanceFromPoints(r, xy[i,]), is.na(r), NA)
plot(d)
points([email protected][i,1],[email protected][i,2],pch=20,col=2,cex=2)
text([email protected][i,1],[email protected][i,2],row.names([email protected])[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([email protected][i,1],[email protected][i,2],pch=20,col=2,cex=2)
text([email protected][i,1],[email protected][i,2],row.names([email protected])[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)
Last modified 2yr ago