3. Extract data within a shapefile using ERDDAP

This tutorial uses the Xtractomatic package to download data from within the boundaries of the Papahanaumokuakea Marine National Monument (PMNM).

This script uses the rerddap version of Xtractomatic: rerddapXtracto.
More information is available here: https://coastwatch.pfeg.noaa.gov/xtracto/
and here: https://cran.r-project.org/web/packages/rerddapXtracto/

This tutorial will teach you how to extract and display SST values for a particular time period or average SST over the whole time-series available within a shapefile.

The shapefile for the PMNM boundaries can be downloaded here: 
http://sanctuaries.noaa.gov/library/imast_gis.html. 
Extract the monument boundaries from the downloaded shapefile.
Save the file https://sanctuaries.noaa.gov/library/imast/pmnm_py.zip on your computer and extract it.
Navigate to the PMNM_py_files folder.
If you don't have the sf package installed on your computer, you will need to install it.
install.packages('sf')
library(sf)

test=st_read('PMNM_py_files/PMNM_py.shp') 
PNMN=st_coordinates(test)
write.csv(PNMN,'pnmn.csv')
Open the pnmn.csv file in Excel and delete rows 494-497 and 2703-2705. This will prevent an extra line from being plotted where the monument boundary crosses the dateline.
Load packages
(use the install.packages('name_of_package') command to install the packages you don't already have).
library(maps) 
library(mapdata) 
library(maptools) 
library(rerddap) 
library(rerddapXtracto)
Load the Monument boundary
We are only interested in column 2 and 3 of pnmn.csv, which correspond to the longitudes and latitudes of the boundary, respectively.
poly=read.csv('pnmn.csv',header=TRUE)[,2:3]
names(poly) <- c("lon","lat")
We need to convert the longitudes to 0-360º to make it easier to work across the dateline.
I=which(poly$lon<0)
poly$lon[I]=poly$lon[I]+360
Data extraction
The example below extracts monthly 5km geopolar blended SST data within the monument boundary.
The time variable passed to xtractogon must be two elements, the start and endpoints of the desired time period.
ERDDAP_Node="https://oceanwatch.pifsc.noaa.gov/erddap/"
dataInfo <- rerddap::info('goes-poes-monthly-ghrsst-RAN', url=ERDDAP_Node)
Let's look at the variables in our dataset:
dataInfo$variable$variable_name
[1] "analysed_sst" "sea_ice_fraction"
We are interested in sea surface temperature data, so we need to select the fist variable:
analysed_sst

parameter=dataInfo$variable$variable_name[1]

xcoord <- poly$lon 
ycoord <- poly$lat
Let's download data for a few months:
tcoord <- c("2015-03-15", "2015-11-16")
Let's extract the data within our polygon.
Note! The rxtractogon function will not work for regions across the dateline (like PNMN...) and datasets with longitudes between -180 and 180º. For regions that cross the dateline, you need to find a dataset on ERDDAP with longitudes between 0 and 360º.
This command can take a while to run:
sst <- rxtractogon (dataInfo, parameter=parameter, xcoord=xcoord, ycoord=ycoord, tcoord=tcoord)
Plotting the data
The extracted data contains several time steps (months) of sst data in the monument boundaries. 
Let's make a plot of the second time step for example:
i=2
We need to transform our extracted SST data for that time step into a matrix for the image function to work properly. This SST dataset is also in Kelvin degrees, so let-s convert it to Celsius degree:
test=as.matrix(sst$analysed_sst[,,i])-273.15
Let's get the date for our time step:
date=as.Date(strsplit(as.character(sst$time[i]),' ')[[1]][1])
Let's set the bounds of our region: 
xlim = c(177,207) 
ylim = c(18,32)
Let's define a custom color scale:
jet.colors <-colorRampPalette(c("blue", "#007FFF", "cyan","#7FFF7F", "yellow", "#FF7F00", "red", "#7F0000")) 
h=hist(test, 50, plot=FALSE) 
breaks=h$breaks 
n=length(breaks)-1 
c=jet.colors(n)
Make sure your graphics window is pretty wide or you might get an error message about the plot margins.
The layout function splits the graphics window. One part for the map, the other for the color scale.
layout(matrix(c(1,2,3,0,4,0), nrow=1, ncol=2), widths=c(5,1), heights=4) 
layout.show(2)
Let's plot a basemap and overlay the Monument boundary:
par(mar=c(3,3,3,1)) 
land <- maps::map("world2", fill = TRUE, col='grey',xlim = xlim, ylim = ylim, plot=FALSE, xaxs='i',yaxs='i',asp=1, cex.axis = 3, proj4string = CRS("+proj=longlat +datum=WGS84"))

ids <- sapply(strsplit(land$names, ":"), function(x) x[1]) 
bPols <- map2SpatialPolygons(land, IDs = ids, proj4string = CRS("+proj=longlat +datum=WGS84")) 

plot(bPols, col = "grey", axes = FALSE, xlim = xlim, ylim = ylim, xaxs='i',yaxs='i',asp=1, cex.axis = 3, main=format(date,'%b %Y')) 
lines(poly$lon,poly$lat,lwd=2)
Let's make pretty axes:
x = seq(175, 205, 5) 
axis(1, x, sapply(paste(-x + 360, "ºW", sep = ""), function(x) x[1])) 

y = seq(10, 45, 5) 
axis(2, las = 1, y, sapply(paste(y, "ºN", sep = ""), function(x) x[1])) 

box()
We'll add the color plot on top of the map:
par(new=TRUE)
image(sst$longitude,sst$latitude,test,col=c,breaks=breaks,xlab='',ylab='', axes=FALSE, xlim = xlim, ylim = ylim,xaxs='i',yaxs='i',asp=1, las=1)
Plot color scale using the scale.R file:
source('scale.R')
par(mar=c(3,1,3,3))
image.scale(test, col=c, breaks=breaks, horiz=FALSE, yaxt="n",xlab='',ylab='',main='SST') 
axis(4,las=1) 
box()
On your own!
Plot the average SST for the period we downloaded:
Hint: here is how to compute the average SST
mean_sst=as.matrix(apply(sst$analysed_sst,c(1,2),mean,na.rm=TRUE))
Now use mean_sst to plot instead of test.
2. Comparison of chlorophyll data from different sensors
4. Extract data along a turtle track
Last modified 2yr ago