# 2. Comparison of chlorophyll data from different sensors

This tutorial will showcase the use of the `rerddap` package, which was developed to make it easier to interact with ERDDAP servers from R.\
\
More information about the package can be found here:\
<https://cran.r-project.org/web/packages/rerddap/index.html>\
\
and here: <https://cran.r-project.org/web/packages/rerddap/vignettes/Using_rerddap.html>

As an example, we are going to plot time-series of mean chlorophyll a concentration from various sensors from 1997 to 2018 to look at the periods of overlap. \
We are going to download data from Seawifs (1997-2010), MODIS (2002-present) and VIIRS (2012-present) and compare it to the ESA-CCI data which combines all 3 sensors into a homogeneous time-series.

Load packages:\
`require(lubridate)` \
`require(plyr)` \
`library(rerddap)` \
`library(ncdf4)`

First we define the longitude-latitude boundaries of the box:

`xcoord<-c(198,208)` \
`ycoord<-c(15,25)`

Next we define the URL of the ERDDAP we will be using:

`ERDDAP_Node="`[`https://oceanwatch.pifsc.noaa.gov/erddap/`](https://oceanwatch.pifsc.noaa.gov/erddap/)`"`

## **Get monthly seawifs data, which starts in 1997.**

Go to ERDDAP to find the name of the dataset for monthly SeaWIFS data: `sw_chla_monthly_2018_0`

You should always examine the dataset in ERDDAP to check the date range, names of the variables and dataset ID, to make sure your `griddap` calls are correct:\
<https://oceanwatch.pifsc.noaa.gov/erddap/griddap/sw_chla_monthly_2018_0.html>

First we need to know what our variable is called:

`dataInfo <- rerddap::info('sw_chla_monthly_2018_0', url=ERDDAP_Node)` \
`var=dataInfo$variable$variable_name`\
\
`var`

`[1] "chlor_a" "palette"`

We are interested in the `chlor_a` variable, which contains the values of chlorophyll-a concentration. \
This is `var[1]`.

`griddap` is a function from the `rerddap` package. It grabs the data from ERDDAP based on the parameters we give it.

We are grabbing a lot of data so all the `griddap` commands might take a while.

`sw <- griddap(url=ERDDAP_Node, 'sw_chla_monthly_2018_0', time = c('1997-12-01', '2010-12-01'), latitude = ycoord, longitude = xcoord, fields = var[1] )`

Spatially average all the data within the box

`swAVG <- ddply( sw$data, .(time), function(x) mean(x$chlor_a, na.rm =TRUE) )`

## Get monthly MODIS data, which starts in 2002.

`dataInfo <- rerddap::info('aqua_chla_monthly_2018_0', url=ERDDAP_Node)` \
`var=dataInfo$variable$variable_name`\
\
`MOD <- griddap(url=ERDDAP_Node, 'aqua_chla_monthly_2018_0', time = c('2002-07-16', '2018-12-16'), latitude = ycoord, longitude = xcoord, fields = var[1])`

Spatially average all the data within the box:

`MODAVG <- ddply( MOD$data, .(time), function(x) mean(x$chlor_a, na.rm =TRUE) )`

## Get monthly VIIRS data, which starts in 2012.

`dataInfo <- rerddap::info('noaa_snpp_chla_monthly', url=ERDDAP_Node)` \
`var=dataInfo$variable$variable_name`\
\
`VIIRS <- griddap(url=ERDDAP_Node, 'noaa_snpp_chla_monthly', time = c('2012-01-02', '2018-12-01'), latitude = ycoord, longitude = xcoord, fields = var)`

Spatially average all the data within the box:

`VIIRSAVG <- ddply( VIIRS$data, .(time), function(x) mean(x$chlor_a, na.rm =TRUE) )`

## Get OC-CCI data (September 1997 to Dec 2018)

`dataInfo <- rerddap::info('esa-cci-chla-monthly-v4-2', url=ERDDAP_Node)` \
`var=dataInfo$variable$variable_name`\
\
`CCI <- griddap(url=ERDDAP_Node, 'esa-cci-chla-monthly-v4-2', time = c('1997-09-04', '2018-12-01'), latitude = ycoord, longitude = xcoord, fields = var[1] )`

Spatially average all the data within the box:

`CCIAVG <- ddply( CCI$data, .(time), function(x) mean(x$chlor_a, na.rm =TRUE) )`

## Plot time series result

`plot(as.Date(swAVG$time), swAVG$V1, type='l', col=2,lwd=2, xlab="", xlim=as.Date(c("1997-12-01","2018-12-01")), ylim=c(0.035,0.10), ylab="CHL")` \
`axis(2)` \
`points(as.Date(swAVG$time), swAVG$V1,pch=20,col=2)`\
\
`lines(as.Date(MODAVG$time), MODAVG$V1, col=4, lwd=2)`\
`points(as.Date(MODAVG$time), MODAVG$V1,pch=20,col=4)`\
\
`lines(as.Date(VIIRSAVG$time), VIIRSAVG$V1, col=3, lwd=2)`\
`points(as.Date(VIIRSAVG$time), VIIRSAVG$V1,pch=20,col=3)`\
\
`legend('bottomleft',legend=c('sw','mod','viirs'),cex=0.6,col=c(2,4,3),lwd=2)`

![](https://2946503769-files.gitbook.io/~/files/v0/b/gitbook-legacy-files/o/assets%2F-LylLNCSXaUER_FiqDSx%2F-MK6Rf0eG1-YzhRZmJke%2F-MK6UiXuTr3pbY_7m06S%2Fimage.png?alt=media\&token=853208b8-8e00-4a2e-bc89-67eebf45d0eb)

You can see that the values of chl-a concentration don't match between sensors.

## Make another plot with CCI as well to compare

`plot(as.Date(swAVG$time), swAVG$V1, type='n', col=2,lwd=2, xlab="", xlim=as.Date(c("1997-12-01","2018-12-01")), ylim=c(0.035,0.10), ylab="CHL")` \
`axis(2)` \
\
`points(as.Date(swAVG$time), swAVG$V1,pch=20,col=2)`\
`points(as.Date(MODAVG$time), MODAVG$V1,pch=20,col=4)`\
`points(as.Date(VIIRSAVG$time), VIIRSAVG$V1,pch=20,col=3)`\
`lines(as.Date(CCIAVG$time),CCIAVG$V1)` \
\
`legend('bottomleft',legend=c('sw','mod','viirs','cci'), cex=0.6,col=c(2,4,3,1),pch=c(20,20,20,NA),lty=c(NA,NA,NA,1),lwd=2)`

![](https://2946503769-files.gitbook.io/~/files/v0/b/gitbook-legacy-files/o/assets%2F-LylLNCSXaUER_FiqDSx%2F-MK6Rf0eG1-YzhRZmJke%2F-MK6VhpK63FobHGeDrOz%2Fimage.png?alt=media\&token=2289bd36-d185-456a-b4c1-65ba8ef0b9c7)