1. How to work with satellite data in R - Great Lakes example

This tutorial will show the steps to grab data in ERDDAP from R, how to work with NetCDF files in R and how to make some maps and time-series of sea surface temperature (SST) around the Great Lakes.

If you do not have the ncdf4 and httr packages installed in R, you will need to install them:

install.packages('ncdf4') install.packages('httr')

Downloading data in R

Because ERDDAP includes RESTful services, you can download data listed on any ERDDAP platform from R using the URL structure.

For example, the following page allows you to subset daily SST data:

Select your region and date range of interest, then select the '.nc' (NetCDF) file type and click on "Just Generate the URL".

In this specific example, the URL we generated is :

​https://coastwatch.glerl.noaa.gov/erddap/griddap/GLSEA_GCS.htmlTable?sst[(2021-07-21T12:00:00Z):1:(2021-07-28T12:00:00Z)][(38.8749871947229):1:(50.6059751976437)][(-92.4199507342304):1:(-75.8816402880577)]​

You can also edit this URL manually. In R, run the following to download the data using the generated URL (you need to copy it from your browser):

library(ncdf4) library(httr)

junk <- GET('https://oceanwatch.pifsc.noaa.gov/erddap/griddap/CRW_sst_v1_0_monthly.nc?analysed_sst[(2018-01-01T12:00:00Z):1:(2018-12-01T12:00:00Z)][(17):1:(30)][(195):1:(210)]', write_disk("sst.nc", overwrite=TRUE))

Importing the downloaded data in R

Now that we've downloaded the data locally, we can import it and extract our variables of interest:

• open the file

nc=nc_open('sst.nc')

• examine which variables are included in the dataset:

names(nc$var)

[1] "sst"

• Extract sst:

v1=nc$var[[1]] sst=ncvar_get(nc,v1)

• examine the structure of sst:

dim(sst)

[1] 1181 838 7

Our dataset is a 3-D array with 1181 rows corresponding to longitudes, 838 columns corresponding to latitudes for each of the 7 time steps that we downloaded.

• get the dates for each time step:

dates=as.POSIXlt(v1$dim[[3]]$vals,origin='1970-01-01',tz='GMT') dates

[1] "2021-07-21 12:00:00 GMT" "2021-07-22 12:00:00 GMT" [3] "2021-07-23 12:00:00 GMT" "2021-07-24 12:00:00 GMT" [5] "2021-07-25 12:00:00 GMT" "2021-07-26 12:00:00 GMT" [7] "2021-07-27 12:00:00 GMT"

• get the longitude and latitude values

lon=v1$dim[[1]]$vals lat=v1$dim[[2]]$vals

• Close the netcdf file and remove the data and files that are not needed anymore.

nc_close(nc) rm(junk,v1) file.remove('sst.nc')

Creating a map for one time step

Let's create a map of SST for 07/21/2021 (our first time step). You will need to download the scale.R file and copy it to your working directory to plot the color scale properly.

• set some color breaks

breaks=seq(10,25.5,0.05) n=length(breaks)-1

• define a color palette

jet.colors <-colorRampPalette(c("blue", "#007FFF", "cyan","#7FFF7F", "yellow", "#FF7F00", "red", "#7F0000"))

• set color scale using the jet.colors palette

c=jet.colors(n)

• prepare graphic window : left side for map, right side for color scale

layout(matrix(c(1,2,3,0,4,0), nrow=1, ncol=2), widths=c(5,1), heights=4) layout.show(2) par(mar=c(3,3,3,1))

• plot the SST map

image(lon,lat,sst[,,1],col=c,breaks=breaks,xlab='',ylab='',axes=TRUE,xaxs='i',yaxs='i',asp=1, main=paste("Monthly SST", dates[1]))

• example of how to add points to the map

points(seq(-83,-81,0.5),rep(45,5), pch=20, cex=1)

• example of how to add a contour (this is considered a new plot, not a feature, so you need to use par(new=TRUE)) to overlay it on top of the SST map

par(new=TRUE) contour(lon,lat,sst[,,1],levels=15,xaxs='i',yaxs='i',labcex=0.8,vfont = c("sans serif", "bold"),axes=FALSE,asp=1)

• plot color scale using 'image.scale' function from 'scale.R' script)

par(mar=c(3,1,3,3)) source('scale.R') image.scale(sst[,,1], col=c, breaks=breaks, horiz=FALSE, yaxt="n",xlab='',ylab='',main='SST') axis(4, las=1) box()

Plotting a time series

Let's pick a box encompassing Lake Superior. We are going to generate a time series of mean SST within that box.

I=which(lon<=-84) J=which(lat>=46) sst2=sst[I,J,] n=dim(sst2)[3] res=rep(NA,n) for (i in 1:n) res[i]=mean(sst2[,,i],na.rm=TRUE) plot(1:n,res,axes=FALSE,type='o',pch=20,xlab='',ylab='SST (ºC)') axis(2) axis(1,1:n,format(dates,'%m/%d')) box()

Creating a map of average SST over a week

sst.mean=apply(sst[,,1:7],c(1,2),mean,na.rm=TRUE) breaks=seq(10,25.5,0.05) n=length(breaks)-1 c=jet.colors(n) layout(matrix(c(1,2,3,0,4,0), nrow=1, ncol=2), widths=c(5,1), heights=4) layout.show(2) par(mar=c(3,3,3,1)) image(lon,lat,sst.mean,col=c,breaks=breaks,xlab='',ylab='',axes=TRUE,xaxs='i',yaxs='i',asp=1,main=paste("Mean SST", format(dates[1],'%Y/%m/%d'),' - ',format(dates[7],'%Y/%m/%d'))) par(mar=c(3,1,3,3)) image.scale(sst.mean, col=c, breaks=breaks, horiz=FALSE, yaxt="n",xlab='',ylab='',main='SST') axis(4) box()

​