1. How to work with satellite data in R

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 main Hawaiian islands.

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 monthly 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://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)]

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] "analysed_sst"

• Extract analysed_sst:

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

• examine the structure of sst:

dim(sst)

[1] 301 261 12

Our dataset is a 3-D array with 301 rows corresponding to longitudes, 261 columns corresponding to latitudes for each of the 12 time steps.

• get the dates for each time step:

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

[1] "2018-01-01 GMT" "2018-02-01 GMT" "2018-03-01 GMT" "2018-04-01 GMT" "2018-05-01 GMT" "2018-06-01 GMT" "2018-07-01 GMT" [8] "2018-08-01 GMT" "2018-09-01 GMT" "2018-10-01 GMT" "2018-11-01 GMT" "2018-12-01 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')

Working with the extracted data

Creating a map for one time step

Let's create a map of SST for January 2018 (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

h=hist(sst[,,1], 100, plot=FALSE) breaks=h$breaks 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(202:205,rep(26,4), pch=20, cex=2)

• 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=20,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 the following box : 24-26N, 200-206E. We are going to generate a time series of mean SST within that box.

I=which(lon>=200 & lon<=206) J=which(lat>=24 & lat<=26) 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')) box()

Creating a map of average SST over a year

sst.yr=apply(sst[,,1:12],c(1,2),mean,na.rm=TRUE) h=hist(sst.yr, 100, plot=FALSE) breaks=h$breaks 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.yr,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[12],'%Y/%m/%d'))) par(mar=c(3,1,3,3)) image.scale(sst.yr, col=c, breaks=breaks, horiz=FALSE, yaxt="n",xlab='',ylab='',main='SST') axis(4) box()