Getting started | Requirements
To run this tutorial we will need oceanmap version >= 0.13. You can find the newest version on github
Introduction
Visualizing data is a crucial step in analyzing and exploring data. During the last two decades the statistical programming language R has become a major tool for data analyses and visualization across different fields of science. However, creating figures ready for scientific publication can be a tricky and time consuming task. The oceanmap package provides some helpful functions to facilitate and optimize the visualization of geographic and oceanographic data, such as satellite and bathymetric data sets. Its major functions are
- plotmap(), ggplotmap() and gglotmaply() to plot land masks
- set.colorbar(), to define the colorbar (position, ticks, colormap)
- figure(), to facilitate resizing and saving of figures in diverse formats (jpeg, png, eps, pdf and eps)
- v(), an optimized tool for creating image or contour plots from different data formats (matrix, array, raster, netcdf, binary)
These functions were written in a way that they do not require a large amount of their numerous arguments to be specified but still return nice plots. Each of these functions will be subsequently discussed. A closely related problem represents the extraction of spatial data at defined positions/areas. Some examples on this will be discussed in an extra section.
Adding a colorbar to plotmap | The set.colorbar-function
Setting up the colorbar for an image plot in R is particularly tricky. Common plot functions, such as image.plot of the fields-package set the colorbar to a defined margin, but lack manipulation procedures. Setting up a colorbar by hand requires a lot of data input and is hence time consuming. In particular four different arguments are needed:
- the limits of the colorbar, defining its position, height and width
- the axis where to set tick labels (bottom, left, top, right)
- the colormap and its gradient (horizontal or vertical)
- the labels of the colorbar
The set.colorbar-function of the oceanmap-package comes with 3 ways to define the colorbar position:
defining its spatial extent (arguments cbx, cby)
defining argument cbpos with a single letter (“b”, “I”, “t”, "r) that indicates the position of the colorbar (bottom, left, top, right)
set the position of colorbar manually with the mouse cursor (run set.colorbar() without defining cbx, cby or cbpos).
n either case, the spatial extent of the thus defined colorbar position and extent will be returned (as cbx and cby) that can be reused in later function calls. This feature is further included in the add.region- and v-functions to set up the (default) colorbar placement of regions.
Note that the colormap and ticks of the colorbar can be defined through additional arguments of set.colorbar()
Example 1: plot colorbars manually
par(mar=c(5,8,3,8))
plot(0.5,0.5,xlim=c(0,1),ylim=c(0,1))
set.colorbar(cbx=c(0, 1), cby=c(-.3, -.4)) # bottom
set.colorbar(cby=c(0, 1), cbx=c(-.4, -.3)) # left
set.colorbar(cbx=c(0, 1), cby=c(1.2, 1.3)) # top
set.colorbar(cby=c(0, 1), cbx=c(1.2, 1.3)) # right
Example 2: use cbpos
par(mar=c(5,8,3,8))
plot(0.5,0.5,xlim=c(0,1),ylim=c(0,1))
set.colorbar(cbpos='b') # bottom
set.colorbar(cbpos='l') # left
set.colorbar(cbpos='t') # top
set.colorbar(cbpos='r') # right
Example 3: interactive placement
par(mar=c(5,8,3,8))
plot(0.5,0.5,xlim=c(0,1),ylim=c(0,1))
cb <- set.colorbar() # interactive
plot (0.5,0.5,xlim=c(0,1),ylim=c(0,1))
set.colorbar(cbx=cb$cbx, cby=cb$cby) # reuse stored colorbar
Available colormaps
oceanmap comes with a set of pre-installed colormaps. The jet.light colormap is semitransparent and useful when visualizing overlapping data sets. The year.jet colormap starts and ends with similar colors and is therefore particularly useful when analyzing seasonal patterns in data sets. The haxby colormap is the standard colormap for bathymetry data, starting with a white color.
Example 4: available colormaps
[1] "ano" "bathy" "blue"
[4] "chla" "green" "haxby"
[7] "jet" "rainbow" "red"
[10] "sst" "haxbyrev" "orange"
[13] "year.jet" "dark.blue" "jetrev"
[16] "light.jet" "white2black" "black2white"
[19] "white2darkgrey"
Instead of using preinstalled colormaps, you can also define your own one:
Example 4: own colormaps
setwd(system.file("test_files", package="oceanmap"))
gz.files <- Sys.glob('*.gz')
# figure(width=15,height=15)
par(mfrow=c(4,5))
# for(n in names(cmap)) v(gz.files[2], v_area='lion', pal=n, adaptive.vals=TRUE, main=n)
## define new color maps from blue to red to white:
n <- colorRampPalette(c('blue','red','white'))(100)
# v(gz.files[2], v_area='lion', pal=n, adaptive.vals=TRUE, main="own colormap")
The v-function: Creating image-plots of spatial data
The v-function combines the advantages of the plotmap and set.colorbar to visualize 2D (oceanographic) data. To facilitate plotting, it is further capable of handling different input formats. Valid formats are:
- matrix and array objects
- objects of class ‘Raster’ (‘RasterLayer’, ‘RasterStack’ or ‘RasterBrick’),
- netcdf-files ‘nedf4’-objects,
- ‘.gz’-files (Hervé Demarq, IRD)
Note that for plotting the input data will be internally transformed into a Raster-object, no matter which of the valid input format was selected. Examples to visualize matrix and array objects will not be discussed here but are given in ?matrix2raster(). The other input data formats are discussed in the following sections.
v() and RasterLayers:
The standard object type to handle spatial data in R are Raster-objects. Here some examples to visualize such objects with v().
Example 1: load & plot a sample Raster-object
setwd(system.file("test_files", package="oceanmap"))
load("medw4_modis_sst2_4km_1d_20020705_20020705.r2010.0.qual0.Rdata")
dat <- raster::crop(dat,extent(c(0,10,40,44))) ## crop data, xlim/ylim not yet implemented in v()
print(dat)
class : RasterLayer
dimensions : 96, 240, 23040 (nrow, ncol, ncell)
resolution : 0.04166667, 0.04166667 (x, y)
extent : 0, 10, 40, 44 (xmin, xmax, ymin, ymax)
crs : +proj=longlat +ellps=WGS84
source : memory
names : layer
values : 16.8, 25.8 (min, max)v(dat, main="Raster-object", cbpos='r')
v() and nefcdf-files
Example 2: load & plot sample netcdf-file (‘.nc’-file)
option a) load netcdf-file with ncdf4-package and plot it
setwd(system.file("test_files", package="oceanmap"))
ncfiles <- Sys.glob('*.nc') # load list of sample-'.nc'-files
par(cex=.7)
print(ncfiles)
[1] "herring_lavae.nc" "topo.v02.nc" File herring_lavae.nc (NC_FORMAT_CLASSIC):
1 variables (excluding dimension variables):
float Conc[lon,lat,time]
_FillValue: -1
3 dimensions:
lon Size:352
lat Size:310
time Size:4 *** is unlimited ***
units: seconds since 2006-03-01 00:00:00
option b) load and plot netcdf-file as RasterStack object
setwd(system.file("test_files", package="oceanmap"))
nc <- nc2raster(ncfiles[1])
File herring_lavae.nc (NC_FORMAT_CLASSIC):
1 variables (excluding dimension variables):
float Conc[lon,lat,time]
_FillValue: -1
3 dimensions:
lon Size:352
lat Size:310
time Size:4 *** is unlimited ***
units: seconds since 2006-03-01 00:00:00
converting array to RasterStack objectATTENTION: nc2raster()-calls require a “varname” selection if the netcdf-file holds multiple spatial variables.
setwd(system.file("test_files", package="oceanmap"))
## attention you may have to increase margin spacings
## for correct colorbar text rendering
par(cex=.7)
v(nc, cbpos="r") # plot RasterStack object
option c) plot netcdf-file directly
setwd(system.file("test_files", package="oceanmap"))
ncfiles <- Sys.glob('*.nc') # load list of sample-'.nc'-files
par(mfrow=c(2,1),mar=c(2,5,2,5)*1.2,cex=.6)
v(ncfiles[1], cbpos="r")
v(ncfiles[1], cbpos="r", replace.na=TRUE)
plot multiple layers:
setwd(system.file("test_files", package="oceanmap"))
par(mfrow=c(2,2),mar=c(5,5,5,5))
v(ncfiles[1], t=1:4, cbpos="r")
ggplotmap and raster files
The newest version of oceanmap includes now a gglplot and plotly-version of the base-function plotmap(). You can use gglplotmap() to create land masks, which you can add to existing ggplots. This is particularly helpful when you are dealing with rasterized data:
library(dplyr)
library(ggplot2)
data(cmap)
setwd(system.file("test_files", package="oceanmap"))
nc <- nc2raster(ncfiles[1])
File herring_lavae.nc (NC_FORMAT_CLASSIC):
1 variables (excluding dimension variables):
float Conc[lon,lat,time]
_FillValue: -1
3 dimensions:
lon Size:352
lat Size:310
time Size:4 *** is unlimited ***
units: seconds since 2006-03-01 00:00:00
converting array to RasterStack objectrs2df <- nc[[1]] %>% ## take first layer
rasterToPoints() %>% ## convert raster to xyz matrix
as.data.frame() ## convert to data frame
names(rs2df) <- c("Lon","Lat","Conc") ## reset names (important for ggplotmaply hover text)
ggobj <- ggplot() + geom_raster(data = rs2df, aes(x=Lon,y=Lat,fill=Conc))
ggobj_with_land_mask <- ggplotmap(add_to = ggobj) +
scale_fill_gradientn(colours=cmap$jet) # change colorbar
ggobj_with_land_mask
ggplotmaply | The ggplotly version for ggplotmap
Now let’s convert this figure into an interative visualization using the ggplotmaply-function: (Note that the scale bars of ggplotmap() can not be inlcuded in the plotly visualizations.)
ggplotmaply(ggobj_with_land_mask)
Although the remaining tutorial will focus on the base-functions of oceanmap, you can easily visualize the example data below with ggplotmap() by using swapping the rasterized netcdf data in ggplotnmap()-example above by the rasterized example data below (E.g. the get.bathy() from the next section creates directly a raster object.)
v() and get.bathy():
Gathering oceanographic data is not yet implemented in oceanmap. However, bathymetric data can already be obtained from the NOAA ETOPO1! data base https://sos.noaa.gov/datasets/etopo1-topography-and-bathymetry/) thanks to the get.bathy() function of oceanmap. In order to download the “bedrock” bathymetry of a region, this function requires information on the region extent that can be provided as
- longitude and latitude coordinates (lon, lat)
- an extent (raster)-object
- a region-keyword
Given this information get.bathy() downloads the bathymetric data, plots and converts it in a RasterLayer-object. The default resolution of the grid is 4 minutes, but can be changed by through the resolution statement of get.bathy(). Note that for resolutions < 1 min, the NOAA server will interpolate the bathymetry, which is very time consuming and can be done the same way, but much faster with the function disaggregate() of the raster package.
calling get.bathy() with coordinates (lon, lat):
Example 1: load & plot bathymetry of the Baltic Sea, defined by longitudes and latitudes
Example 2: load & plot bathymetry data from the NOAA-ETOPO1 database
par(mfrow=c(2,1),cex=.7)
bathy <- get.bathy("medw4", terrain=T, res=3, keep=T, visualize=T, grid=F)
# load('bathy_medw4_res.3.dat',verbose = T); bathy <- h
v(bathy, param="bathy", terrain=F, levels=c(200 ,2000)) # show contours
How to produce only contour plots from bathymetry data: set v_image=FALSE
b) only contour lines:
par(mfrow=c(2,1),cex=.6,mar=c(3,5,3,5))
h <- get.bathy("lion",terrain=FALSE, res=1, visualize=TRUE,
v_image = FALSE, levels=c(200,2000) )
## use v-function for same plot but on subregion:
v(h,v_area = "survey", param="bathy",
v_image = FALSE, levels=c(200,2000))
v() and gz-files:
Example 3: plot sample-‘.gz’-file
owd <- getwd()
setwd(system.file("test_files", package="oceanmap"))
gz.files <- Sys.glob('*.gz')
# v(gz.files[2]) ## plot content of gz-file
load and plot gz-file as a Raster-object (standard spatial object in R):
Example 4: load sample-‘.gz’-file manually as Raster-object and plot it
# obj <- readbin(gz.files[2],area='lion')
# par (mfrow=c(1,2))
# v(obj ,param="sst")
# v(obj,param="Temp") ## note unset "pal" (colormap) for unkown "param"-values!
available oceanographic parameters
The standard param values required in v() to define the title and colormap of the colorbar are stored in the parameter definitions dataset. Here some examples of pre-installed definitions.
[1] "param" "a" "b"
[4] "c" "log" "name1"
[7] "unit" "pal1" "minv"
[10] "maxv" "min" "max"
[13] "invalid_data_dc" "coast_dc" "land_dc"
[16] "no_data_dc" # ?parameter_ definitions
setwd(system.file("test_files", package="oceanmap"))
# figure (width=12, height=6.2)
par(mfrow=c(2,3))
# v('*sst2*707*' ,v_area="medw4",main="sst")
# v('*chla*531*' ,v_area="medw4" ,main="chla")
# v('*chlagrad*' ,v_area="medw4" ,main="chlagrad")
# v('?*p100*' ,v_area="medw4",main="p100 (oceanic fronts)")
# v('*xsla*',v_area="medw4",main="sla")
# h <- get.bathy("medw4",visualize=TRUE,terrain=F,res=4, main="bathy")
Extra: Extracting data from spatial Raster-objects:
Extracting spatial data from Raster-objects can be easily done with a bunch of functions implemented in the sp and maptools packages (in particular with the function extract). Here some examples, for points, polygons and circles:
printing file: lion_bathy #' extract values
# Points <- locator()
Points <- cbind(lon=c(6,8),lat=c(40,40))
points(Points,pch=19,col=c('black','red'))
extract(bathy, Points) # extract points
[1] 2809 1102#' polygon
library('rgeos')
cds1 <- rbind(c(3,42), c(3,43), c(7,43), c(3,42)) # polygon need to be closed (first point = last point)
poly <- SpatialPolygons(list(Polygons(list(Polygon(cds1)), 1)))
plot(poly[1],add=T,lwd=2,border="blue")
extract(bathy, poly[1])[[1]] %>% head()
[1] NaN NaN 10 7 26 29#
poly2 <- gBuffer(poly[1], width = 10/60) # increase polygon by 10 nm
plot(poly2,add=T,lwd=2,border="blue",lty='dotted')
# source("SpatialCircle.r")
library(maptools)
circ <- SpatialCircle(Points[1,1],Points[1,2],r = 1)
circ_poly <- PolySet2SpatialPolygons(SpatialLines2PolySet(circ)) # convert circle from spatialline to spatialpolygon
plot(circ_poly, add=T)
extract(bathy, circ_poly)[[1]] %>% head()
[1] 2668 2668 2669 2669 2673 2675