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

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 <a href="#load-packages" id="load-packages"></a>

`library(raster)` \
`library(sp)` \
`library(rerddap)` \
`library(lubridate)`

## Load survey data <a href="#load-survey-data" id="load-survey-data"></a>

`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 <a href="#download-bathymetry-and-chl-a-concentration-data-for-the-survey-area" id="download-bathymetry-and-chl-a-concentration-data-for-the-survey-area"></a>

`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 <a href="#convert-bathymetry-and-chl-a-data-to-rasters" id="convert-bathymetry-and-chl-a-data-to-rasters"></a>

`rWB=raster(WakeBathy$summary$filename)`

`rVI=raster(WakeVIIRS$summary$filename,varname="chlor_a")`&#x20;

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

![](https://gblobscdn.gitbook.com/assets%2F-LylLNCSXaUER_FiqDSx%2F-MFpnA4ciLLD9ty71bOg%2F-MFpnhWAI6PV6ZAf79RN%2Fimage.png?alt=media\&token=9f1a0bc8-f26a-429b-85de-a544037d50f5)

![](https://gblobscdn.gitbook.com/assets%2F-LylLNCSXaUER_FiqDSx%2F-MFpnA4ciLLD9ty71bOg%2F-MFpnt9ceeHs5CSDh58c%2Fimage.png?alt=media\&token=80a26e03-0068-4a20-888f-0a2421100da4)

## Depth at survey sites <a href="#depth-at-survey-sites" id="depth-at-survey-sites"></a>

`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.&#x20;

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.&#x20;

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

![](https://gblobscdn.gitbook.com/assets%2F-LylLNCSXaUER_FiqDSx%2F-M-lUvWSXcdHu8tA1iTZ%2F-M-lVN3kzTfasKBhLW5N%2Fimage.png?alt=media\&token=bd6a2b07-d988-4f76-b311-307ed2867b24)

`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:

![](https://gblobscdn.gitbook.com/assets%2F-LylLNCSXaUER_FiqDSx%2F-M-lVPHtjq6PPWlcbR8E%2F-M-lVb-BsUygFy0pMqRI%2Fimage.png?alt=media\&token=42bce899-9278-4abe-a177-e63e5836a22f)

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){` \
&#x20;  `return(length(which(depths>threshold)))` \
`}`&#x20;

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

![](https://gblobscdn.gitbook.com/assets%2F-LylLNCSXaUER_FiqDSx%2F-M-lWh6D2XT17zHotGM3%2F-M-lWtrI4hwTd6ut-jFQ%2Fimage.png?alt=media\&token=92944eef-d20d-492b-a6f0-eb06477d307b)

## Masking <a href="#masking" id="masking"></a>

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

![](https://2946503769-files.gitbook.io/~/files/v0/b/gitbook-legacy-files/o/assets%2F-LylLNCSXaUER_FiqDSx%2F-MG9lrDf34nxSIYQ3FAp%2F-MG9oFIeft7nmDC-E66s%2Fimage.png?alt=media\&token=1eaeca3d-bf0f-4cad-bd91-05056fb880d3)

## Values for survey locations in masked pixels <a href="#values-for-survey-locations-in-masked-pixels" id="values-for-survey-locations-in-masked-pixels"></a>

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.

![](https://2946503769-files.gitbook.io/~/files/v0/b/gitbook-legacy-files/o/assets%2F-LylLNCSXaUER_FiqDSx%2F-MG9oVcTHhvNHsfjOdbS%2F-MG9oyvA2Rp4kvuSJENq%2Fimage.png?alt=media\&token=70aee390-65b9-48f0-b788-51b11e601f61)

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.](https://2946503769-files.gitbook.io/~/files/v0/b/gitbook-legacy-files/o/assets%2F-LylLNCSXaUER_FiqDSx%2F-MG9oVcTHhvNHsfjOdbS%2F-MG9pawIUSzT-wrREmjH%2Fimage.png?alt=media\&token=c9feb346-5774-41d2-9d55-4ba405d4a6c1)

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]) {` \
&#x20;  `d=replace(distanceFromPoints(r, xy[i,]), is.na(r), NA)`\
&#x20;  `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)`

![](https://2946503769-files.gitbook.io/~/files/v0/b/gitbook-legacy-files/o/assets%2F-LylLNCSXaUER_FiqDSx%2F-MGFSP6AlztcdQ_hsAQd%2F-MGFWKQGnW3YeaZIOHSK%2Fimage.png?alt=media\&token=d36526b5-5556-475b-b30c-cdb5677cc5b1)