# 5. Extract ocean color data in optically-shallow waters 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:\ ## 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)` ![](https://2946503769-files.gitbook.io/~/files/v0/b/gitbook-legacy-files/o/assets%2F-LylLNCSXaUER_FiqDSx%2F-MFpnA4ciLLD9ty71bOg%2F-MFpnhWAI6PV6ZAf79RN%2Fimage.png?alt=media\&token=9f1a0bc8-f26a-429b-85de-a544037d50f5) ![](https://2946503769-files.gitbook.io/~/files/v0/b/gitbook-legacy-files/o/assets%2F-LylLNCSXaUER_FiqDSx%2F-MFpnA4ciLLD9ty71bOg%2F-MFpnt9ceeHs5CSDh58c%2Fimage.png?alt=media\&token=80a26e03-0068-4a20-888f-0a2421100da4) ## 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)` ![](https://2946503769-files.gitbook.io/~/files/v0/b/gitbook-legacy-files/o/assets%2F-LylLNCSXaUER_FiqDSx%2F-M-lUvWSXcdHu8tA1iTZ%2F-M-lVN3kzTfasKBhLW5N%2Fimage.png?alt=media\&token=bd6a2b07-d988-4f76-b311-307ed2867b24) `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](https://2946503769-files.gitbook.io/~/files/v0/b/gitbook-legacy-files/o/assets%2F-LylLNCSXaUER_FiqDSx%2F-M-lVPHtjq6PPWlcbR8E%2F-M-lVb-BsUygFy0pMqRI%2Fimage.png?alt=media\&token=42bce899-9278-4abe-a177-e63e5836a22f) 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)` ![](https://2946503769-files.gitbook.io/~/files/v0/b/gitbook-legacy-files/o/assets%2F-LylLNCSXaUER_FiqDSx%2F-M-lWh6D2XT17zHotGM3%2F-M-lWtrI4hwTd6ut-jFQ%2Fimage.png?alt=media\&token=92944eef-d20d-492b-a6f0-eb06477d307b) ## 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)` ![](https://2946503769-files.gitbook.io/~/files/v0/b/gitbook-legacy-files/o/assets%2F-LylLNCSXaUER_FiqDSx%2F-MFpo3839MTTSAWDxmQ0%2F-MFpoNFnr6GhozBlWlep%2Fimage.png?alt=media\&token=6c3050f3-907d-4124-a889-d4a77611b094) ## 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](https://oceanwatch.pifsc.noaa.gov/files/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()` ![](https://2946503769-files.gitbook.io/~/files/v0/b/gitbook-legacy-files/o/assets%2F-LylLNCSXaUER_FiqDSx%2F-MFpo3839MTTSAWDxmQ0%2F-MFposnTkEuaslg5Wh8V%2Fimage.png?alt=media\&token=70856871-77d9-49cd-a587-d84b48e5b896)