R script for reading netcdf files.
# map_ereefs.R by Barbara Robson, CSIRO Land and Water
#
# This script plots the surface layer of a variable from an eReefs output file, optionally
# overlain on a Google Earth map.
#
# The script should work under Linux or MacOS. The R netcdf4 library doesn’t handle opendap
# served files under Windows, unfortunately.
library(‘ncdf4’)
library(ggplot2)
library(ggmap)
source(“plume_class.R”) # needed only if you want to map flood plume optical classes
# User settings here:
# Date to map:
year <- 2011
month <- 1
day <- 30
# Variable you want to map:
var_name <- “FineSed”
#var_name <- “plume_rms”
# Location of opendap-served netcdf file:
input_stem <- “http://dapds00.nci.org.au/thredds/dodsC/fx3/gbr4_bgc_GBR4_H2p0_B2p0_Chyd_Dcrt/gbr4_bgc_simple_”
Google_map_underlay <- TRUE
# We take x_grid and y_grid from a hydrodynamic model netcdf file as they aren’t present in the latest bgc output
input_grid <- “http://dapds00.nci.org.au/thredds/dodsC/fx3/gbr4_hydro_all/gbr4_all_2016-09.nc”
#input_grid <- ‘http://cmar-02-mel.it.csiro.au:8080/opendap/scratch/rob575/gbr_924/tricho_924_2011-06.nc’
nc <- nc_open(input_grid)
x_grid <- ncvar_get(nc, “x_grid”)
y_grid <- ncvar_get(nc, “y_grid”)
nc_close(nc)
if (var_name==”plume_rms”) {
inputfile <- paste0(input_stem, format(as.Date(paste(year, month, 1, sep=’-‘)), ‘%Y-%m’),
‘.nc?R_412,R_443,R_488,R_531,R_547,R_667,R_678,time’)
nc <- nc_open(inputfile)
R_412 <- ncvar_get(nc, “R_412”, start=c(1,1,day), count=c(-1,-1,1))
R_443 <- ncvar_get(nc, “R_443”, start=c(1,1,day), count=c(-1,-1,1))
R_488 <- ncvar_get(nc, “R_488”, start=c(1,1,day), count=c(-1,-1,1))
R_531 <- ncvar_get(nc, “R_531”, start=c(1,1,day), count=c(-1,-1,1))
R_547 <- ncvar_get(nc, “R_547”, start=c(1,1,day), count=c(-1,-1,1))
R_667 <- ncvar_get(nc, “R_667”, start=c(1,1,day), count=c(-1,-1,1))
R_678 <- ncvar_get(nc, “R_678”, start=c(1,1,day), count=c(-1,-1,1))
rsr <- list(R_412, R_443, R_488, R_531, R_547, R_667, R_678)
ems_var <- plume_class(rsr)
} else {
inputfile <- paste0(input_stem, format(as.Date(paste(year, month, 1, sep=’-‘)), ‘%Y-%m’),
‘.nc?time,’, var_name)
nc <- nc_open(inputfile)
ems_var <- ncvar_get(nc, var_name, start=c(1,1,47,day), count=c(-1,-1,1,1))
# Note that for 2D variables, you need start=c(1,47,day) and count=c(-1,-1,1)
# Also note layer 47 is the surface layer, so adjust that if you want a different layer.
}
#ds <- as.Date(ncvar_get(nc, “time”), origin = as.Date(“1990-01-01”))
nc_close(nc)
a <- dim(ems_var)[1]
b <- dim(ems_var)[2]
# Set up the polygon corners. 4 per polygon.
gx <- c(x_grid[1:a, 1:b], x_grid[2:(a+1), 1:b], x_grid[2:(a+1), 2:(b+1)], x_grid[1:a, 2:(b+1)])
gy <- c(y_grid[1:a, 1:b], y_grid[2:(a+1), 1:b], y_grid[2:(a+1), 2:(b+1)], y_grid[1:a, 2:(b+1)])
gx <- array(gx, dim=c(a*b,4))
gy <- array(gy, dim=c(a*b,4))
# Find and exclude points where not all corners are defined
gx_ok <- !apply(is.na(gx),1, any)
gy_ok <- !apply(is.na(gy),1, any)
# (gx_ok and gy_ok should be identical, but let’s be certain)
gx <- c(t(gx[gx_ok&gy_ok,]))
gy <- c(t(gy[gx_ok&gy_ok,]))
# Values associated with each polygon at chosen timestep
#n <- c(ems_var[,,tstep])[gx_ok&gy_ok]
n <- c(ems_var)[gx_ok&gy_ok]
# Unique ID for each polygon
id <- 1:length(n)
id <- as.factor(id)
values <- data.frame(id = id, value = n)
positions <- data.frame(id=rep(id, each=4), x = gx, y = gy)
datapoly <- merge(values, positions, by = c(“id”))
if (Google_map_underlay) {
MapLocation<-c(min(x_grid, na.rm=TRUE)-0.5,
min(y_grid, na.rm=TRUE)-0.5,
max(x_grid, na.rm=TRUE)+0.5,
max(y_grid, na.rm=TRUE)+0.5)
myMap<-get_map(location=MapLocation, source=”google”, maptype=”hybrid”, crop=TRUE)
p <- ggmap(myMap) +
geom_polygon(aes(x=x, y=y, fill=value, group=id), data = datapoly) +
scale_fill_gradient(low=”ivory”, high=”coral4″, na.value=”transparent”, guide=”colourbar”)
} else {
p <- ggplot() +
geom_polygon(aes(x=x, y=y, fill=value, group=id), data = datapoly) +
scale_fill_gradient(low=”tan4″, high=”green2″, na.value=”transparent”, guide=”colourbar”)
# limits=c(30,36))
}
# To get the plot to display on the screen, type “p” after running this script