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.htmlAs 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:
Go to ERDDAP to find the name of the dataset for monthly SeaWIFS data:
sw_chla_monthly_2018_0You 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.htmlFirst 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) )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) )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) )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(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.
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)