---
title: "Mapping large datasets"
output: rmarkdown::html_vignette
vignette: >
  %\VignetteIndexEntry{Mapping large datasets}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

```{r, include = FALSE}
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>",
  fig.width = 7,
  fig.height = 4.5,
  dpi = 96
)
run <- all(vapply(c("sf", "dplyr", "terra", "maps"), requireNamespace, logical(1), quietly = TRUE))
knitr::opts_chunk$set(eval = run)
```

Drawing a hexagon for every cell becomes slow long before you run out of data:
a million cells means a million polygons. globalhex avoids this by keeping
data as cell IDs and building polygons only for the cells that end up on the
screen.

```{r setup}
library(globalhex)
```

## Cell data

`hex_data()` pairs a vector of cells with one or more columns of values. It
is a tibble, and it contains no geometry:

```{r}
storms <- dplyr::storms
cells <- point_to_cell(storms[c("long", "lat")], resolution = 6)

x <- hex_data(cells, wind = storms$wind, observations = 1L)
x
```

A cell can appear many times. Here each row is one storm observation, and
several observations fall in the same cell.

## Plotting

`plot()` draws a `hex_data` object. Use `value` to choose a column (the first
one is used by default):

```{r}
plot(x, value = "wind", resolution = 4)
```

Before drawing, values are summarized from the data's resolution up to the
display resolution. The summary depends on the column type:

| Column type             | Default summary |
|-------------------------|-----------------|
| numeric                 | mean            |
| logical                 | majority        |
| character or factor     | most common value |

Counts and other totals should be added up instead, with `fun = sum`:

```{r}
plot(x, value = "observations", fun = sum, resolution = 4, border = NA)
```

## Choosing the display resolution

If you leave out `resolution`, globalhex picks the finest resolution whose
visible cells fit on the current device at a few pixels each.
`plot_resolution()` reports what it would choose:

```{r}
plot_resolution(x)
```

You can steer the automatic choice instead of fixing it:

- `max_cells` caps the number of cells drawn.
- `cell_pixels` sets the smallest cell size, in screen pixels.
- `lod = FALSE` turns off automatic aggregation and draws the data at its own
  resolution.

```{r}
plot(x, value = "wind", max_cells = 2000)
```

## Zooming in

`extent` (or the familiar `xlim` and `ylim`) limits the map to a
longitude/latitude window. Cells outside it are dropped before aggregation, so
zooming in allows finer detail. Use `reset = FALSE` to keep drawing on the same
map afterwards, here to add coastlines:

```{r}
world <- sf::st_as_sf(maps::map("world", plot = FALSE, fill = TRUE))

plot(x, value = "observations", fun = sum, resolution = 5,
     xlim = c(-98, -76), ylim = c(17, 32), reset = FALSE)
plot(sf::st_geometry(world), add = TRUE, border = "black", lwd = 0.6)
```

For a window that crosses the 180° meridian, give an `extent` whose `xmin` is
larger than its `xmax`, such as `c(170, -170, -10, 10)`.

## Map projections

Set `crs` to any projection understood by `sf::st_crs()`. Values are summarized
first and only the final cells are projected:

```{r}
temperature <- terra::rast(nrows = 180, ncols = 360)
xy <- terra::xyFromCell(temperature, seq_len(terra::ncell(temperature)))
terra::values(temperature) <- 30 * cos(xy[, 2] * pi / 180) - 10
names(temperature) <- "temperature"

global <- raster_to_cell(temperature, resolution = 4, progress = FALSE)
global

plot(global, crs = "+proj=eqearth")
```

## Large maps

Drawing thousands of separate polygons is slow on most graphics devices and
produces very large PDF or SVG files. When more than 20,000 cells would be
drawn and terra is installed, `plot()` therefore draws them as a single image
instead. The `render` argument controls this:

- `"auto"` (default) chooses based on `raster_threshold`.
- `"polygon"` always draws separate polygons, which you need for cell borders.
- `"raster"` always draws an image.

## Summarizing without plotting

`cell_aggregate()` performs the same summary as `plot()` and returns the
result, so you can analyze it or plot it another way:

```{r}
coarse <- cell_aggregate(x, resolution = 3, fun = sum, value = "observations")
coarse[order(-coarse$observations), ]
```

To use other plotting tools, such as ggplot2, convert the cells to polygons with
`sf::st_as_sf()`. Aggregate first so that you only build the polygons you need:

```{r, eval = FALSE}
polygons <- sf::st_as_sf(coarse$cell)
polygons$observations <- coarse$observations

library(ggplot2)
ggplot(polygons) +
  geom_sf(aes(fill = observations), color = NA)
```

`sf::st_as_sf()` prepares the polygons for flat longitude/latitude maps. Cells
that cross the 180° meridian are split into two parts, one at each edge of the
map, and cells containing a pole are closed along the top or bottom edge. Areas
from `sf::st_area()` are unchanged. For a Pacific-centred map, shift the
longitudes to 0–360° with `sf::st_shift_longitude()`:

```{r, eval = FALSE}
pacific <- sf::st_shift_longitude(polygons)
```

Use `sf::st_as_sf(cells, wrap = FALSE)` to get the unsplit spherical polygons,
for example to project them to a polar map yourself.