README.Rmd

---
output: github_document
---

<!-- README.md is generated from README.Rmd. Please edit that file -->

```{r, include = FALSE}
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>",
  fig.path = "man/figures/README-",
  out.width = "100%"
)
```

# sdf

<!-- badges: start -->
[![Lifecycle: experimental](https://img.shields.io/badge/lifecycle-experimental-orange.svg)](https://lifecycle.r-lib.org/articles/stages.html#experimental)
<!--[![CRAN status](https://www.r-pkg.org/badges/version/sdf)](https://CRAN.R-project.org/package=sdf)-->
<!-- badges: end -->

The goal of `{sdf}` is to provide an extensible spatial data frame. sdf acts as a common interface to a variety of different geometry types. sdf is inspired in part by dplyr where the dplyr syntax acts as an interface that can be utilized by a number of different of data sources without changing the syntax—e.g. dbplyr, sparklyr, dtplyr, among others. 

It is inspired by Rust traits where if you implement the correct generics, you get a ton of sweet benefits. By implementing the required generics your geometry type can become an `sdf` (spatial data frame).

If you don't know what is meant by "generic" [start here](https://adv-r.hadley.nz/s3.html).

## Installation

You can install the development version of sdf like so:

``` r
remotes::install_github("josiahparry/sdf")
```

## Example

The `sdf` generics have been implemented for the `sf`, `geos`, and `s2` packages. `sdf` will work with those three packages out of the box.


First lets grab an sf object to try converting from. 

```{r}
library(sf) # so sf methods work
library(sdf)

# pull some spatial data
data(guerry, package = "sfdep")
g <- dplyr::select(guerry, code_dept, crime_pers, region)
```

We can then cast the sf object to a spatial data frame. 

```{r}
# using sf, cast to a spatial data frame
sdf_sf <- as_sdf(g)

sdf_sf
```


And the same thing goes for geos. 

```{r}
# using geos
sdf_geos <- as_sdf(
  # convert geometry to rust geometry
  dplyr::mutate(g, geometry = geos::as_geos_geometry(geometry))
)

sdf_geos
```

The spatial data frame is compatible with dplyr out of the box and it doesn't care about your geometry class. 

```{r}
library(dplyr, warn.conflicts = FALSE)

sdf_sf |> 
  group_by(region) |> 
  summarise(total_crime = sum(crime_pers))

sdf_geos |> 
  group_by(region) |> 
  summarise(total_crime = sum(crime_pers))
```

## Philosophy

There are a number of different geometry types among low-level libraries. Geometry types should be interchangeable and the type you use should not be dictated by the libraries you use. All geometries should be able to convert to and from one another. Some libraries have better performance, accuracy, or algorithms that can and should be used when appropriate. 

Geospatial and data science more broadly are quickly moving towards a unified memory format à la arrow. Being able to represent our geometries in a common memory format enables us to convert to and from anything. Today, in R, that approach is handled by `{wk}`. In the future, this will be `{geoarrow}`. 

This package is intended to serve as an example of what a spatial data frame of the future might look like. The geometry column is independent of the R package itself. All that is required is a number of s3 generic functions to be implemented. 

This is how I hope to see the `{sf}` package. I want to see it as an enabler of a variety of geometry backends some of which may even be lazily evaluated (such as a PostGIS backend).


## Required method implementations

The following generics must have an implementation for your geometry class.
The geometry class must be an S3 object / compatible with a data frame.

- `c()` 
- `is_geometry()`
- `bounding_box()`
- `combine_geometry()`

Implementing these generics provides access to dplyr functionality only. 

### Default implementations 

Since the `wk` package acts as a middle-man between memory formats, we can use it for default implementations. Any object which is identified as a vector by `vctrs` and is handleable by `wk` will be minimally functional.

- `is_geometry(x)` has a default implementation which ensures that `x` is a vector and is useable by `wk` via `vctrs::vec_is()` and `wk::is_handleable()`. 

- `bounding_box()` has default implementation of `wk::wk_bbox()`

- `combine_geometry()` has a default implementation of `wk::wk_collection()`

## Spatial joins and filters

The functions `sdf_join()` and `sdf_filter()` have an argument `predicate` which takes a predicate function—by default `sdf_intersects()`. For two spatial data frames `x` and `y` the predicate must return a list with the same length as `x`. Each element must be a list of integer vectors where each element indicates the corresponding elements in `y` where the predicate function is true. 

The following predicate generic functions are exported:

- `sdf_intersects()`
- `sdf_contains()`
- `sdf_within()`
- `sdf_crosses()` 
- `sdf_covers()`
- `sdf_covered_by()`
- `sdf_equals()`
- `sdf_disjoint()`
- `sdf_touches()`
- `sdf_overlaps()`

It is recommended that you implement a few of generics, though not mandatory. Any function that returns a sparse matrix as described above for two geometry columns can be used in `sdf_join()` and `sdf_filter()`

## Other generics 

There are a number of other generics that are exports. The following generics are considered to be part of the "main" interface. These are the functions that are strongly recommended to implement for your geometry vector class. 

- `union_geometry()`
- `simplify_geometry()`
- `buffer_geometry()`
- `centroid()`
- `convex_hull()`
- `concave_hull()`

## Implementing your geometry class

`sdf` provides templates for you to implement your own geometry class. There are two components:

1. writing the methods for your class and
2. registering your method with `sdf`.

To this end there are two templates for you. `use_sdf_template_class()` will open a new R script with the scaffolding needed to implement your geometry class. 

The first 30 lines of the template are below to illustrate what the template looks like:

```{r, echo = FALSE}
cat(readLines(system.file("templates/new-class-compat.R", package = "sdf"))[1:30], sep = "\n")
```
Secondly the function `use_sdf_template_zzz()` will create the required `.onLoad()` function call in a `zzz.R` file for you to register your methods and class. 

The first 15 lines are below for an example.

```{r, echo = FALSE}
cat(readLines(system.file("templates/zzz.R", package = "sdf"))[1:15], sep = "\n")
```


### Implementing `{geos}` generics

While there already is a geos implementation, it is worth walking through how to do so. 

```{r}
library(geos)
geo <-  as_geos_geometry(guerry$geometry)
```

The first step is to define the `is_geometry()` generic which must return `TRUE` for your geometry class. 

```{r}
is_geometry.geos_geometry <- function(x) inherits(x, "geos_geometry")

is_geometry(geo)
```

Next, we need to implement the bounding box generic. There is not any bounding box function in geos that I am aware of so we have to implement our own. The `bounding_box()` generic _must_ return a numeric vector of length 4 with the names `xmin`, `ymin`, `xmax`, and `ymax` in that order. The bounding box is checked with the internal function `validate_bbox()`.

```{r}
bounding_box.geos_geometry <- function(x) {
  extents <- geos::geos_extent(x)
  apply(extents, 2, min)
}

bounding_box(geo)
```

For minimal functionality, we need to lastly provide a `combine_geometry()` method. This is used in the `summarise()` function. `combine_geometry()` is intended to combine single geometries into the multi variety or, at minimum, a geometry collection as done here. 

```{r}
combine_geometry.geos_geometry <- function(x) {
  geos::geos_make_collection(x)
}

combine_geometry(geo)
```

Now we can create an sdf object with only this alone.

```{r}
library(sf)
library(dplyr)

geosdf <- as_tibble(guerry) |> 
  mutate(geometry = as_geos_geometry(geometry)) |> 
  as_sdf()

geosdf
```

With only those three methods we can now use dplyr functionality. 

```{r}
geosdf |> 
  group_by(code_dept) |> 
  summarise(total_crime = sum(crime_pers))
```

This alone is super useful! But we may also want to be able to implement spatial joins and filters. To do this we need to implement some predicate functions methods. In this case we will create only the method for `sdf_intersects()`. 

The below function creates a GEOS STRTree which is a type of search tree that makes speeds up join candidate searches. Then we use the tree to actually find the intersecting neighbors.

```{r}
sdf_intersects.geos_geometry <- function(x, y, ...) {
  rtree <- geos::geos_strtree(y)
  geos::geos_intersects_matrix(x, rtree)
}
```

To use this function let's sample 25 random points from the original sf geometry and join on the region data to it. 

```{r}
pnt_sample <- sf::st_sample(guerry, 25)

pnts <- as_sdf(
  list(geometry = as_geos_geometry(pnt_sample))
)

sdf_join(pnts, geosdf)
```