Example of downloading and plotting ACCESS met model data in R

access.r

# Demonstration - Downloading and animating ACCESS forecast data
# grant.williamson@utas.edu.au

# General functions, including str_replace, map
library(tidyverse) 

# threddscrawler is a useful package for parsing THREDDS data servers
# Install it with:
# devtools::install_github("BigelowLab/threddscrawler")
library(threddscrawler)

# Dealing with raster data
library(raster)

# Dealing with NetCDF files
library(ncdf4)

# Thematic mapping
library(tmap)

# Spatial data
library(sf)

# Load the simple world map that comes with the tmap package
data("World")

# Define the ACCESS domain you wish to plot
# New models appear every six hours
# vt = Victoria/Tasmania
# sy = Sydney
# ph = Perth
# ad = Adelaide
# bn = Brisbane
# r and g are regional/global models, but they have different
# variable names 
access_domain = "vt"

# Root URL for the ACCESS Data on BOM Opendap server
root = paste0("http://opendap.bom.gov.au:8080/thredds/catalog/bmrc/access-",access_domain,"-fc/ops/surface/latest/catalog.xml")

# Get catalog of current latest model run output
Top = get_catalog(root)
hours = Top$get_datasets()

# Number of hours avaiable varies
nhours = length(hours)

# Function to calculate the timestamp of each model file
tfn = function(x){
  # Different domains have different numbers of characters ("vt" vs "g)
  # And this changes the filename length
  shft = nchar(access_domain)
  # Name of the file
  this_name = x$name
  # Extract the date the model was run, the hour the model was run
  # And the hourly offset for each file after the run time
  model_date = substr(this_name,9+shft,16+shft)
  model_hour = substr(this_name,7+shft,18+shft)
  offset_hour = as.numeric(substr(this_name,20+shft,22+shft))
  # Cram them all together into a date format
  datestr = paste0(substr(model_date,1,4),"-",substr(model_date,5,6),"-",substr(model_date,7,8)," ",model_hour,":00:00")
  # Original file dates are in UTC, but after conversion to POSIX
  # They will appear in local time
  datepct = as.POSIXct(datestr,tz="UTC")
  # Add the hour offset for each file
  dateoffset = datepct + offset_hour * 60 * 60
  dateoffset
}

# Apply that function to calculate the timestamps over our list of filenames
times = map(hours,tfn) %>% reduce(c)

# Download files to the temp directory
# Will take a few minutes, the server seems pretty slow!
# I should parallelize this to see if it would speed it up
for(i in 1:nhours){
  # Files are stored in a different hierarchy than the THREDDS catalog
  gname=str_replace(hours[[i]]$url,"catalog","fileServer")
  download.file(gname,paste0(tempdir(),"/",hours[[i]]$name))
}  

# Get a list of all the filenames we have so we can
# stack them and clean them up
nc_list = list.files(tempdir(),pattern=".nc",full.names = TRUE)

# Make a cube (raster stack) of the rain data
# Stack the "accumulated precipitation" variable
accum_prcp = stack(nc_list,varname="accum_prcp")

# Give the raster layers names based on their timestamp
names(accum_prcp) = as.character(times)

# Save accumulated precipitation to a GeoTiff
writeRaster(accum_prcp,"accum_prcp.tif",overwrite=TRUE)

# As its name suggests, accum_prcp contains rainfall accumulating over time
# To calculate the actual rainfall amount each hour, we can calculate
# the differential
prcp = calc(accum_prcp, fun = diff)

# Give each layer in the raster a name based on the time
# This is useful for animation plotting
# Note that because we took a differential, we have one less
# layer in this stack, so we need one less timestamp
# So we use the timestamp at the end of the hour of accumulation
names(prcp) = as.character(times[2:length(times)])

# Save that rain stack to a GeoTiff file 
writeRaster(prcp,"prcp.tif",overwrite=TRUE)

# Use tmap to make a map of a single timestep
p = tm_shape(prcp[[1]],is.master = TRUE) + 
  tm_raster(style="fixed",breaks=c(0,1,2.5,5,10,25,50),palette=colorRampPalette(c("white","blue","red","black"))(7)) +
  tm_shape(World) + 
  tm_borders() 

# Plot the map
p

# Make an animation - to do this we define the facet, but make it only
# a single row and column - this will push the facets into the time
# dimension.  Note tm_facets seems to automatically use the raster layers
p = tm_shape(prcp,is.master = TRUE) + 
  tm_raster(style="fixed",breaks=c(0,1,2.5,5,10,25,50),palette=colorRampPalette(c("white","blue","red","black"))(7)) +
  tm_shape(World) + 
  tm_borders() + 
  tm_facets(nrow=1,ncol=1,free.scales.raster = FALSE)

# Use tmap_animation to save the animated map to a GIF
# Note that this requires the ImageMagick software to be installed
# - not the R "magick" package, this seems to use the actual
# command line program.  But I suspect animation using the magick package
# itself would be pretty simple
tmap_animation(p, filename="rain.gif", width=800, delay=20)

# Run this to delete the input files
# unlink(hours)

# Lots of other variables are available, eg.
# soil_mois - soil moisture content in four layers
# temp_scrn - near-surface temperature
# rh2m - relative humidity
# u10, v10 - zonal and meridional wind