lec19.Rmd

---
title: "High-resolution and high-dimension data"
output: 
  html_document:
    fig_caption: no
    number_sections: yes
    toc: yes
    toc_float: false
    collapsed: no
---

```{r set-options, echo=FALSE}
options(width = 105)
knitr::opts_chunk$set(dev='png', dpi=300, cache=FALSE)
pdf.options(useDingbats = TRUE)
klippy::klippy(position = c('top', 'right'))
```
```{r load, echo=FALSE, cache=FALSE, warning=FALSE}
load(".Rdata")
```
<p><span style="color: #00cc00;">NOTE:  This page has been revised for Winter 2021, but may undergo further edits.</span></p>
# Introduction #

There has been a rapid increase in the availability of data sets that may contain thousands of observations (high-resolution, i.e. high-density or high-frequency data sets) and tens or hundreds of variables or attributes (high-dimension).  Such data sets arise from:

- satellite remote sensing and other remote-sensing methods (e.g. hyperspectral methods)
- comprehensive censuses and other polling activities
- synthesis databases of existing data sets
- "reanalysis" sets of climatological variables that may consist of multiple observations per day over spans of multiple decades
- climate-model simulations of past, present and future climates
- grid-point and other forms of data exported from GIS-managed data sets.

Although computing resources now permit the rapid production of such data sets, there is still a need for visualizing and understanding the results, and in many contexts, the analysis of the data sets is lagging the production of the data.  Several approaches are evolving::

- methods for visualizing many hundreds or thousands (or more) data points using "conventional" approaches, with about a million points now (2010’s) considered to represent a "large" data set.
- methods for examining "linked" displays (usually interactively) that allow individual observations to be tracked across multiple views of a data
- special-purpose analysis methods for dealing with high-resolution (density or frequency) data sets
- traditional and novel methods for dealing with high-dimensional data sets (i.e. those with many variables, or more variables than observations)

In `R`, some of these approaches are implemented in the `R` `mvtsplots` and `tabplots` packages, and by appropriate modification of commonplace displays like scatterplots and parallel coordinate plots included in packages available on [`CRAN`](http://cran.us.r-project.org/), plus others available in the [`Bioconductor`](http://www.bioconductor.org/) packages as extensions of traditional mulitivariate approaches motivated by large data sets produced in the various biomedical fields.

# Enhanced standard displays #

The techniques for the display and interactive analysis of high-resolution and high-dimension data sets are rapidly developing.  The approaches illustrated here make use of the `iplots` and `tabplot` packages, plus some enhanced variations of basic plots.

## Visualizing high-resolution data sets using enhanced versions of standard displays ##

The display of high resolution data (in addition to the standard approach of simply mapping it) can be illustrated using a set of precipitation-ratio data for the western United States, consisting of 3728 observations of 15 variables.  Although these data are not of extremely high resolution (or high dimension), they illustrate the general ideas.

The data are mappend in the lecture on *Maps in R*.

In what follows, we'll want to examine the large-scale patterns of the seasonality (summer-wet vs. winter-wet) of precipitation.  The data consist of monthly precipitation ratios, or the average precipitation for a particular month and station divided by the average annual total precipitation.  This has the effect of removing the very large imprint of elevation on precipitation totals.  The ratio of July to January precipitation provides a single overall description of precipitation seasonality.

## Transparency ##

A simple scatter plot (below, left) showing the relationship between January and July precipitation ratios illustrates how the crowding of points makes interpretation difficult.  The crowding can be overcome by plotting transparent symbols specified using the "alpha channel" of the color for individual points (below, right).


```{r hd05, fig.width=10}
# plot January vs. July precipitation ratios
opar <- par(mfcol=c(1,2))

# opaque symbols
plot(wus_pratio$pjanpann, wus_pratio$pjulpann, pch=16, cex=1.25, col=rgb(1,0,0))

# transparent symbols
plot(wus_pratio$pjanpann, wus_pratio$pjulpann, pch=16, cex=1.25, col=rgb(1,0,0, .2))
par <- opar
```

When there are a lot of points, sometimes the graphics capability of the GUIs are taxed, and it is more efficient to make a .pdf image directly.

```{r hd0, echo=TRUE, eval=FALSE, message=FALSE}
# transparent symbols using the pdf() device 
pdf(file="highres_enhanced_plot01.pdf")
plot(wus_pratio$pjanpann, wus_pratio$pjulpann, pch=16, cex=1.25, col=rgb(1,0,0, .2))
dev.off()
```
```{r hd07, include=FALSE, fig.width=10, echo=FALSE, eval=FALSE}
# transparent symbols using the pdf() device 
pdf(file="highres_enhanced_plot01.pdf")
plot(wus_pratio$pjanpann, wus_pratio$pjulpann, pch=16, cex=1.25, col=rgb(1,0,0, .2))
dev.off()
```

It's easy to see how the transparency of the symbols provides a visual measure of the density of points in the various regions in the space represented by the scatter plot.  A set of stripcharts illustrates the same principle.


```{r hd08, fig.width=10}
# stripcharts -- opaque symbols
opar <- par(mfcol=c(1,2))
stripchart(wus_pratio$pjanpann, xlab="PJan/Pann", method="overplot", pch=15, col=rgb(0,0,0))
stripchart(wus_pratio$pjanpann, xlab="PJan/Pann", method="stack", pch=15, col=rgb(0,0,0))
opar <- par
```

```{r hd09, fig.width=10}
# stripcharts -- alpha-channel transparency
opar <- par(mfcol=c(1,2))
stripchart(wus_pratio$pjanpann, xlab="PJan/Pann", method="overplot", pch=15, col=rgb(0,0,0,0.1))
stripchart(wus_pratio$pjanpann, xlab="PJan/Pann", method="stack", pch=15, col=rgb(0,0,0,0.1))
par <- opar
```

Over the region as a whole, the interesting question is the roles location and elevation may play in the seasonality of precipitation.  The following plots show that the dependence of precipitation seasonality on elevation is rather complicated.

```{r hd10, fig.width=10, fig.height=4.5}
# seasonal precipitation vs. elevation
opar <- par(mfcol=c(1,3))
plot(wus_pratio$elev, wus_pratio$pjanpann, pch=16, col=rgb(0,0,1, 0.1))
plot(wus_pratio$elev, wus_pratio$pjulpann, pch=16, col=rgb(0,0,1, 0.1))
plot(wus_pratio$elev, wus_pratio$pjulpjan, pch=16, col=rgb(0,0,1, 0.1))
par <- opar
```
[[Back to top]](lec19.html)

# Parallel coordinate plots #

Parallel coordiate plots in a sense present an individual axis for each variable in a dataframe along which the individual values of the variable are plotted (usually rescaled to lie between 0 and 1), with the points for a particular observation joined by straight-line segments.  Usually the points themselves are *not* plotted to avoid clutter, but all data values (or "cells" in the dataframe) explictly appear on the plot, nearly a million of them in the example here.  Parallel coordinate plots can be generated using the `GGally` package.

## Global fire data ##

The global fire data set provides a good instructional example of the development and use of parallel coordinate plots.  The input data consist of the GFEDv3.1 burned fraction (of half-degree grid cells, `gfed`), several bioclimatic variables (mean annual temperature (`mat`), mean temperature of the warmest and coldest months (`mtwa` and `mtco` respectively), growing degree days (5-deg C base, `gdd5`), annual precipitation (`map`), annual equilibrium evapotranspiration (`aet`), precipitation minus evapotranspiration (`pme`), `alpha` (P-E),lightning-strike frequency (`hmrc`), population density (`gpw`), net primary productivity (`npp`), treecover (`treecov`), and two categorical variables describing "megabiomes" and potential natural vegetation types (`megabiome` and `vegtype` respectively).

Load the necessary packages:

```{r  message=FALSE}
# load packages
library(sf)
library(ggplot2)
library(tidyr)
library(GGally)
#library(gridExtra)
```

 Read the data from an existing `.csv` file.
 
```{r, echo=TRUE, eval=TRUE, results='hide'}
# read the data -- output suppressed
csvfile <- "c:/Users/bartlein/Documents/geog495/data/csv_files/global_fire.csv"
gf <- read.csv(csvfile)
str(gf)
summary(gf)
```

The two categorical (i.e. factor) variables have levels arranged in the default alphabetical order.  The following code reorders the levels into something more ecologically and climatically sensible:

```{r, echo=TRUE, eval=TRUE, results='hide'}
# reorder megabiomes -- output suppressed
megabiome_name <- c("TropF", "WarmF", "SavDWd", "GrsShrb", "Dsrt", "TempF", "BorF", "Tund", "None", "Ice")
gf$megabiome <- factor(gf$megabiome, levels=megabiome_name)

# reorder vegtypes
vegtype_name <- c("TrEFW","TrDFW","TeBEFW","TeNEFW","TeDFW","BrEFW","BrDFW","EDMFW",
                  "Savan","GrStp","ShrbD","ShrbO","Tund","Dsrt","PDRI")
gf$vegtype <- factor(gf$vegtype, levels=vegtype_name)

# check the new ordering of factor levels
str(gf[16:17])
```

Drop the last two categories

```{r}
# drop last two categories
mb2 <- c("TropF", "WarmF", "SavDWd", "GrsShrb", "Dsrt", "TempF", "BorF", "Tund")
gf <- droplevels(gf[gf$megabiome %in% mb2, ])
table(gf$megabiome)
```

Do the usual transformations:

```{r, eval=TRUE, echo=FALSE, results='hide'}
gf$hrmc <- sqrt(gf$hrmc)
min(log10(gf$gpw[gf$gpw > 0]))
gf$gpw <- log10(gf$gpw + 1e-10)
gf$map <- sqrt(gf$map)
gf$pme <- sqrt(gf$pme - min(gf$pme))
```
```{r, eval=FALSE, echo=TRUE}
# transformations (output suppressed)
hist(gf$hrmc)
gf$hrmc <- sqrt(gf$hrmc)
hist(gf$hrmc)

hist(log10(gf$gpw[gf$gpw > 0]))
min(log10(gf$gpw[gf$gpw > 0]))
gf$gpw <- log10(gf$gpw + 1e-10)
hist(gf$gpw)

hist(gf$map)
gf$map <- sqrt(gf$map)
hist(gf$map)

hist(gf$pme)
gf$pme <- sqrt(gf$pme - min(gf$pme))
hist(gf$pme)
```

Get a shapefile too:

```{r, echo=TRUE, eval=FALSE, results='hide'}
# read a world outlines shape file
shp_path <- "/Users/bartlein/Documents/geog495/data/shp/"
shp_name <- "ne_110m_admin_0_countries.shp"
shp_file <- paste(shp_path, shp_name, sep="")
world_sf <- read_sf(shp_file)
```
```{r}
ggplot() +
  geom_sf(data = world_sf, fill="grey70", color="black") +
  scale_x_continuous(limits = c(-180, 180), breaks = seq(-180, 180, by=30)) +
  scale_y_continuous(limits = c(-90, 90), breaks = seq(-90, 90, by=30)) +
  labs(x = "Longitude", y =  "Latitude") +
  coord_sf(expand = FALSE) + 
  theme_bw()
```

## Plots ##

Parallel coordinate plots can be made using the `ggparcoord()` function from the `GGally` package:

```{r, pcp01x, eval=FALSE, echo=TRUE}
# parallel coordinate plot
ggparcoord(data = gf[1:16],
  scale = "uniminmax", alphaLines = 0.01) + 
  ylab("PCP of Global Fire Data") +
  theme_bw() +
  theme(axis.text.x  = element_text(angle=315, vjust=1.0, hjust=0.0, size=14),
        axis.title.x = element_blank(), 
        axis.text.y = element_blank(), axis.ticks.y = element_blank() )
```
```{r, pcp01, echo=FALSE}
# parallel coordinate plot
pngfile <- "pcp01.png"
png(pngfile, width=600, height=600) # open the file
ggparcoord(data = gf[1:16],
  scale = "uniminmax", alphaLines = 0.01) + 
  ylab("PCP of Global Fire Data") +
  theme_bw() +
  theme(axis.text.x  = element_text(angle=315, vjust=1.0, hjust=0.0, size=14),
        axis.title.x = element_blank(), 
        axis.text.y = element_blank(), axis.ticks.y = element_blank() )
#dev.off()
```
![](images/pcp01.png)

Note the use of transparency, specified by the `alphaLines` argument.  Individual observations (gridpoints) that have similar values for each variable trace out dense bands of lines, and several distinct bands, corresponding to different "pyromes" (like biomes) can be observed.

Parallel coordinate plots are most effective when paired with another display with points highlighted in one display similarly highlighted in the other.  In the case of the global fire data set the logical other display is a map. There are (at least) two ways of doing the highlighting.  Here, a particular megabiome ("BorF" -- Boreal forest/Taiga) is used.  First create a new variable that will be set to zero, except where `gf$megabiome == "BorF"`.  Then, the data are sorted so that the non-"BorF" observations get plotted first.

```{r}
# select those points with megabiome = "BorF"
gf$select_points <- rep(0, dim(gf)[1])
gf$select_points[gf$megabiome == "BorF"] <- 1
gf$select_points <- as.factor(gf$select_points)
table(gf$select_points)
```

```{r pcp02x, eval=FALSE, echo=TRUE}
# pcp
ggparcoord(data = gf[order(gf$select_points),],
columns = c(1:16), groupColumn = "select_points",
scale = "uniminmax", alphaLines=0.01) + ylab("")  +
theme_bw() + 
theme(axis.text.x  = element_text(angle=315, vjust=1.0, hjust=0.0, size=14),
      axis.title.x = element_blank(),
      axis.text.y = element_blank(), axis.ticks.y = element_blank(),
      legend.position = "none") +
scale_color_manual(values = c(rgb(0, 0, 0, 0.2), "red"))
```

```{r pcp02, echo=FALSE}
# pcp
pngfile <- "pcp02.png"
png(pngfile, width=600, height=600) # open the file
ggparcoord(data = gf[order(gf$select_points),],
columns = c(1:16), groupColumn = "select_points",
scale = "uniminmax", alphaLines=0.01) + ylab("")  +
theme_bw() + 
theme(axis.text.x  = element_text(angle=315, vjust=1.0, hjust=0.0, size=14),
      axis.title.x = element_blank(),
      axis.text.y = element_blank(), axis.ticks.y = element_blank(),
      legend.position = "none") +
scale_color_manual(values = c(rgb(0, 0, 0, 0.2), "red"))
#dev.off()
```
![](images/pcp02.png)

Plot the points.

```{r, eval=FALSE, echo=TRUE}  
# map
ggplot()  +
  geom_point(data = gf, aes(x = lon, y = lat, color = select_points), size = 0.8 ) +
  geom_sf(data = world_sf, aes(geometry=geometry), fill=NA, color="black") +
  theme_bw() + 
  theme(legend.position = "none") +
  coord_sf() + scale_color_manual(values = c("gray80", "red"))
```
```{r, echo=FALSE}  
# map
pngfile <- "map02.png"
png(pngfile, width=600, height=300) # open the file

ggplot()  +
  geom_point(data = gf, aes(x = lon, y = lat, color = select_points), size = 0.8 ) +
  geom_sf(data = world_sf, aes(geometry=geometry), fill=NA, color="black") +
  theme_bw() + 
  theme(legend.position = "none") +
  coord_sf() + scale_color_manual(values = c("gray80", "red"))
#dev.off()
```
![](images/map02.png)

Before moving on, delete the variable `select_points` from the `gf` dataframe:

```{r}
gf$select_points <- NULL
```

# The iplots package:  transparency and linked plots #

The `iplots` package provides another way of producing "alpha-channel" plots (plus more, as will be illustrated below).

## Installation ## 

The `iplots` package was designed to run in a "third-party" GUI called `JGR` which handles the execution of the Java-language plotting routines.  The installation of `JGR` on Windows is pretty straightforward, but a little tricky on the Mac.  

### Windows ###

On Windows, download and install the `iplots` package as usual.  (Also install the `rJava` package; other necessary packages should be automatically installed).  Then load the `rJava` and `iplots` packages as usual.

```{r hd11, message=FALSE, eval=FALSE}
library(rJava)
library(iplots)
```

### Mac OS X ###

**NOTE:  As of this writing (Feb 2021), JGR does not run correctly on macOS Catalina.  The last good instructions follow:**

On the Mac, the installation is a little more involved, `iplots` is run through an "IDE/script editor"" called `JGR`, which combines a script editor and console window, and resembles the usual Mac OS X R GUI.  `JGR` was built using an older version of Java (Java 6), and the appropriate "runtime environment" (currently Java for OS X 2015-001) has be downloaded from Apple and installed. 

Go to [https://support.apple.com/kb/DL1572?locale=en_US](https://support.apple.com/kb/DL1572?locale=en_US), and click on the "Download" button, and follow instructions.

(Note:  Apple recommends keeping Java up to date.  To download the current Java run-time environment (JRE) from Oracle.  

- Go to [[http://www.oracle.com/technetwork/java/javase/downloads/index.html]](http://www.oracle.com/technetwork/java/javase/downloads/index.html).  Click on the JRE Download button, follow the instructions, and eventually download and install `jre-8u121-macosx-x64.dmg` (or the most current version).  Then  
- Open a terminal window and type or copy and paste the following (including the quotes `"`):  
`export JAVA_HOME="/Library/Internet Plug-Ins/JavaAppletPlugin.plugin/Contents/Home"`
- Check the installation by typing `java -version` and you should see `java version "1.8.0_121-b13"`)

Then download `JGR` and the appropriate libraries:

- Go to [[http://rforge.net/JGR/files/]](http://rforge.net/JGR/files/) and download `JGR-1.8.dmg` (at the time of this writing), and then open the disk image and install it by copying it to your Applications folder. (Note that there is a `JGR` application (`JGR.app`) and the `JGR` `R` package; they’re not the same thing, and both need to be installed.)
- install `JGR-1.8.dmg` (you may have to deal with Mac security restrictions)
- Install the necessary packages the usual way (i.e. `install.packages(c("rJava", "JavaGD", "iplots"))`)

Then, to run `JGR`:

- load the `JGR` package by typing `library(JGR)`; some warnings will be produced.
- Then start `JGR` by typing `JGR()` in the `R` or `RStudio` console window. (It's ok to quit `R` or `RStudio` at this point.)
- switch to the `JGR` app.
- use the `JGR` `Workspace` menu to load a workspace, or paste the following into the bottom pane of the `JGR` app:

```{r loadWorkspace, eval=FALSE, echo=TRUE}
con <- url("https://pjbartlein.github.io/GeogDataAnalysis/data/Rdata/geog495.RData")
load(file=con)
close(con)
```

(`JRG` maintains its own workspace (and `.RData` file), distinct from the one maintained by RStudio.)

**end of last good macOS instructions**

## Scatter-plot maps ##

Here's the `iplots` version of the scatter-plot map.  Note the size of the symbols can be controlled by the up and down arrow keys, and the transparency by the right and left arrow keys.

```{r hd12, echo=TRUE, eval=FALSE}
attach(wus_pratio)
# simple scatter plot iplots version with transparency
# use arrow keys to control transparency and symbol size
iplot(elev,pjanpann, xlab="elev", ylab="pjanpann") 
detach(wus_pratio)
```
![](images/iplot011.png)

[[Back to top]](lec19.html)

## Linked interactive plots and "brushing" ##

Linked plots are those in which multiple plots are generated and viewed, and particular observations are flagged or "brushed" to call attention to them in each of the individual plots.  The iplots package provides an easy way to create a small number of linked plots.  The data set used here (precipitation ratios in the region around Yellowstone National Park) is a subset of the larger data set, and facilitates a more rapid demonstration.  The function `ipcp()` creates a parallel-coordinates plot, while the function `iplot()` in this instance creates a crude map. After the plots are created, the mouse can be used to draw a rectangular selection region that "lights up" points in both plots.  

Note:  If you get a "Data length mismatch" dialog box warning when executing the code below, reply "Yes".

### Linked plots with continuous variables ##

```{r hd15}
# linked interactive plots -- brushing
# example with Yellowstone region pratios iplots version
attach(yellpratio)
names(yellpratio)
```

```{r hd16, echo=TRUE, eval=FALSE}
# parallel coordinate plot plus scatter plot map
ipcp(yellpratio) # use arrow keys to control transparency
iplot(Lon, Lat) # use arrow keys to control symbol size
```

Here, stations with low winter precipitation ratios were selected on the parallel coordinates plot, and the corresponding locations are shown on the scatter-plot map.

![](figures_iplots/iplot021.png) ![](images/iplot022.png)

```{r hd17}
detach(yellpratio)
```

Here's a second example using the higher-density western U.S. precipitation ratios:

```{r hd18}
# example with the higher-density western U.S. pratios
attach(wus_pratio) 
```

```{r hd19, echo=TRUE, eval=FALSE}
# parallel coordiate plot plus scatter plot map
ipcp(wus_pratio) # use arrow keys to control transparency
iplot(lon, lat) # use arrow keys to control symbol size
```

In these linked plots, the stations in western Oregon and Washington were selected on the map, and the parallel coordinates plot shows their precipitation "regimes".

![](images/iplot031.png) ![](images/iplot032.png)

```{r hd20}
detach(wus_pratio)
```

### Linked plots with categorical variables ###

Data sets that contain a mixture of categorical or factor variables and ordinal or ratio-scale "continuous" variables can also be visualized, as illustrated by the Summit Cr. data set.  Here the "continuous" variables (`CumLen` and `WidthWS`) appear in boxplots and scatter plots (right), with `WidthWS` on the y-axis, and `Reach` or `CumLen` on the x-axis, while the categorical variables, `Reach` and `HU` are summarized by traditional bar plots and mosaic plots.


```{r hd22, echo=TRUE, eval=FALSE}
# linked plots -- Summit Cr. data set, including categorical variables
attach(sumcr)
imosaic(data.frame(Reach,HU))
ibar(HU)
ibar(Reach)

# add a scatter plot
iplot(CumLen, WidthWS)
iplot.opt(col=unclass(Reach)+3)

# add a boxplot
ibox(WidthWS, Reach)
```

On the plots below, pools (`P`) in the exclosure reach (`B`) were highlighted in the mosaic plot, and appear in red on the other plots.

![](images/iplot041.png) ![](images/iplot042.png) ![](images/iplot043.png)  
![](images/iplot044.png) ![](images/iplot045.png) 

```{r hd23, eval=FALSE}
detach(sumcr)
```

[[Back to top]](lec19.html)
<hr>
# High-density plots with the tabplot package #

The "tableplot" is an interesting visualization method that can be useful for discovering large-scale patterns in the covariations among a small (< 15) set or subset of variables, which can be either continuous or categorical.  The basic idea is to organize the rows of a data set along an axis of increasing or decreasing values of a continuous variable or the different "levels" of a categorical variable.  A reduction in the density of the data is achieved by binning multiple rows of the data and plotting bin averages (or frequencies) of the variables.  As in all compositing techniques, if there a consistent response in the other variables to the variable selected to organize the rows, then this will emerge, otherwise inconsistent reponses will wash each other out, and no apparent pattern will be evident in the plots of the other variables.

## North American modern pollen data ##

This example makes use of a data set of 4360 modern pollen samples for North America, from Whitmore et al. (2005, *Quaternary Science Reviews* [https://doi.org/10.1016/j.quascirev.2005.03.005](https://doi.org/10.1016/j.quascirev.2005.03.005).  The variables include the location and elevation of each sample, the relative abundances of the 12 most-abundant pollen types (transformed by taking square roots), and two categorical variables, Federova et al. (1994) vegetation types and USGS land-cover types from satellite remote-sensing data (see Whitmore et al. 2005 for details).  This data set provides an interesting example, inasmuch as the resulting tableplots are reminiscent of "downcore"-plotted pollen diagrams.

```{r, read napollen}
# read the North American modern pollen data
csv_path <- "c:/Users/bartlein/Documents/geog495/data/csv_files/"
csv_name <- "NAmodpol.csv"
csv_file <- paste(csv_path, csv_name, sep="")
napol <- read.csv(csv_file) # takes a while
str(napol)
```

Get a quick map of the pollen site locations:
  
```{r, hd31, fig.width=6, fig.height=5}
plot(st_geometry(world_sf), xlim=c(-180, -10), ylim=c(20, 80), axes=TRUE)
points(napol$Lat ~ napol$Lon, pch=16, cex=0.5, col="darkgreen")
```

## Plots ##
Load the `tabplot` package.  Note that `tabplot` loads some other packages (`ffbase`, `ff`, `bit`, etc.) that allow for efficient processing of large data sets and replace some of the functions in the `base` package.  It might be good to restart R after using the `tabplot` package.

```{r, hd32, message=FALSE}
library(tabplot)
```

Here's a simple tableplot, with the observations arranged by latitude (`sortCol="Lat"`) and linear scales specified for all continuous variables (`scales="lin"`, the `tabplot()` function attempts to find a suitable scale for each variable, but here the pollen variables were all transformed prior to the analysis, and so we want to suppress further transformation here).

(Note that it just happens that the resulting plot looks like a standard pollen diagram--this isn't designed in any way.)

```{r, hd33, echo=TRUE, fig.width=10, eval=FALSE}
tableplot(napol, sortCol="Lat", scales="lin", fontsize = 8)
```
```{r, hd34, out.height="90%", out.width="90%", message=FALSE, echo=FALSE}
suppressWarnings(tableplot(napol, sortCol="Lat", scales="lin", fontsize = 8))
```

It's possible to see broad trends in all of the pollen types, as well as in the frequency of different vegetationtypes as a function of latitude.

Each bin the above plot represents a bin that, in this case, holds 46 observations, and what is plotted as bars in each column  is the bin-average or frequency for those observations. Here's a second plot like the first, using more bins:
  
  ```{r, hd35, echo=TRUE, eval=FALSE}
tableplot(napol, sortCol="Lat", scales="lin", nBins=200, fontsize = 8)
```
```{r, hd36, out.height="90%", out.width="90%", message=FALSE, echo=FALSE}
suppressWarnings(tableplot(napol, sortCol="Lat", scales="lin", nBins=200, fontsize = 8))
```

The same broad-scale patterns are evident in this plot, suggesting that the approach is robust with respect to the number of bins.  The next plot using the abundance of *Picea* (spruce to organize the data):
  
  ```{r, hd37, echo=TRUE, eval=FALSE}
tableplot(napol, sortCol="Picea", scales="lin", nBins=100, , fontsize = 8)
```
```{r, hd38, out.height="90%", out.width="90%", message=FALSE, echo=FALSE}
suppressWarnings(tableplot(napol, sortCol="Picea", scales="lin", nBins=100, fontsize = 8))
```

The plot shows that *Aremisia* (sagebrush) pollen has a similar profile to *Picea* and both tend to be high at high elevation and low (western) longitudes.  The next plot using one of the categorical variables to organize the data.

```{r, hd39, echo=TRUE, eval=FALSE}
tableplot(napol, sortCol="Fed", scales="lin", nBins=100, fontsize = 8)
```
```{r, hd40, out.height="90%", out.width="90%", message=FALSE, echo=FALSE}
suppressWarnings(tableplot(napol, sortCol="Fed", scales="lin", nBins=100, fontsize = 8))
```

The blockiness of the different profiles shows that the individual vegetation types have distinctive mixtures of the different pollen types.

## Global fire data

The global fire data set has roughly the same number of variables as the pollen data, but ten times the number of observations.  The number of individual observations in each bin is therefore larger, increasing the likelihood that the contents of individual bins will be more heterogeneous, with an consequent loss of overall information.

## Plots ##

To start, produce a tableplot, arranged by mean-annual temperature (`mat`), an obvious first-order predictor of biomass burning, as well as overall climate (through the dependence of the water balance on temperature via the Clausius-Clapyron relationship):

```{r, out.height="90%", out.width="90%"}
attach(gf)
tableplot(gf, sortCol="mat", scales="lin", fontsize = 8, pals=list("Set8", "Set8"))
```

Overall, the trends of individual variables are relatively smooth, and it's easy to see distiguish the temperature and moisture variables (as well as the population variable `gpw`).  The vegetation-type variables `megabiome` and `vegtype` are also well sorted out along a mean annual temperature gradient.  Burned fraction (`gfed`) obviously displays a threshold-type relationshion, where values are small at low temperatures, and higher at high values of `mat`.

Next, the tableplot is sorted by precipitation minus evapotranspiration (`pme`, or P-E):  

```{r, out.height="90%", out.width="90%"}
tableplot(gf, sortCol="pme", scales="lin", fontsize = 8)
```

The relationship between burned fraction and P-E evident in univariate plots--the arch-shaped pattern, where burned fraction is highest at intermediate values of effective moisture--is also evident here:  as P-E monotonically incrases, burned fraction (`gfed`) rises, then falls.

Finally, the tableplot can be sorted by burned fraction (`gfed`):

```{r, out.height="90%", out.width="90%"}
tableplot(gf, sortCol="gfed", scales="lin", fontsize = 8)
```

This plot confirms the existing relationships, where burned fraction broadly tracks temperature and P-E.  However the plots also reveal a distinct threshold that is apparently related to a distinct discontinuity in vegetation:  the transition beteen desert and tundra vegetation types, and grassland, shubland and forested ones.

[[Back to top]](lec19.html)


# Plotting multivariate/multiple time series #

Although the tableplots of the North American pollen data above happened to look like conventional pollen diagrams, which were a comparatively early (over a century ago) attempt to visualize multivariate data, the data displayed by the `tableplot()` function does not necessarily need to be a series.  The next example illustrates a display method in which the data *do* form a series, in this case pollen data again, only this time for a single site, which form a multivariate time series.

## Carp L. pollen data ##

The data here come from a long pollen record from Carp L., in Washinton State.  (Whitlock, C.L. et al., 2000, *Palaeogeography Palaeoclimatology Palaeoecology* 155(1-2):7-29 [https://doi.org/10.1016/S0031-0182(99)00092-9] (https://doi.org/10.1016/S0031-0182\(99\)00092-9) and Whitlock, C. and P.J. Bartlein, 1997, *Nature* 388:57-61 [https://doi.org/10.1038/40380](https://doi.org/10.1038/40380).  Pollen data are "long-tailed" and so the data were transformed earlier by taking the square root of the pollen percentages.

```{r, read carp}
# read Carp L. pollen data (transformed to square-roots)
csvpath <- "c:/Users/bartlein/Documents/geog495/data/csv_files/"
csvname <- "carp96t.csv" # square-root transformed data
carp <- read.csv(paste(csvpath, csvname, sep=""))
summary(carp)
```

Pollen data (expressed as proportions or percentages of individual pollen types) vary over an order of magnitude in value, and also have long-tailed distributions.  Consequently they should be tranformed, and a useful transformation is to simply take the square root of the proportions.

```{r}
# transform 
tcarp <- carp
tcarp[4:25] <- sqrt(tcarp[4:25])
names(tcarp)
```

As it appears in the spreadsheet, the data are not particularly well organized, but instead the columns are arranged according to the traditional style of pollen diagrams.  An alternative approach is to arrange the indiviual variables (pollen types) according to their affinities as they vary over the core, which can be accomplished by doing a cluster analysis on a dissimilary matrix among pollen types: taxa that vary in a similar (or inverse) fashion downcore appear adjacent to one another.

```{r}
# cluster rearrangement
carp_heat <- heatmap(as.matrix(tcarp[4:25]), Rowv = NA, scale="none")
attributes(carp_heat)
carp_heat$colInd
sort_cols <- carp_heat$colInd + 3
```

The "heatmap" (which in itself provides an interesting way to look at the data), can be used to rearrange the columns of the dataframe to put similar (or inversely similar, when one type is high the other is low) next to one another:

```{r}
# rearrange matrix
stcarp <- tcarp
names(stcarp)
stcarp[4:25] <- tcarp[sort_cols]
names(stcarp) <- c("Depth", "Age", "IS", names(tcarp)[sort_cols])
names(stcarp)
```

## Plots ##

The idea behind the `mvtsplots()` package (and function) is to use color to reveal underlying patterns in the data.  The first plot shows the "raw" data, simply coded for each pollen taxon independently from low (magenta) to high (green).  The plots on the margins are of boxplots for the individual series (right), and a time series of the median values across all series (bottom), which attempt to summarize the overall level of the series.

```{r, mvtsplots01}
library(mvtsplot)
mvtsplot(stcarp[, 4:25], xtime=carp$Age, levels=9, margin=TRUE)
```

PLotted this way, one begins to get an idea of some organization to the data along the time axis.  The plot is a little noisy, and can be simplified by smoothing the individual series.

```{r, mtvsplot02}
mvtsplot(stcarp[, 4:25], xtime=carp$Age, smooth.df=25, levels=9, margin=TRUE)
```

Further simplification can be gotten by "normalizing" the data, transforming and assigning colors based on the range of values in each series individually.

```{r, mtvsplot03}
mvtsplot(stcarp[, 4:25], xtime=carp$Age, smooth.df=2, norm="internal", levels=9, margin=TRUE)
```

The color field now shows that the data are well organized in time, which can be verified by adding a plot of the variable `IS` i.e. the MIS (Marine Isotopic Stage) to the plot (on the last row).  The progression of plots can also be seen to have proceeded from a fairly chaotic pattern to one with coherent patches:

```{r, mvtsplot04}
mvtsplot(stcarp[, 3:25], xtime=carp$Age, smooth.df=2, norm="internal", levels=9, margin=TRUE)
```

[[Back to top]](lec19.html)

# Interactive plots

R has a long history of attempts to develop interactive plots, usually consisting of the basic plots developed from low-dimensional data. RStudio's *Shiny* package [https://shiny.rstudio.com](https://shiny.rstudio.com) and the *Plotly R Open Source Graphics Library* [https://plot.ly/r/](https://plot.ly/r/) are examples that are rapidly developing.  An interesting plot for displaying high-dimensional data is the *cubeplot*, which is easiest to explain by looking at one.  The data here come from the *TraCE-21k* transient climate simulation from 21,000 years ago to present (at monthly resolution) (Liu, Z., et al. 2009, Science 325(5938): 310-314).  A data set of of decadal-average "anomaly" values (the difference between past and present) is examined here.

## Data ##

The data are read and plotted using the `raster` and `rasterVis` framework.  First, load the packages:

```{r cube, echo=FALSE, eval=TRUE, include=FALSE, results="hide", messages=FALSE, cache=FALSE}
# load packages
library(ncdf4)
library(raster)
library(rasterVis)
library(maptools)
library(latticeExtra)
library(RColorBrewer)
library(cubeview)
```

Read the data:

```{r read tr21, echo=TRUE, eval=TRUE, results="hide", warning=FALSE, cache=FALSE}
# read TraCE21-ka transient climate-model simulation decadal data
tr21dec_path <- "c:/Users/bartlein/Documents/geog495/data/nc_files/"
tr21dec_name <- "Trace21_TREFHT_anm2.nc"
tr21dec_file <- paste(tr21dec_path, tr21dec_name, sep="")
tas_anm_ann <- brick(tr21dec_file, varname="tas_anm_ann", band=seq(1:2204))
tas_anm_ann <- rotate(tas_anm_ann)
tas_anm_ann
```

Even though the data is decadal resolution, it's still a lot of data, and so for this example, the data are (actually) decimated to produce century time-step decadal average values.

```{r tr21 subset, echo=TRUE, cache=FALSE}
tas2 <- subset(tas_anm_ann, subset=seq(1,2201, by=10))
```

A levelplot of one slice of the the data can be used to become oriented to the geographic layout of the data:

```{r plot setup, cache=FALSE}
mapTheme <- rasterTheme(region=(rev(brewer.pal(10,"RdBu"))))
cutpts <- c(-40,-10,-5,-2,-1,0,1,2,5,10,20)
col <- rev(brewer.pal(10,"RdBu"))
plt <- levelplot(subset(tas2,1), at=cutpts, margin=FALSE, par.settings=mapTheme)
plt + latticeExtra::layer(sp.lines(as_Spatial(world_sf), col="black", lwd=1))
```

## Cube plot ## 

Then, the cubeplot can be produced using the `cubeview()` function from the `cubeview` package.  The plot is interactively controlled as follows (from the `cubeview()` help file:

The location of the slices can be controlled by keys:

- x-axis: LEFT / RIGHT arrow key
- y-axis: DOWN / UP arrow key
- z-axis: PAGE_DOWN / PAGE_UP key

Other controls:

- Press and hold left mouse-button to rotate the cube
- Press and hold right mouse-button to move the cube
- Spin mouse-wheel or press and hold middle mouse-button and move mouse down/up to zoom the cube
- Press space bar to show/hide slice position guides

The cubeplot, as rendered here, displays a map on the "front (or rear) face" of the brick, and Hovmöller plots on the top (or bottom) and sides.

```{r tr21 cube, echo=TRUE, eval=FALSE, cache=FALSE}
cubeView(tas2, at=cutpts, col.regions=col)
```

![](images/tr21cube.png)

Here is a link to the plot on the server:  [https://pjbartlein.github.io/REarthSysSci/tr21cube.html](https://pjbartlein.github.io/REarthSysSci/tr21cube.html)  

[[Back to top]](lec19.html)

# Readings #

- Section 6.4 in James, G., et al. (2013) *An Introduction to Statistical Learning*, Springer. (Search by title for the eBook version on the UO Library page [[http://library.uoregon.edu]](http://library.uoregon.edu).)
- Chapter 18 in Hastie, T. et al., (2009) *The Elements of Statistical Learning*, Springer. (Search by title of the eBook version on the UO Library page [[http://library.uoregon.edu]](http://library.uoregon.edu).)  This is a technical discussion of the problems that arise in analyzing high-dimension data.  (Note, the link seems to be currently slightly broken--you may have to search again for the title once on the Springer Link web page.  Look for the 2009 edition.)