reprojecting_satellite_buoy_data.Rmd

---
title: Reprojecting satellite and buoy data
author: NOAA CoastWatch
date: May 20, 2021
output:
  md_document:
    variant: gfm
---

# Reprojecting satellite and buoy data 
>notebook filename | reprojecting_satellite_buoy_data.Rmd  
history | Created Mar 2020 | Updated Jan 2022

The CoastWatch program provides ocean satellite data and also non-satellite data that may be useful for ocean scientists. In exercise, we will obtain cruise track data from a research ship, in-situ buoy data, and satellite sea surface temperature data from the PolarWatch ERDDAP data server and then plot all three types of the data on a map in an Alaska Albers projection.  

The exercise demonstrates the following techniques:  

* Finding and accessing in-situ buoy data in ERDDAP  
* Accessing satellite data from ERDDAP  
* Making projected maps of ship tracks, buoy locations and satellite data  

Key R packages/functions used:

* **rerddap::info** - to retrieve information about a dataset from ERDDAP   
* **rerddapXtracto::xtracto_3D** - to extract satellite data from ERDDAP  
* **rerddap::tabledap** - to extract ship track and buoy data from ERDDAP    
* **oce** - to plot data from ERDDAP on a projected map 

Datasets used:  

* NOAA Geo-polar Blended Sea Surface Temperature daily satellite data, 'nesdisGeoPolarSSTN5NRT'   
* NDBC Buoys, 'cwwcNDBCMet'  
* A Ship Track from the R/V Healy 'fsuResearchShipNEPP' 

These datasets are all provided in ERDDAP with lat-lon coordinates. For an example of working with datasets provided in projected (ie. polar stereographic) coordinates see the "Mapping projected datasets" exercise in this tutorial.

## Install required packages and load libraries  

```{r install, message=FALSE, warning=FALSE}
# Function to check if pkgs are installed, install missing pkgs, and load
pkgTest <- function(x)
{
  if (!require(x,character.only = TRUE))
  {
    install.packages(x,dep=TRUE)
    if(!require(x,character.only = TRUE)) stop(x, " :Package not found")
  }
}

list.of.packages <- c( "ncdf4","parsedate", "rerddap", "plotdap", "sp", "httr",
                       "lubridate", "gridGraphics",  "mapdata", "cmocean",
                       "ggplot2", "RColorBrewer", "grid", "PBSmapping", "oce", 
                       "ocedata", "rerddapXtracto")
# create list of installed packages
pkges = installed.packages()[,"Package"]
for (pk in list.of.packages) {
  pkgTest(pk)
}
```

## Ship Track
Get a track from the research ship Healy in Alaska. This dataset is in the PolarWatch ERDDAP with dataset id 'fsuResearchShipNEPP'. Here we pull the track for a segment of data from June 2011.

```{r}
# View variable names to help form the track request
shipDataInfo <- rerddap::info('fsuResearchShipNEPP')

# Use details reurned from the info request to form a request for the track
ship.df <- tabledap(shipDataInfo, fields = c('latitude',  'longitude', 'time','airTemperature'), 'time>=2011-06-21T00:00:00Z', 'time<=2011-06-29T00:00:00Z')

```
__Visualize the ship track__   

Here we demonstrate viewing the track in Alaska Albers projection which is commonly used by local scientists. See the Fish and Wildlife reference link for more info on projections in Alaska.

The **oce** package is great for making maps in a variety of projections. The **oce::mapPlot** function accepts projections in a number of formats, see the mapPlot reference below for more details. Here we specify the projection with the proj4 string for Alaska Albers EPSG:6393. 

```{r plot1}
# define map extents
ylim <- c(53,69)
xlim <- c(-179,-130)
data("coastlineWorldMedium") # included in ocedata

# set plot margins
par(mar=c(2, 2, 3, 1))

# define projection with proj4 string for Alaska Albers EPSG:6393
epsg6393 = '+proj=aea +lat_1=55 +lat_2=65 +lat_0=50 +lon_0=-154 +x_0=0 +y_0=0 +ellps=GRS80 +units=m +no_defs' 

## make a base
mapPlot(coastlineWorldMedium, 
        projection=epsg6393,
        col="lightgray", 
        longitudelim=xlim, 
        latitudelim=ylim, 
        clip = TRUE,
        main="         Ship Track  \n   Alaska Albers Projection      "
        )

## add ship track
mapLines(longitude = as.numeric(ship.df$longitude), 
         latitude = as.numeric(ship.df$latitude), 
         col = "black", lwd = 4)
```

## Add Nearby Buoy locations
See if any NDBC buoys have data available in this area and add them to the map. Request NDBC buoy data with the **rerddap::tabledap** function, using the same temporal and spatial coverages as our satellite data request. Here we request the station, latitude, longitude, and time and then put the data into a data frame.

```{r getBuoys}

buoysDatasetId <- 'cwwcNDBCMet'

# Request the NDBC buoy data
buoys <- tabledap(
  buoysDatasetId, 
  fields=c('station', 'latitude',  'longitude', 'time', 'wtmp'), 
  'time>=2011-06-21',   'time<=2011-06-29', 
  'latitude>=50','latitude<=70', 
  'longitude>=-179','longitude<=-150'
)
buoys$wtmp <- as.numeric(buoys$wtmp)

# Create dataframe 
buoys.df <-data.frame(
  station=buoys$station,
  longitude=as.numeric(buoys$longitude),
  latitude=as.numeric(buoys$latitude),
  time=strptime(buoys$time, "%Y-%m-%dT%H:%M:%S"),
  sst=as.numeric(buoys$wtmp),
  col=as.numeric((buoys$wtmp-min(buoys$wtmp,na.rm=T))/max(buoys$wtmp,na.rm=T))
  )
# Subset buoys to only show those with data at night, to correspond with the satellite data. 
# Night hours are 2000 - 0700 local = (0400 - 1500 GMT)
nightbuoy.df <- subset(buoys.df,hour(time)==8)
#nightbuoy.df <- subset(nightbuoy.df,hour(time)<15)

#set sst colors to scale of 0 to 1
nightbuoy.df$col[nightbuoy.df$col <= 0] = 0
nightbuoy.df$col[nightbuoy.df$col > 1] = 1

# See how many buoys were found
length(unique(nightbuoy.df$station))

```

**Make a map with both the ship track and buoy locations**
```{r plot2}

## draw scalebar
Zlim = c(0,12)

drawPalette(zlim = Zlim, 
            breaks = seq(0,12,2),
            col= oceColorsTemperature, 
            at = seq(0,12,2), 
            pos = 4, 
            drawTriangles = FALSE)

## make a base
mapPlot(coastlineWorldMedium, 
        projection=epsg6393,
        col="lightgray", 
        longitudelim=xlim, 
        latitudelim=ylim, 
        clip = TRUE,
        main="Healy Track and NDBC Buoys \n   Alaska Albers Projection      "
        )

## add ship track
mapLines(longitude = as.numeric(ship.df$longitude), 
         latitude = as.numeric(ship.df$latitude), 
         col = "black", lwd = 4)

numColors <- 6 
mypalette <- rev(cmocean('thermal')(numColors))

## add the locations of the NDBC buoys

# mapPoints(longitude = buoys.df$longitude, 
#           latitude = buoys.df$latitude, 
#           col = "#0099cc", pch = 20, cex = 1.25)

mapPoints(longitude = nightbuoy.df$longitude, 
          latitude = nightbuoy.df$latitude, 
          col = mypalette[as.integer(1+(numColors-1)*nightbuoy.df$col)], 
          pch = 20, cex = 1.25)

```

## Add Satellite SST Data

* Request metadata for the Geo-Polar Blended SST dataset using the **rerddap::info** function. 
* Use the returned info about dimensions and variable names to form a data request.
* The satellite dataset id is **nesdisGeoPolarSSTN5NRT**

```{r satdataInfo}

# Request and view info about satellite dataset
WCN_ERDDAP <- "https://coastwatch.pfeg.noaa.gov/erddap/"
satdataInfo <- rerddap::info('nesdisGeoPolarSSTN5SQNRT', url = WCN_ERDDAP)
satdataInfo
```

**Extract Satellite SST Data**

Using information returned from the **rerddap::info** request, form a data request using **rerddapXtracto::rxtracto_3D**. Here we request SST data for June 24, 2011 using the same extents that we used to define our map.

```{r getSatelliteGrid}

# Set the parameter and spatial extent for the satellite data request 
parameter='analysed_sst'

satMap <- rxtracto_3D(satdataInfo, xcoord=xlim, ycoord=ylim, parameter=parameter, tcoord=c('2011-06-24','2011-06-24'))

satMap$sat_sst <- drop(satMap$sat_sst)

```

**Make Map with Ship Track, Buoys and Satellite SST**

```{r plot3}

# set plot margins
par(mar=c(2, 2, 4, 1))

## draw scalebar
Zlim = c(0,12)

drawPalette(zlim = Zlim, 
            breaks = seq(0,12,2),
            col= oceColorsTemperature, 
            at = seq(0,12,2), 
            pos = 4, 
            drawTriangles = FALSE)

## make a base
mapPlot(coastlineWorldMedium, 
        projection=epsg6393,
        col="lightgray", 
        longitudelim=xlim, 
        latitudelim=ylim, 
        clip = TRUE,
        main="   Healy Ship Track, NDBC Buoys and Satellite SST \n   Alaska Albers Projection      "
        )

## overlay gridded sst layer
## The oce package incorporates cmocean colormaps
mapImage(longitude = satMap$longitude,
         latitude = satMap$latitude, 
         z = satMap$analysed_sst[,,1], 
         zlim = Zlim,
         col = oceColorsTemperature)

## add a land polygon on top of the data
mapPolygon(coastlineWorldMedium, 
        col="lightgray")

## add ship track
mapLines(longitude = as.numeric(ship.df$longitude), 
         latitude = as.numeric(ship.df$latitude), 
         col = "black", lwd = 4)

# create points colormap
numColors <- 6 
pointsPalette <- rev(cmocean('thermal')(numColors))

mapPoints(longitude = nightbuoy.df$longitude, 
          latitude = nightbuoy.df$latitude, 
          col = pointsPalette[as.integer(1+(numColors-1)*nightbuoy.df$col)], 
          pch = 20, 
          cex = 1.25)

```



## Related Materials 

[Projected Datasets R Notebook](https://coastwatch.pfeg.noaa.gov/projects/r/Projected.html) - Access and plot polar stereographic ice datasets from ERDDAP

[ERDDAP Training Course Book](https://coastwatch.pfeg.noaa.gov/projects/erddap/) - Using the features of ERDDAP

## References
Dan Kelley, OCE Vignettes - Using Map Projections,Feb 21, 2020, https://cran.r-project.org/web/packages/oce/vignettes/map_projections.html

Dan Kelley, R Documentation from oce v1.2.0, mapPlot function, accessed March 11, 2020, https://www.rdocumentation.org/packages/oce/versions/1.2-0/topics/mapPlot

Masumbuko Semba, Mapping with oce, 2/13/2019,  https://semba-blog.netlify.com/02/13/2019/mapping-with-oce/

U.S. Fish and Wildlife Service Alaska Region, A Geodesy Primer, November 08, 2019, https://www.fws.gov/r7/nwr/Realty/data/LandMappers/Public/Help/HTML/R7-Public-Land-Mapper-Help.html?Datumsprojectionsandcoordinatesy.html