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
  • Downloading data in R
  • Importing the downloaded data in R
  • Working with the extracted data
  • Creating a map for one time step
  • Plotting a time series
  • Creating a map of average SST over a year

Was this helpful?

  1. Tutorials
  2. R tutorial

1. How to work with satellite data in R

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

Last updated 3 years ago

Was this helpful?

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 :

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)

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

  • 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()

junk <- GET('', write_disk("sst.nc", overwrite=TRUE))

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

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)]
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)]
scale.R