---
title: "What a grid is"
author: "Michael D. Sumner"
date: "`r Sys.Date()`"
output:
  rmarkdown::html_vignette:
    fig_width: 7
    fig_height: 5
vignette: >
  %\VignetteIndexEntry{What a grid is}
  %\VignetteEncoding{UTF-8}
  %\VignetteEngine{knitr::rmarkdown}
---

```{r setup, include = FALSE}
knitr::opts_chunk$set(collapse = TRUE, comment = "#>")
```

Every interpolation in this package produces a grid, and a grid here is a list
of four things. There is no fifth thing and no class hierarchy to learn.

```{r}
library(guerrilla)
g <- grid_spec(dimension = c(8, 6), extent = c(0, 4, 0, 3))
str(unclass(g))
```

`dimension` is how many cells across and down, in that order. `extent` is the
outer boundary as `xmin, xmax, ymin, ymax`. `crs` is what the coordinates mean,
or `NULL` when nobody said. `values` is one number per cell, filled in later by
whatever produced it.

That is genuinely all a raster is. Every raster class in R holds these same few
numbers plus a vector, and the rest is convenience. Keeping them in a plain list
means the interpolation code can be read without knowing anything about raster
classes, which is the point of the package.

# Where the cells are

Cell size is not stored, because it is `extent` divided by `dimension` and
storing it too would let the two disagree:

```{r}
grid_res(g)
grid_ncell(g)
```

Cell centres come from `grid_xy()`. That one function is everything the
interpolation code needs to know about grids:

```{r}
head(grid_xy(g))
```

Cells run left to right along the top row first, then down. That is the raster
convention everywhere, and it is why every conversion in this package and every
comparison with another package involves a `rev()` somewhere.

```{r}
plot(g_filled <- local({ g$values <- seq_len(grid_ncell(g)); g }))
```

The first cell is the top left one and holds 1.

All of the arithmetic comes from **vaster**, which does grid index work in
base R with no dependencies at all:

```{r}
vaster::cell_from_xy(g$dimension, g$extent, cbind(0.5, 2.8))
vaster::xy_from_cell(g$dimension, g$extent, 1)
```

# Making one

`grid_spec()` takes points and covers them:

```{r}
set.seed(1)
xy <- cbind(runif(50), runif(50))
grid_spec(xy)
```

The extent is exactly the range of the points, so the outermost samples sit on
the boundary. Often you want a margin, and `pad` is a proportion of the span:

```{r}
grid_spec(xy, pad = 0.1)
```

`dimension` and `extent` can be given directly, and either can come from
somewhere else entirely. `as_grid()` reads a `RasterLayer` or a `SpatRaster`,
so a target grid from another package works as-is:

```{r}
if (requireNamespace("terra", quietly = TRUE)) {
  as_grid(terra::rast(nrows = 10, ncols = 20, xmin = 0, xmax = 2))
}
```

# Handing it to something else

The converters exist because other packages do things this one does not, and
each is a handful of lines because there is nothing to translate:

```{r}
g$values <- runif(grid_ncell(g))

if (requireNamespace("raster", quietly = TRUE)) as_raster(g)
```

```{r}
if (requireNamespace("terra", quietly = TRUE)) as_terra(g)
```

`as_gdalraster()` has to write a file, because that is how GDAL works. It hands
back the open dataset, which you close when you are done with it:

```{r}
if (requireNamespace("gdalraster", quietly = TRUE)) {
  ds <- as_gdalraster(g, tempfile(fileext = ".tif"))
  print(ds$dim())
  print(ds$getGeoTransform())
  ds$close()
}
```

Note the geotransform: GDAL puts the origin at the top left with a negative y
step, which is the same top-down row order the grid uses, spelled differently.

`as.matrix()` and `as.data.frame()` are there for the cases where you just want
the numbers:

```{r}
dim(as.matrix(g))
head(as.data.frame(g), 3)
```

# What is not here

A grid has no time dimension, no bands, no attributes, no units and no
resampling. Those are all real and useful things and they are what **terra**,
**stars** and **gdalraster** are for. This package stops at the point where
you would need one of them, and hands you the numbers.