Scanfish from OpenEarth OPeNDAP

Willem Stolte

January 8, 2016

require("RNetCDF")
FALSE Loading required package: RNetCDF
library(rworldmap)
FALSE Loading required package: sp
FALSE ### Welcome to rworldmap ###
FALSE For a short introduction type :    vignette('rworldmap')
library(ggplot2)
library(scales)
library(plyr)

Download netcdf scanfish data

Scanfish data are available via the OpenEarth OPeNDAP server at Deltares. The NetCDF file is first downloaded to a local directory before it is read into memory. Direct acces to NetCDF files using R is possible but this is not yet a standard procedure under MS-Windows.

setwd("d:/stolte/temp/")

url_scanfish <-"http://opendap.deltares.nl/thredds/fileServer/opendap/rijkswaterstaat/scanfish/scanfish_temperature.nc" 

# download the netcdf file
download.file(url_scanfish, "scanfish_temperature.nc" , method = "auto",
              quiet = FALSE, mode="wb", cacheOK = TRUE)

# check directory
list.files()
FALSE [1] "1"                       "2"                      
FALSE [3] "scanfish_temperature.nc" "scanfishtemperature.csv"
# open connection to nc file
con <- open.nc("scanfish_temperature.nc")

Before plotting, the data are converted into a dataframe. the RNetCDF package function read.nc is used for this in combination with the standard function data.frame.

# make into dataframe
df <- data.frame(read.nc(con))
# turn time array into something useful
df$date <- as.Date(df$TIME, origin = "1970-01-01 00:00:00", tz = "UTC")
df$datetime <- as.POSIXct(df$date)
df$months <- format(df$date, "%m")

The data are now plotted on a static map, using a background map from the rworldmap package.

# load country polygon data
data(countriesLow)
# turn into dataframe for plotting
world <- fortify(countriesLow)
FALSE Regions defined for each Polygons
# make subset of scanfish data
df2 <- subset(df, df$date >= as.Date("2011-05-01"))
# define bbox for plotting
xxlim <- expand_range(range(df2$lon), mul = 0.1)
yylim <- expand_range(range(df2$lat), mul = 0.1)

p <- ggplot(df2, aes(lon, lat))
p + 
  geom_polygon(data = world,
               aes(x=long, y=lat, group=group),
               color = "lightgrey", fill = "darkgrey") +
  geom_jitter(aes(color = temperature), order = desc(df$z)) +
  facet_wrap(~ months) +
  coord_map(xlim = xxlim, ylim = yylim)

time plot

p <- ggplot(df2, aes(datetime, -z))
p + geom_point(aes(color = temperature))

xy plot

p <- ggplot(df2, aes(-z, temperature))
p + geom_point(aes(color = months)) + 
  coord_flip()

  • No labels