1. How to work with satellite data in R
Last updated
Last updated
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')
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 :
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))
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')
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()
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()
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()
