---
title: Satellite SST and CHL data extractions from selected locations - COESSING Ghana
  School
output: 
  html_document:
    df_print: paged
---

## Marine Biodiversity Observation Network Pole to Pole of the Americas ([MBON Pole to Pole](https://marinebon.org/p2p/))

Written by E. Montes (emontesh@usf.edu) and Eduardo Klein (eklein@usb.ve) on July 12, 2021.

This code pulls data from NOAA's [ERDDAP](https://coastwatch.pfeg.noaa.gov/erddap/index.html) servers and creates time series plots of sea surface temperature (SST) and chlorophyll-a concentration (CHL), and maps showing the latest available data for the selected region.

# Step 1
First, let's load required libraries
```{r, message=FALSE, warning=FALSE}
library(readr)
library(rerddap)
library(lubridate)
library(dplyr)
library(flexdashboard)
library(reshape2)
library(leaflet)
library(ggplot2)
library(vegan)
library(xts)
library(dygraphs)
library(plotly)
library(mapdata)

library(RColorBrewer)
palette(brewer.pal(8, "Set2"))
```

# Step 2
Query sea surface temperature (SST) data from [ERDDAP](https://coastwatch.pfeg.noaa.gov/erddap/griddap/jplMURSST41mday.html). These are monthly SST means at ~ 1km pixel resolution starting in 2002.
```{r, message=FALSE, warning=FALSE}
## remove all spaces from string
NoSpaces = function(x){
  return(gsub(" ", "", x))
}

## set site coordinates and time for SST extraction. 
SSTSiteName = "Ghana Shelf"   ## for the resulting file name
SSTcoords.lon = -1.2  ## this is decimal longitude
SSTcoords.lat = 4.5  ## this is decimal latitude

SSTstartDate = "2015-01-01"  ## define start date of your time series 

## set dataset source (monthly SST)
SSTsource = info("jplMURSST41mday")

##
## Get sst 
SST <- griddap(SSTsource, 
              time=c(SSTstartDate, "last"),
              longitude = c(SSTcoords.lon,SSTcoords.lon),
              latitude = c(SSTcoords.lat,SSTcoords.lat),
              fields = "sst",
              fmt = "csv")

SST = SST[,c(1,4)]
names(SST) = c("time", "SST")

## convert time to a Data object
SST$time = as.Date(ymd_hms(SST$time))

```

# Step 3
Plot SST time series
```{r, message=FALSE, warning=FALSE}
SST.xts = as.xts(SST$SST, SST$time)
dygraph(SST.xts, 
        ylab = "Sea Surface Temperature (Deg C)") %>% 
  dySeries("V1", label ="SST (Deg C)", color = "steelblue") %>%
  dyHighlight(highlightCircleSize = 5, 
              highlightSeriesBackgroundAlpha = 0.2,
              hideOnMouseOut = FALSE) %>% 
  dyOptions(fillGraph = FALSE, fillAlpha = 0.4) %>% 
  dyRangeSelector(dateWindow = c(max(SST$time) - years(5), max(SST$time)))

```

# Step 4
Save SST time series data as a CSV table
```{r, message=FALSE, warning=FALSE}
write_csv(SST, path = paste0(NoSpaces(SSTSiteName), "_SST.csv"))
```

# Step 5
Create a map of the latest SST data. The blue dot indicates the location of the site where the time series was extracted in the code chunck above.
```{r, message=FALSE, warning=FALSE}
# define geographic domain
lat_lims<- c(-5., 10.)  ## change the first and second value for desired minimum and maximum latitude.
lon_lims <- c(-8., 12.)  ## change the first and second value for desired minimum and maximum longitude.

sstInfo <- info('jplMURSST41mday')
# get latest composite sst
GHRSST <- griddap(sstInfo, latitude = lat_lims, longitude = lon_lims, time = c('last','last'), fields = 'sst')

mycolor <- colors$temperature
w <- map_data("worldHires", ylim = lat_lims, xlim = lon_lims)
ggplot(data = GHRSST$data, aes(x = lon, y = lat, fill = sst)) + 
  geom_polygon(data = w, aes(x = long, y = lat, group = group), fill = "grey80") +
  geom_raster(interpolate = FALSE) +
  geom_point(x = SSTcoords.lon, y = SSTcoords.lat, colour = "blue", size = 3) +
  scale_fill_gradientn(colours = mycolor, na.value = NA) +
  theme_bw() + ylab("latitude") + xlab("longitude") +
  coord_fixed(1.1, xlim = lon_lims,  ylim = lat_lims) + ggtitle("Latest SST")
```


# Step 6
Now, let's query surface chlorophyll-a concentration (CHL; proxy for phytoplankton biomass) data from [ERDDAP](https://coastwatch.pfeg.noaa.gov/erddap/griddap/erdMH1chla8day.html).  These are weekly CHL means from the MODIS Aqua sensor at ~ 4km pixel resolution starting in 2003.
```{r, message=FALSE, warning=FALSE}
## remove all spaces from string
NoSpaces = function(x){
  return(gsub(" ", "", x))
}

## set site coordinates and time for CHL extraction
CHLSiteName = "Gulf of Guinea"   ## for the resulting file name
CHLcoords.lon = -1.2 ## this is decimal longitude
CHLcoords.lat = 4.5 ## this is decimal latitude

CHLstartDate = "2015-01-01"

## set dataset source
CHLsource = info("erdMH1chla8day")

##
## Get CHL 
CHL <- griddap(CHLsource, 
               time=c(CHLstartDate, "last"),
               longitude = c(CHLcoords.lon,CHLcoords.lon),
               latitude = c(CHLcoords.lat,CHLcoords.lat),
               fields = "chlorophyll", fmt = "csv")

CHL = CHL[,c(1,4)]
names(CHL) = c("time", "CHL")
CHL = na.omit(CHL)

## convert time to a Data object
CHL$time = as.Date(ymd_hms(CHL$time))

```

# Step 7
Plot CHL time series
```{r, message=FALSE, warning=FALSE}
CHL.xts = as.xts(CHL$CHL, CHL$time)
dygraph(CHL.xts, 
        ylab = "Chlorophyll a (mg m-3)") %>% 
  dySeries("V1", label ="CHL", color = "steelblue") %>%
  dyHighlight(highlightCircleSize = 5, 
              highlightSeriesBackgroundAlpha = 0.2,
              hideOnMouseOut = FALSE) %>% 
  dyOptions(fillGraph = FALSE, fillAlpha = 0.4) %>% 
  dyRangeSelector(dateWindow = c(max(CHL$time) - years(5), max(CHL$time)))

```

# Step 8
Save CHL time series data
```{r, message=FALSE, warning=FALSE}
write_csv(CHL, path = paste0(NoSpaces(CHLSiteName), "_CHL.csv"))
```

# Step 9
Create a map of the latest CHL data. The blue dot indicates the location of the site where the time series was extracted in the code chunck above.
```{r, message=FALSE, warning=FALSE}
require("rerddap")
require("ggplot2")
require("mapdata")

# define geographic domain
lat_lims<- c(-5., 10.)  ## change the first and second value for desired minimum and maximum latitude.
lon_lims <- c(-8., 12.)  ## change the first and second value for desired minimum and maximum longitude.

chlaInfo <- info('nesdisVHNSQchlaMonthly')
CHLA <- griddap(chlaInfo, latitude = lat_lims, longitude = lon_lims, time = c('last','last'), fields = 'chlor_a')

# Map monthly chl (VIIRS)
mycolor <- colors$chlorophyll
w <- map_data("worldHires", ylim = lat_lims, xlim = lon_lims)
ggplot(data = CHLA$data, aes(x = lon, y = lat, fill = log(chlor_a))) + 
  geom_polygon(data = w, aes(x = long, y = lat, group = group), fill = "grey80") +
  geom_raster(interpolate = FALSE) +
  geom_point(x = CHLcoords.lon, y = CHLcoords.lat, colour = "blue", size = 3) +
  scale_fill_gradientn(colours = mycolor, na.value = NA) +
  theme_bw() + ylab("latitude") + xlab("longitude") +
  coord_fixed(1.3, xlim = lon_lims,  ylim = lat_lims) + ggtitle("Last month")
```
