How to plot AIMS eReefs data with R

This tutorial is based on the article How to open and work with NetCDF data in R.

You can download the source code for this example from the GitHub repository. If you are interested in accessing the raw eReefs model data you can use Barbara Robson's R package eReefs (see R package to facilitate extracting and visualising eReefs model output data for more information).

Load the required R packages

First of all, load the necessary libraries.

library(ncdf4) # package for netcdf manipulation
library(raster) # package for raster manipulation
library(ggplot2) # package for plotting
library(rgdal) # package for geospatial analysis

Define which data to plot

In this section we define which data we want to read and plot.

• inputFile
  The netCDF input file. This can either be a downloaded file (see How to manually download derived data from THREDDS) or a OPeNDAP URL from the AIMS THREDDS server. For this tutorial we are using the OPeNDAP URL for the file “EREEFS_AIMS-CSIRO_gbr4_v2_hydro_daily-monthly-2010-09.nc”.

• selectedVariable
  The name of the variable in the netCDF file.

• selectedTimeIndex
  The time slice in the netCDF file. For example, in the netCDF file “EREEFS_AIMS-CSIRO_gbr4_v2_hydro_daily-monthly-2010-09.nc” the time steps are “days”. This means if you select selectedTimeIndex=2 it refers to the second day in this file, which is 02/09/2010.

• selectedDepthIndex
  The depth slice in the netCDF file. See the following table for a mapping of index to value:

  Index (k)	Hydrodynamic 1km model	Hydrodynamic and BioGeoChemical 4km model
  1	-140.00 m	-145.00 m
  2	-120.00 m	-120.00 m
  3	-103.00 m	-103.00 m
  4	-88.00 m	-88.00 m
  5	-73.00 m	-73.00 m
  6	-60.00 m	-60.00 m
  7	-49.00 m	-49.00 m
  8	-39.50 m	-39.50 m
  9	-31.00 m	-31.00 m
  10	-24.00 m	-23.75 m
  11	-18.00 m	-17.75 m
  12	-13.00 m	-12.75 m
  13	-9.00 m	-8.80 m
  14	-5.25 m	-5.55 m
  15	-2.35 m	-3.00 m
  16	-0.50 m	-1.50 m
  17	n/a	-0.50 m

# OPeNDAP URL to file "EREEFS_AIMS-CSIRO_gbr4_v2_hydro_daily-monthly-2010-09.nc". Hydrodynamic 4km model, daily data for the month September 2010
input_file <- "http://thredds.ereefs.aims.gov.au/thredds/dodsC/s3://aims-ereefs-public-prod/derived/ncaggregate/ereefs/gbr4_v2/daily-monthly/EREEFS_AIMS-CSIRO_gbr4_v2_hydro_daily-monthly-2010-09.nc"
input_file

[1] "http://thredds.ereefs.aims.gov.au/thredds/dodsC/s3://aims-ereefs-public-prod/derived/ncaggregate/ereefs/gbr4_v2/daily-monthly/EREEFS_AIMS-CSIRO_gbr4_v2_hydro_daily-monthly-2010-09.nc"

# The "temperature" variable
selectedVariable <- "temp"
selectedVariable

[1] "temp"

# 2nd of September 2010
selectedTimeIndex <- 2
selectedTimeIndex

[1] 2

# -1.50 m depth
selectedDepthIndex <- 16
selectedDepthIndex

[1] 16

Read in the netCDF file contents

Read in the netCDF file contents and store the latitude, longitude, time and depth data in variables for later use.

nc_data <- nc_open(input_file)

lon <- ncvar_get(nc_data, "longitude")
head(lon) # look at the first few entries in the longitude vector

[1] 142.1688 142.1988 142.2288 142.2588 142.2888 142.3188

lat <- ncvar_get(nc_data, "latitude")
head(lat) # look at the first few entries in the latitude vector

[1] -28.69602 -28.66602 -28.63602 -28.60602 -28.57602 -28.54602

time <- ncvar_get(nc_data, "time")
head(time) # look at the first few entries in the time vector

[1] 7548 7549 7550 7551 7552 7553

depth <- ncvar_get(nc_data, "zc")
head(depth) # look at the first few entries in the depth vector

[1] -145 -120 -103  -88  -73  -60

Read in the data from the selected variable and select the time and depth slice.

variableData.slice <- ncvar_get(nc=nc_data,
                                varid=selectedVariable,
                                start = c(1, 1, selectedDepthIndex, selectedTimeIndex),
                                count = c(length(lon), length(lat), 1, 1))
dim(variableData.slice) 

[1] 491 723

Now that all data is stored in memory, close the netCDF file.

nc_close(nc_data)
print("file closed")

[1] "file closed"

Prepare data

Save the data in a raster. We also need to transpose and flip to orient the data correctly.

r <- raster(t(variableData.slice), xmn=min(lon), xmx=max(lon), ymn=min(lat), ymx=max(lat), crs=CRS("+init=epsg:4326"))
r <- flip(r, direction='y')
dim(r)

[1] 723 491   1

Plotting

Finally plot the data.

plot(r)