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/"

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)

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)