# Using Quadtrees

#### Derek Friend

#### Last compiled on August 28, 2023

## Vignette content

This vignette covers ways of interacting with a `Quadtree` object after it has been created - for example, extracting values and modifying values.

    library(terra)
    library(quadtree)

    habitat <- terra::rast(system.file("extdata", "habitat.tif", package="quadtree"))
    rast <- habitat

    qt <- quadtree(rast, .15)

## Retrieving basic info about a quadtree

There are a number of functions for retrieving basic info about a `Quadtree`.

### Getting a summary of a quadtree

A basic summary of a quadtree can be shown by just typing the variable name or using the `summary()` function:

    qt
    #> class         : Quadtree
    #> # of cells    : 4189
    #> min cell size : 250
    #> extent        : 0, 64000, 0, 64000 (xmin, xmax, ymin, ymax)
    #> crs           : PROJCRS["unknown",
    #>     BASEGEOGCRS["unknown",
    #>         DATUM["Unknown based on Airy 1830 ellipsoid",
    #>             ELLIPSOID["Airy 1830",6377563.396,299.3249646,
    #>                 LENGTHUNIT["metre",1,
    #>                     ID["EPSG",9001]]]],
    #>         PRIMEM["Greenwich",0,
    #>             ANGLEUNIT["degree",0.0174532925199433,
    #>                 ID["EPSG",9122]]]],
    #>     CONVERSION["Transverse Mercator",
    #>         METHOD["Transverse Mercator",
    #>             ID["EPSG",9807]],
    #>         PARAMETER["Latitude of natural origin",49,
    #>             ANGLEUNIT["degree",0.0174532925199433],
    #>             ID["EPSG",8801]],
    #>         PARAMETER["Longitude of natural origin",-2,
    #>             ANGLEUNIT["degree",0.0174532925199433],
    #>             ID["EPSG",8802]],
    #>         PARAMETER["Scale factor at natural origin",0.9996012717,
    #>             SCALEUNIT["unity",1],
    #>             ID["EPSG",8805]],
    #>         PARAMETER["False easting",400000,
    #>             LENGTHUNIT["metre",1],
    #>             ID["EPSG",8806]],
    #>         PARAMETER["False northing",-100000,
    #>             LENGTHUNIT["metre",1],
    #>             ID["EPSG",8807]]],
    #>     CS[Cartesian,2],
    #>         AXIS["easting",east,
    #>             ORDER[1],
    #>             LENGTHUNIT["metre",1,
    #>                 ID["EPSG",9001]]],
    #>         AXIS["northing",north,
    #>             ORDER[2],
    #>             LENGTHUNIT["metre",1,
    #>                 ID["EPSG",9001]]]]
    #> values        : 0.109, 0.974 (min, max)

### Getting the projection of a quadtree

`projection()` can be used to get and set the projection of the `Quadtree`. In this case there is no projection, so an empty string is returned.

    projection(qt)
    #> [1] "PROJCRS[\"unknown\",\n    BASEGEOGCRS[\"unknown\",\n        DATUM[\"Unknown based on Airy 1830 ellipsoid\",\n            ELLIPSOID[\"Airy 1830\",6377563.396,299.3249646,\n                LENGTHUNIT[\"metre\",1,\n                    ID[\"EPSG\",9001]]]],\n        PRIMEM[\"Greenwich\",0,\n            ANGLEUNIT[\"degree\",0.0174532925199433,\n                ID[\"EPSG\",9122]]]],\n    CONVERSION[\"Transverse Mercator\",\n        METHOD[\"Transverse Mercator\",\n            ID[\"EPSG\",9807]],\n        PARAMETER[\"Latitude of natural origin\",49,\n            ANGLEUNIT[\"degree\",0.0174532925199433],\n            ID[\"EPSG\",8801]],\n        PARAMETER[\"Longitude of natural origin\",-2,\n            ANGLEUNIT[\"degree\",0.0174532925199433],\n            ID[\"EPSG\",8802]],\n        PARAMETER[\"Scale factor at natural origin\",0.9996012717,\n            SCALEUNIT[\"unity\",1],\n            ID[\"EPSG\",8805]],\n        PARAMETER[\"False easting\",400000,\n            LENGTHUNIT[\"metre\",1],\n            ID[\"EPSG\",8806]],\n        PARAMETER[\"False northing\",-100000,\n            LENGTHUNIT[\"metre\",1],\n            ID[\"EPSG\",8807]]],\n    CS[Cartesian,2],\n        AXIS[\"easting\",east,\n            ORDER[1],\n            LENGTHUNIT[\"metre\",1,\n                ID[\"EPSG\",9001]]],\n        AXIS[\"northing\",north,\n            ORDER[2],\n            LENGTHUNIT[\"metre\",1,\n                ID[\"EPSG\",9001]]]]"

### Getting the number of cells in a quadtree

`n_cells()` returns the number of cells in the quadtree. It has one optional parameter, `terminal_only` - if `TRUE`, the number of cells of the quadtree that are terminal (i.e. have no children) is returned. If `FALSE`, the total number of nodes in the quadtree is returned.

    n_cells(qt)
    #> [1] 4189
    n_cells(qt, terminal_only = FALSE)
    #> [1] 5585

### Getting the extent of a quadtree

`extent()` can be used to return the extent of the quadtree. It has one optional parameter `original` - if `FALSE` (the default), it returns the total extent covered by the quadtree. If `TRUE`, the function returns the extent of the original raster used to create the quadtree, before `NA` rows/columns were added to pad the dimensions. You may need to preface `extent()` with `quadtree::` to avoid conflicts with the `raster` package.

    quadtree::extent(qt)
    #> SpatExtent : 0, 64000, 0, 64000 (xmin, xmax, ymin, ymax)
    quadtree::extent(qt, original = TRUE)
    #> SpatExtent : 0, 40250, 0, 44500 (xmin, xmax, ymin, ymax)

## Retrieving cell-level data

### Extracting values

Values can be ‘extracted’ at point locations using the `extract()` function. This function has one optional parameter - `extents`. If `extents` is `FALSE` (the default), then only the values at the point locations are returned. If `extents` is `TRUE`, then a matrix is returned that also returns the x and y limits of each cell in addition to the cell value.

    pts <- cbind(x = c(5609, 3959, 20161, 27662, 32763),
                 y = c(10835, 29586, 31836, 10834, 36337))

    plot(qt, crop = TRUE, border_lwd = .3, na_col = NULL)
    points(pts, pch = 16)

    quadtree::extract(qt, pts)
    #> [1]       NaN 0.1451094 0.8202500 0.9556426 0.6800000
    quadtree::extract(qt, pts, extents = TRUE)
    #>        id  xmin  xmax  ymin  ymax     value
    #> [1,] 4315     0  8000  8000 16000       NaN
    #> [2,] 1814  2000  4000 28000 30000 0.1451094
    #> [3,] 3378 20000 21000 31000 32000 0.8202500
    #> [4,] 4545 24000 32000  8000 16000 0.9556426
    #> [5,] 1555 32500 33000 36000 36500 0.6800000

### Retrieving neighbor relationships

In addition to extracting values from cells, it is possible to determine neighbor relationships between cells via the `get_neighbors()` function. Given a point, `get_neighbors()` returns the cells that neighbor the cell that the point falls in. Note that cells that are diagonal from each either (i.e. only touch at the corner) are considered to be neighbors.

    get_neighbors(qt, as.numeric(pts[3,]))
    #>        id  xmin  xmax  ymin  ymax     value
    #> [1,] 1038 20000 21000 32000 33000 0.8236875
    #> [2,] 1039 21000 22000 32000 33000 0.7893750
    #> [3,]  967 19500 20000 32000 32500 0.8557500
    #> [4,] 3380 20000 21000 30000 31000 0.8688750
    #> [5,] 3381 21000 22000 30000 31000 0.8125625
    #> [6,] 3379 21000 22000 31000 32000 0.8061875
    #> [7,] 3279 19000 20000 31000 32000 0.9098125
    #> [8,] 3285 19000 20000 30000 31000 0.9466250

## Copying a quadtree

A key concept to understand is that `Quadtree` objects are pointers to C++ objects, and copying these objects simply by assigning them to a new variable makes a shallow copy rather than a deep copy. This is important to note because this differs from how R objects normally behave. For example, if I have a data frame named `df1`, I can make a copy of it by doing `df2 <- df1`. This makes a deep copy of `df1` - if I were to modify `df2`, `df1` would remain unchanged. This is not the case with `Quadtree` objects. The R object that users interact with is a pointer to a C++ object. If I were to attempt to make copy of a quadtree `qt1` by doing `qt2 <- qt1`, this simply copies the pointer to the C++ object, so `qt1` and `qt2` point to the same object. Thus, if `qt2` is modified (for example, by using a function like `set_values()` or `transform_values()`), `qt1` will also be modified.

The `copy()` function can be used to make a deep copy of a `Quadtree` object. This is useful if the user wishes to make changes to a quadtree without modifying the original.

    qt_copy <- copy(qt)
    qt
    #> class         : Quadtree
    #> # of cells    : 4189
    #> min cell size : 250
    #> extent        : 0, 64000, 0, 64000 (xmin, xmax, ymin, ymax)
    #> crs           : PROJCRS["unknown",
    #>     BASEGEOGCRS["unknown",
    #>         DATUM["Unknown based on Airy 1830 ellipsoid",
    #>             ELLIPSOID["Airy 1830",6377563.396,299.3249646,
    #>                 LENGTHUNIT["metre",1,
    #>                     ID["EPSG",9001]]]],
    #>         PRIMEM["Greenwich",0,
    #>             ANGLEUNIT["degree",0.0174532925199433,
    #>                 ID["EPSG",9122]]]],
    #>     CONVERSION["Transverse Mercator",
    #>         METHOD["Transverse Mercator",
    #>             ID["EPSG",9807]],
    #>         PARAMETER["Latitude of natural origin",49,
    #>             ANGLEUNIT["degree",0.0174532925199433],
    #>             ID["EPSG",8801]],
    #>         PARAMETER["Longitude of natural origin",-2,
    #>             ANGLEUNIT["degree",0.0174532925199433],
    #>             ID["EPSG",8802]],
    #>         PARAMETER["Scale factor at natural origin",0.9996012717,
    #>             SCALEUNIT["unity",1],
    #>             ID["EPSG",8805]],
    #>         PARAMETER["False easting",400000,
    #>             LENGTHUNIT["metre",1],
    #>             ID["EPSG",8806]],
    #>         PARAMETER["False northing",-100000,
    #>             LENGTHUNIT["metre",1],
    #>             ID["EPSG",8807]]],
    #>     CS[Cartesian,2],
    #>         AXIS["easting",east,
    #>             ORDER[1],
    #>             LENGTHUNIT["metre",1,
    #>                 ID["EPSG",9001]]],
    #>         AXIS["northing",north,
    #>             ORDER[2],
    #>             LENGTHUNIT["metre",1,
    #>                 ID["EPSG",9001]]]]
    #> values        : 0.109, 0.974 (min, max)
    qt_copy
    #> class         : Quadtree
    #> # of cells    : 4189
    #> min cell size : 250
    #> extent        : 0, 64000, 0, 64000 (xmin, xmax, ymin, ymax)
    #> crs           : PROJCRS["unknown",
    #>     BASEGEOGCRS["unknown",
    #>         DATUM["Unknown based on Airy 1830 ellipsoid",
    #>             ELLIPSOID["Airy 1830",6377563.396,299.3249646,
    #>                 LENGTHUNIT["metre",1,
    #>                     ID["EPSG",9001]]]],
    #>         PRIMEM["Greenwich",0,
    #>             ANGLEUNIT["degree",0.0174532925199433,
    #>                 ID["EPSG",9122]]]],
    #>     CONVERSION["Transverse Mercator",
    #>         METHOD["Transverse Mercator",
    #>             ID["EPSG",9807]],
    #>         PARAMETER["Latitude of natural origin",49,
    #>             ANGLEUNIT["degree",0.0174532925199433],
    #>             ID["EPSG",8801]],
    #>         PARAMETER["Longitude of natural origin",-2,
    #>             ANGLEUNIT["degree",0.0174532925199433],
    #>             ID["EPSG",8802]],
    #>         PARAMETER["Scale factor at natural origin",0.9996012717,
    #>             SCALEUNIT["unity",1],
    #>             ID["EPSG",8805]],
    #>         PARAMETER["False easting",400000,
    #>             LENGTHUNIT["metre",1],
    #>             ID["EPSG",8806]],
    #>         PARAMETER["False northing",-100000,
    #>             LENGTHUNIT["metre",1],
    #>             ID["EPSG",8807]]],
    #>     CS[Cartesian,2],
    #>         AXIS["easting",east,
    #>             ORDER[1],
    #>             LENGTHUNIT["metre",1,
    #>                 ID["EPSG",9001]]],
    #>         AXIS["northing",north,
    #>             ORDER[2],
    #>             LENGTHUNIT["metre",1,
    #>                 ID["EPSG",9001]]]]
    #> values        : 0.109, 0.974 (min, max)

In the next section (“Modifying cell values”), this function is used in the examples to create copies of a quadtree before modifying it so that the original remains unchanged.

## Modifying cell values

Two functions are provided for modifying cell values: `set_values()` and `transform_values()`.

### Setting values of cells

`set_values()` takes a matrix of points and associated values and changes the values of the cells the points fall in. Note that only the values of the terminal cells that the point falls in are changed.

    qt2 <- copy(qt)
    set_values(qt2, pts, rep(2, nrow(pts)))
    plot(qt2, crop = TRUE, border_lwd = .3)

### Transforming all cell values

`transform_values()` uses a function to change all cell values. For example, we could use this function to cube all of the cell values:

    qt3 <- copy(qt)
    transform_values(qt3, function(x) x^3)

    par(mfrow = c(1,2))
    plot(qt, crop = TRUE, na_col = NULL, border_lwd = .3, zlim = c(0, 1),
         legend = FALSE, main = "original quadtree")
    plot(qt3, crop = TRUE, na_col = NULL, border_lwd = .3, zlim = c(0, 1),
         legend = FALSE, main = "values cubed")

## Reading and writing quadtrees

One of the disadvantages of using C++ classes via Rcpp is that the objects are not preserved between sessions. For example, we could save a `Quadtree` to a `RData` file using `save()`, but when we loaded it back in we would get an error.

    qt_temp <- copy(qt)
    filepath <- tempfile()
    save(qt_temp, file = filepath)
    load(filepath)
    qt_temp
    #> Error in (function (cond) : error in evaluating the argument 'x' in selecting a method for function 'ext': external pointer is not valid

The `quadtree` package provides read and write functionality via the `cereal` C++ library, which is used to serialize C++ objects. `read_quadtree()` and `write_quadtree()` can be used to read and write `Quadtree` objects.

    filepath <- tempfile()
    write_quadtree(filepath, qt)
    qt_read <- read_quadtree(filepath)
    qt_read
    #> class         : Quadtree
    #> # of cells    : 4189
    #> min cell size : 250
    #> extent        : 0, 64000, 0, 64000 (xmin, xmax, ymin, ymax)
    #> crs           : PROJCRS["unknown",
    #>     BASEGEOGCRS["unknown",
    #>         DATUM["Unknown based on Airy 1830 ellipsoid",
    #>             ELLIPSOID["Airy 1830",6377563.396,299.3249646,
    #>                 LENGTHUNIT["metre",1,
    #>                     ID["EPSG",9001]]]],
    #>         PRIMEM["Greenwich",0,
    #>             ANGLEUNIT["degree",0.0174532925199433,
    #>                 ID["EPSG",9122]]]],
    #>     CONVERSION["Transverse Mercator",
    #>         METHOD["Transverse Mercator",
    #>             ID["EPSG",9807]],
    #>         PARAMETER["Latitude of natural origin",49,
    #>             ANGLEUNIT["degree",0.0174532925199433],
    #>             ID["EPSG",8801]],
    #>         PARAMETER["Longitude of natural origin",-2,
    #>             ANGLEUNIT["degree",0.0174532925199433],
    #>             ID["EPSG",8802]],
    #>         PARAMETER["Scale factor at natural origin",0.9996012717,
    #>             SCALEUNIT["unity",1],
    #>             ID["EPSG",8805]],
    #>         PARAMETER["False easting",400000,
    #>             LENGTHUNIT["metre",1],
    #>             ID["EPSG",8806]],
    #>         PARAMETER["False northing",-100000,
    #>             LENGTHUNIT["metre",1],
    #>             ID["EPSG",8807]]],
    #>     CS[Cartesian,2],
    #>         AXIS["easting",east,
    #>             ORDER[1],
    #>             LENGTHUNIT["metre",1,
    #>                 ID["EPSG",9001]]],
    #>         AXIS["northing",north,
    #>             ORDER[2],
    #>             LENGTHUNIT["metre",1,
    #>                 ID["EPSG",9001]]]]
    #> values        : 0.109, 0.974 (min, max)

## Converting quadtrees to other data types

Functions are also provided for converting a `Quadtree` to different data types.

### Get all cell values as a vector

`as_vector()` returns the cell values as a numeric vector. It has one optional parameter - `terminal_only`. If `TRUE` (the default), only the values of the terminal cells are returned. If `FALSE`, the values of all cells (including those with children) are returned.

    vec1 <- as_vector(qt)
    length(vec1)
    #> [1] 4189
    summary(vec1)
    #>    Min. 1st Qu.  Median    Mean 3rd Qu.    Max.    NA's
    #>  0.1090  0.4843  0.6279  0.6364  0.8142  0.9740     832

    vec2 <- as_vector(qt, FALSE)
    length(vec2)
    #> [1] 5585
    summary(vec2)
    #>    Min. 1st Qu.  Median    Mean 3rd Qu.    Max.    NA's
    #>  0.1090  0.4910  0.6330  0.6447  0.8320  0.9740     832

### Convert to a data frame

`as_data_frame()` converts the quadtree to a data frame. In the resulting data frame, each row represents a single cell and has columns describing the cell ID, x and y limits, value, level (i.e. node depth), the size of its smallest child, the ID of its parent, and whether it has children. As with `as_vector()`, there is one optional parameter called `terminal_only` - when `TRUE` (the default) only the terminal cells are returned, and when `FALSE`, all cells are returned.

    df1 <- as_data_frame(qt)
    dim(df1)
    #> [1] 4189   10
    head(df1)
    #>    id hasChildren level  xmin  xmax  ymin  ymax value smallestChildLength
    #> 3   2           0     2     0 16000 48000 64000   NaN               16000
    #> 4   3           0     2 16000 32000 48000 64000   NaN               16000
    #> 7   6           0     4     0  4000 44000 48000   NaN                4000
    #> 9   8           0     5  4000  6000 46000 48000   NaN                2000
    #> 10  9           0     5  6000  8000 46000 48000   NaN                2000
    #> 11 10           0     5  4000  6000 44000 46000   NaN                2000
    #>    parentID
    #> 3         1
    #> 4         1
    #> 7         5
    #> 9         7
    #> 10        7
    #> 11        7

    df2 <- as_data_frame(qt, FALSE)
    dim(df2)
    #> [1] 5585   10
    head(df2)
    #>   id hasChildren level  xmin  xmax  ymin  ymax     value smallestChildLength
    #> 1  0           1     0     0 64000     0 64000 0.7559473                 250
    #> 2  1           1     1     0 32000 32000 64000 0.7309282                 250
    #> 3  2           0     2     0 16000 48000 64000       NaN               16000
    #> 4  3           0     2 16000 32000 48000 64000       NaN               16000
    #> 5  4           1     2     0 16000 32000 48000 0.8131046                 250
    #> 6  5           1     3     0  8000 40000 48000 0.9623516                 250
    #>   parentID
    #> 1       -1
    #> 2        0
    #> 3        1
    #> 4        1
    #> 5        1
    #> 6        4

### Convert to a raster

A `Quadtree` can also be converted to a raster using the `as_raster()` function. This function has one optional parameter named `rast`, which is used as the template for the output raster. If `NULL` (the default), a raster is automatically created, where the quadtree extent is used as the raster extent and the smallest cell in the quadtree is used to determine the resolution of the raster. Note that the value of a raster cell is determined by the value of the quadtree cell located at the centroid of the raster cell - thus, if a raster cell overlaps several quadtree cells, whichever quadtree cell the centroid of the raster cell falls in will determine the raster cell’s value. If the default value for `rast` is used, the raster cells will never overlap with more than one quadtree cell.

    rast <- as_raster(qt)
    plot(rast)