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 :
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()
Last updated
