library(ggplot2)
library(mapdata)
## Loading required package: maps
library(plotdap)
library(plotly)
## Registered S3 method overwritten by 'httr':
##   method           from  
##   print.cache_info hoardr
## 
## Attaching package: 'plotly'
## The following object is masked from 'package:ggplot2':
## 
##     last_plot
## The following object is masked from 'package:stats':
## 
##     filter
## The following object is masked from 'package:graphics':
## 
##     layout
library(rerddap)
library(rerddapXtracto)

I was asked at my seminar on using ‘R’ to extract data from ERDDAP servers if it was possible to make the ‘plotdap’ related maps interactve. It turns out it can be done readily using the ‘plotly’ package without having to add anything to ‘plotdap’, by understanding the structure of a ‘plotdap’ object, and if needed using the bbox_set() function in ‘plotdap’, and a trick I will show on how to switch the order of layers in a ggplot object.

Be warned that the ‘plotly’ rendering can be quite slow.

Let’s look at one of the examples in the ‘plotdap’ vignette.

sstInfo <- rerddap::info('erdVHsstaWS3day')
# get latest 3-day composite sst
viirsSST <- rerddap::griddap(sstInfo, 
                             latitude = c(41., 31.), 
                             longitude = c(-128., -115), 
                             time = c('last','last'), 
                             fields = 'sst')
 
w <- map("worldHires", xlim = c(-140., -114), ylim = c(30., 42.), fill = TRUE, plot = FALSE)
my_plot <- add_griddap(plotdap(mapData = w), viirsSST, ~sst, fill = "thermal" )

my_plot 

Part of a ‘plotdap’ object is a ‘ggplot2’ object as well as information on bounds and projections. For now we will finesse the question of projections.

The ‘plotly’ package has a function ‘ggplotly()’ that can workj directly with a ‘ggplot2’ object to produce an interactive graphic. It will work with the plot we just made, but perhaps with an unexpected result.

ggplotly(my_plot$ggplot)

The resulting plot is interactive but the bounding box is off. This is because ‘plotdap’ has it’s own print method that calculates a good bounding box before displaying the map. But the function bbox_set() in ‘plotdap’ will set the bounding box in the ggplot object.

my_plot1 <- bbox_set(my_plot, c(-128., -115), c(31, 41.))
ggplotly(my_plot1$ggplot)

Now the bounding box is correct. the ‘plotdap’ print method is also what is used to have the land overlay the data grid, rather than vice-versa. But the code from that method cna be extracted and used here for the same purpose.

my_plot2 <- my_plot1
gg <- my_plot2$ggplot
gg$layers <- rev(gg$layers)
my_plot2$ggplot <- gg
ggplotly(my_plot2$ggplot)

The same can be done with output from the ‘rerddapXtracto’ functions plotTrack() and plotBBox(), as they just produce ‘plotdap’ objects. Here it is done for the Marlin track example, including setting a good bounding box.

tagData <- Marlintag38606
xpos <- tagData$lon
ypos <- tagData$lat
tpos <- tagData$date
zpos <- rep(0., length(xpos))
swchlInfo <- rerddap::info('erdSWchla8day')
swchl1 <- rxtracto(swchlInfo, parameter = 'chlorophyll', xcoord = xpos, ycoord = ypos, tcoord = tpos, zcoord = zpos, xlen = .2, ylen = .2)
#
myPlot <- plotTrack(swchl1, xpos, ypos, tpos, plotColor = 'algae')
myPlot <- bbox_set(myPlot, c(min(xpos - 360), max(xpos - 360)), c(min(ypos), max(ypos)))
ggplotly(myPlot$ggplot)