CoastWatch Satellite Course
  • Preface
    • CoastWatch Learning Portal
  • LECTURES
    • Introduction
    • Remote Sensing Basics
    • Ocean Color
    • Sea Surface Temperature
    • Altimetry
    • Ocean Surface Winds
    • Sea Surface Salinity
    • What Dataset to Choose?
    • Applications of Satellite Data
  • Tutorials
    • NetCDF and Panoply tutorial
    • ERDDAP tutorial
    • R tutorial
      • 1. How to work with satellite data in R
      • 1. How to work with satellite data in R - Great Lakes example
      • 2. Comparison of chlorophyll data from different sensors
      • 3. Extract data within a shapefile using ERDDAP
      • 4. Extract data along a turtle track
      • 5. Extract ocean color data in optically-shallow waters - new version
    • Python tutorial
      • 1. How to work with satellite data in Python
      • 2. Comparison of chlorophyll data from different sensors
      • 3. Extract data within a shapefile using ERDDAP
      • 4. Extract data along a turtle track
    • ArcGIS tutorial
      • Additional ArcGIS Tutorials
  • Archive
    • Remote Sensing Basics - Shorter version
    • 5. Extract ocean color data in optically-shallow waters
    • Voyager tutorial
    • index
Powered by GitBook
On this page
  • Extract the monument boundaries from the downloaded shapefile.
  • Load packages
  • Load the Monument boundary
  • Data extraction
  • Plotting the data
  • On your own!

Was this helpful?

  1. Tutorials
  2. R tutorial

3. Extract data within a shapefile using ERDDAP

Previous2. Comparison of chlorophyll data from different sensorsNext4. Extract data along a turtle track

Last updated 4 years ago

Was this helpful?

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: and here: 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: .

Extract the monument boundaries from the downloaded shapefile.

Save the file 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

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)

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.

The example below extracts data within the monument boundary. The time variable passed to xtractogon must be two elements, the start and endpoints of the desired time period.

Plot color scale using the file: source('scale.R')

https://coastwatch.pfeg.noaa.gov/xtracto/
https://cran.r-project.org/web/packages/rerddapXtracto/
http://sanctuaries.noaa.gov/library/imast_gis.html
https://sanctuaries.noaa.gov/library/imast/pmnm_py.zip
monthly 5km geopolar blended SST
scale.R