---
title: "Introduction to rasterly"
author: "Zehao Xu"
date: "`r Sys.Date()`"
output:
  html_vignette:
    toc: true
vignette: >
  %\VignetteIndexEntry{Introduction to rasterly}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
  %\VignetteDepends{grid}
  \usepackage[utf8]{inputenc}
---
  
```{r library, eval = TRUE, echo = TRUE, fig.align="center", fig.width = 6, fig.height = 4, out.width = "75%", warning=FALSE, message=FALSE}
library(rasterly)
library(data.table)
library(lubridate)
library(grid)
library(plotly)
```

`rasterly` makes it easy to rapidly generate raster images for large datasets. Although the package is inspired by the Datashader library (http://datashader.org/getting_started/index.html) available for Python, `rasterly` does not attempt to reproduce all the features of Datashader.

Rather, `rasterly` offers comparable performance to Datashader when generating rasters from source data. `rasterly` attempts to provide a flexible, convenient interface which should feel familiar to users of ggplot2 and its aesthetics-based approach to customizing plots and figures.

## Data set

The dataset used in this vignette describes Uber trips taken in New York City from April 1st to September 30th of 2014.

```{r data}
# Load data
ridesRaw_1 <- "https://raw.githubusercontent.com/plotly/datasets/master/uber-rides-data1.csv" %>%
  data.table::fread(stringsAsFactors = FALSE)
ridesRaw_2 <- "https://raw.githubusercontent.com/plotly/datasets/master/uber-rides-data2.csv" %>% 
  data.table::fread(stringsAsFactors = FALSE)
ridesRaw_3 <- "https://raw.githubusercontent.com/plotly/datasets/master/uber-rides-data3.csv"  %>% 
  data.table::fread(stringsAsFactors = FALSE)
ridesDf <- list(ridesRaw_1, ridesRaw_2, ridesRaw_3) %>% 
  data.table::rbindlist()

# Extract hour of trip taken
time <- lubridate::ymd_hms(ridesDf$`Date/Time`)
ridesDf <-  ridesDf[, 'Date/Time':=NULL][, list(Lat, 
                                                Lon,
                                                hour = lubridate::hour(time), 
                                                month = lubridate::month(time),
                                                day = lubridate::day(time))]
head(ridesDf)
```

This dataset has 4,533,327 observations, and includes the variables "latitude", "longitude", "hour", "month" and "day". 

## Basic

If we were to use `graphics::plot()`, it would take several minutes to render the image. What if we "rasterized" the image instead?

```{r basic, warning=FALSE, message=FALSE, fig.width = 4, fig.height = 3}
start_time <- Sys.time()
p <- ridesDf %>% 
  rasterly(mapping = aes(x = Lat, y = Lon)) %>% 
  rasterize_points()
p
end_time <- Sys.time()
end_time - start_time
```

## `rasterly` Structure

#### Subsetting

`rasterly()` generates a parent layer containing initial settings to generate the raster, which include `plot_height`, `plot_width` among others; child layers such as `rasterize_points()` can be piped in as well. Note that "p" above is a list of environments.

```{r list return}
# A list of environments
str(p)
```

The elements in "p" can be easily extracted or replaced by `[` and `[<-`.

```{r subsetting, warning=FALSE, message=FALSE, fig.width = 4, fig.height = 3}
p["background"]
# Replace the background in child layer `rasterly_points()`
p["background", level = 2] <- "black"
p["background"]
# color_maps in both `rasterly()` and `rasterly_points()` are replaced
## fire_map is a vector of colors (as character strings) with length 256
## see `rasterly::fire_map`
p["color_map", level = 1:2] <- fire_map
p
```

  * `level` helps to define which layer to replace; the default is `1` (the parent layer generated by `rasterly()`).
  * Available states which can be extracted or replaced are listed here:
    1. Aggregation: `data`, `mapping`, `plot_width`, `plot_height`, `range`, `x_range`, `y_range`, `xlim`, `ylim`, `aesthetics`, `reduction_func`, `glyph`, `max_size`, `group_by_data_table`, `drop_data`, `variable_check`
    2. Display: `background`, `color_map`, `color_key`, `alpha`, `span`, `show_raster`, `layout`
    
#### Build rasterly by `rasterly_build()`

To retrieve display info, use `rasterly_build()`:

```{r rasterly_build}
build <- rasterly_build(p)
str(build)
```

It contains:

  * agg: aggregation matrices, a list of numerical matrices
  * image: a raster matrix (has the same dimension with aggregation matrices)
  * lims: a list of x limits and y limits for each layer
  * x_range: the range of x over all layers
  * y_range: the range of y over all layers
  * plot_height: plot height, number of rows in aggregation matrix
  * plot_width: plot width, number of columns in aggregation matrix 
  * variable_names: names of variables
  * background: background color
  * colors: color_map (color used to map in each pixel) or color_key (used for categorical variable. In general, 'color_key' would be called when "color" is set in `aes()`)

## Display

`rasterly` does not provide any functionality to display the raster image data it generates, but instead relies on other packages.

#### `plotly` graphics

  + `add_rasterly_heatmap()`: Layers are added to Plotly objects via `add_trace(...)`; `rasterly` provides the `add_rasterly_heatmap()` function which also leverages `add_heatmap()` to generate single channel heatmap overlays for Plotly figures. Multi-channel heatmaps are not currently supported; this feature will be available in an upcoming release.
    
    ```{r add_rasterly_heatmap, fig.width = 4, fig.height = 3}
    plotly::plot_ly(ridesDf, x = ~Lat, y = ~Lon) %>%
      add_rasterly_heatmap() %>% 
      layout(
        title = "Uber drives",
        x = list(
          title = "Lat"
        ),
        y = list(
          title = "Lon"
        )
      )
    ```
  
  + `plotly.rasterly()`: `plotly.rasterly` takes a `rasterly` object and returns a "plotly" object:
    
    ```{r plotly_rasterly, warning=FALSE, message=FALSE, fig.width = 4, fig.height = 3, eval = FALSE}
    # plotly
    ply <- p %>% 
      plotly.rasterly(sizing = "contain")
    ply
    ```

## API

`rasterly` application programming interface

```{r}
r <- rasterly(data = ridesDf, 
                mapping = aes(x = Lat, y = Lon))
```

#### Mapping system 

  * Set `color`

```{r set color, fig.width = 4, fig.height = 3}
r %>% 
  rasterize_points(
     mapping = aes(color = hour),
     color_key = hourColors_map,
     background = "black"
  ) -> g
g
```

Different colors represent different hours:

```{r legend, fig.width = 4, fig.height = 3}
# rasterly doesn't currently support legends, though this feature is forthcoming
plot(1:24, y = rep(1,24), col = hourColors_map, pch = 19, cex = 3)
```

The number of aggregation matrices is equivalent to the number of categories:

```{r number of aggregation matrices}
build_g <- rasterly_build(g)
# the object has only one layer, so we index into the first element
length(build_g$agg[[1]])
# 24
```

The colors attribute in "image" within `build_g` is generated via weighted arithmetic means (default) computed from the aggregation matrices. We can choose the "cover" layout to display multiple aggregation matrices:

```{r set color cover, fig.width = 4, fig.height = 3, eval = FALSE}
r %>% 
  rasterize_points(
     mapping = aes(color = hour),
     color_key = hourColors_map,
     background = "black",
     layout = "cover"
  )
```

The resulting raster will be overlaid onto the plotting surface.

  * Set `on`
  
  `reduction_func` is implemented `on` which variable

```{r set on, fig.width = 4, fig.height = 3, eval = FALSE}
r %>% 
  rasterize_points(
    reduction_func = "mean", # take the "mean" reduction function
    mapping = aes(on = -Lat)
  )  
```

  * Set `size`
  
  To control the number of pixels allocated to an observation, we can set the `size` aesthetic; when specified, the `max_size` argument provides the upper bound of the number of pixels a single observation is allocated:
  
```{r set size, fig.width = 4, fig.height = 3, eval = FALSE}
r %>% 
  rasterize_points(
    mapping = aes(size = month),
    max_size = 4
  )  
```

Currently, only `x`, `y`, `color`, `on` and `size` can be set using `aes()`.

#### Reduction function

A reduction operator function is used when aggregating data points within each bin. One option is to reduce using the mean of the points.

  * `mean` reduction function:

```{r reduction on mean, fig.width = 4, fig.height = 3, eval = FALSE}
r %>% 
  rasterize_points(
    reduction_func = "mean", # process the data points using the mean reduction function
    background = "black",    # change background to "black" from right to left (from dark to light)
    color_map = fire_map # provide a custom color_map
  )  
```

The `mean` reduction function averages the y column (default setting) for every observation. It's also possible to average over other features using the `on` aesthetic; consult the list of available reduction functions below for additional details.

`on` 

  * `any` reduction function:
  
```{r reduction on any, fig.width = 4, fig.height = 3, eval = FALSE}
r %>% 
  rasterize_points(
    reduction_func = "any",
    color_map = c("white", "black")
  )
```

Currently supported reduction functions:

+ `sum`: If `on` is not provided within `aes()`, the default is to take the sum within each bin. When `on` is specified, the function reduces by taking the sum of all elements within the variable named in `on`.

+ `any`: When `on` is provided within `aes()`, the `any` reduction function specifies whether any elements in `on` should be mapped to each bin.

+ `mean`: If `on` is not provided in mapping `aes()`, `on` would be set as variable "y" by default. When `on` is given, the `mean` reduction function takes the mean of all elements within the variable specified by `on`.

The following functions require that `on` is first provided via `aes()`:

+ `m2`: The `m2` function computes the sum of square differences from the mean of all elements in the variable specified by `on`.

+ `var`: The `var` function computes the variance over all elements in the vector specified by `on`.

+ `sd`: The `sd` function computes the standard deviation over all elements in the vector specified by `on`.

+ `first`: The `first` function returns the first element in the vector specified by `on`.

+ `last`: The `last` function returns the last element in the vector specified by `on`.

+ `min`: The `min` function returns the minimum value in the vector specified by `on`.

+ `max`: The `min` function returns the maximum value in the vector specified by `on`.

`rasterly` is in active development; please report issues and request features via https://github.com/plotly/rasterly/issues.

Future work: provide support for `ggplot2` via `geom_rasterly()`, and `loon` via `l_rasterly()`.
