7. Trajectory objects and storm analysis

The esd package includes functions for statistical analysis and visualisation of trajectory data, e.g., the paths of extra-tropical cyclones or ice bergs. Many of the standard tools and methods are applicable to ‘trajectory’ objects, e. g., limiting the range of data by ‘subset’, visualising the trajectories by ‘map’ or ‘plot’, as well as principal component analysis and downscaling (‘PCA’ and ‘DS’).

The trajectory methods have been designed with the analysis of storm tracks in mind, in particular cyclone path data produced under the IMILAST (Inter-comparison of mid latitude storm diagnostics) project (Neu et al., 2012). Trajectory data that follow the format specified for IMILAST can be imported into R using the function ‘read.imilast’. The IMILAST data must then be transformed into a trajectory object with the function ‘trajectory’ before other esd methods can be applied.

> x <- read.imilast(filename,path=pathtofile)
> y <- trajectory(x)

The function ‘trajectory’ transforms a data frame containing spatio-temporal information about trajectories of variable length to a matrix with the trajectories interpolated to the same length (by default 10). The input to ‘trajectory’ must for every time step include a ‘Trajectory’ id number, as well as geographical (‘Lat’*itude, ‘Lon’*gitude) and temporal (‘Year’, ‘Month’, ‘Day’, and ‘Time’) information, and optionally a quality flag (‘Code99’). Other parameters describing the evolution of the trajectory may also be included and will then be interpolated and passed on to the ‘trajectory’ object. Meta data can be entered as arguments to the ‘trajectory’ function, e.g. a description of the parameter (‘param’ or ‘longname’), the source of the data (‘src’), or references to a paper, website or method (‘reference’, ‘URL’, ‘method’).

A sample file, ‘imilast.M03’, containing trajectories of deep cyclones in the North Atlantic region comes with esd and is provided in its examples. These storm paths were obtained by cyclone identification in a gridded data set (ERAinterim, 1.5◦ resolution) using a calculus based cyclone identification (CCI) algorithm (Benestad and Chen) (method 3 of IMILAST). In the future, ‘esd’ will be extended to cyclone identification and tracking using this CCI method.

The example below demonstrates how to select a subset of trajectories and plot the annual storm count using the function ‘plot.trajectory’.

# Sample trajectory object
data(imilast.M03)
# Select storm trajectories in region 40-60N
x <- subset(imilast.M03,is=list(lat=c(40,60)))
# Select winter storms
djf <- subset(x,it='djf')
# Calculate the seasonal storm count
n <- count.trajectory(djf,by='year')
# Plot the annual storm count
plot(x,new=TRUE,ylim=c(0,60),col='black')
# Add storm count for winter (djf)
lines(n,col='blue',lty=2)
legend('topright',c('all year','djf'),lty=c(1,2),col=('black','blue'))

The spatial extent of the trajectories can be plotted either as individual tracks on a map by the function ‘map.trajectory’, or as the number density (unit: year−1 1000km−2 ) by ‘map.density.trajectory’, as seen in the following example.

data(imilast.M03)
#Map storm trajectories for the winter season (djf)
map(imilast.M03,it='djf',projection='sphere')
# Map number density of storms (per year and 1000km2)
map.density.trajectory(imilast.M03,it='djf',projection='lonlat',
     xlim=c(-90,60),ylim=c(30,90),dx=4,dy=2)

Principal component analysis of ‘trajectory’ objects is by default applied to the longitude and latitude displacement with regards to the starting point of the individual trajectories.

# Sample trajectory object
data(imilast.M03)
# Perform PCA of storm tracks
pca <- PCA.trajectory(imilast.M03,anomaly=TRUE)
# Plot PCA
plot.pca.trajectory(pca)
# Show PCA on map
map.pca.trajectory(pca)

Below is a demonstration for storm statistics based on ERA5:

data(storms)
map(as.trajectory(storms),col=rgb(1,0,0,0.2),lwd=3,new=FALSE)

## Get the number of storms
N <- as.station(storms)
plot(N)

## Get the central pressure of the storms
slp <- as.station(storms, param="pcent",FUN='min')
plot(slp)
lines(trend(slp))
grid()