title: "Spatial data handling"

output:

rmdshower::shower_presentation:

theme: material

```{r, echo=F,message=F,warning=F}

Packages:

* `sp` First major spatial data package/format

* `rgdal` reading and writing spatial data

* `rgeos` Interface to open-source geometry engine (GEOS)

* `sf` Spatial Features in the 'tidyverse'

* `raster` raster (gridded) data

* and a few others...

Typically an object in the real world, such as a building or a tree.

A forest stand can be a feature, a forest can be a feature, a city can be a feature.

A satellite image pixel can be a feature, a complete image can be

a feature too.

```{r, echo=F}

What information do we need to store in order to define points, lines, polygons in geographic space?

- lat/lon coordinates

- projection

- what type (point/line/poly)

- if polygon, is it a hole or not

- attribute data

- ... ?

Features have a _geometry_ describing _where_ on Earth the

feature is located, and they have attributes, which describe other

properties.

A Tree:

* delineation of its crown

* its stem

* point indicating its centre

Attributes:

* species,

* height

* diameter at breast height at a particular date, and so on.

"_A simple feature is defined by the OpenGIS Abstract specification to have both spatial and non-spatial attributes. Spatial attributes are geometry valued, and simple features are based on 2D geometry with linear interpolation between vertices._"

All geometries composed of points. Points are coordinates in a

2-, 3- or 4-dimensional space. All points in a geometry have the

same dimensionality. In addition to X and Y coordinates, there are

two optional additional dimensions:

* a Z coordinate, denoting altitude

* an M coordinate (rarely used), denoting some _measure_ that is associated with the point, rather than with the feature as a whole (in which case it would be a feature attribute); examples could be time of measurement, or measurement error of the coordinates

The four possible cases then are:

1. two-dimensional points refer to x and y, easting and northing, or longitude and latitude, we refer to them as XY

2. three-dimensional points as XYZ

3. three-dimensional points as XYM

4. four-dimensional points as XYZM (the third axis is Z, fourth M)

The following seven simple feature types are the most common, and are for instance the only ones used for [GeoJSON](https://tools.ietf.org/html/rfc7946):

| type | description |

| ---- | ----------- |

| `POINT` | a single point |

| `LINESTRING` | sequence of points connected by lines |

| `POLYGON` | sequence of points form a closed ring |

| `MULTIPOINT` | set of points |

| `MULTILINESTRING` | set of linestrings |

| `MULTIPOLYGON` | set of polygons |

| `GEOMETRYCOLLECTION` | set of geometries |

10 more geometries 10 are rare:

* `CIRCULARSTRING`

* `COMPOUNDCURVE`

* `CURVEPOLYGON`

* `MULTICURVE`

* `MULTISURFACE`

* `CURVE`

* `SURFACE`

* `POLYHEDRALSURFACE`

* `TIN`

* `TRIANGLE`

Coordinates can only be placed on the Earth's surface when their coordinate

reference system (CRS) is known; this may be an elipsoidal CRS such as WGS84, a projected, two-dimensional (Cartesian) CRS such as a UTM zone or Web Mercator, or a CRS

in three-dimensions or [including time](http://www.faculty.jacobs-university.de/pbaumann/iu-bremen.de_pbaumann/Papers/acmgis-2012_crs-nts.pdf). Similarly, M-coordinates need an attribute reference

system, e.g. a [measurement unit](https://CRAN.R-project.org/package=units).

Package `sf` represents simple features as native R objects.

Similar to [PostGIS](http://postgis.net/), all functions and methods

in `sf` that operate on spatial data are prefixed by `st_`, which

refers to _spatial and temporal_; this makes them easily findable

by command-line completion. Simple features are implemented as

R native data, using simple data structures (S3 classes, lists,

matrix, vector). Typical use involves reading, manipulating and

writing of sets of features, with attributes and geometries.

Attributes typically stored in `data.frame` objects (or the

very similar `tbl_df`), we will also store feature geometries in

a `data.frame` column. Since geometries are not single-valued,

they are put in a list-column, a list of length equal to the number

of records in the `data.frame`, with each list element holding the simple

feature geometry of that feature. The three classes used to

represent simple features are:

* `sf`, the table (`data.frame`) with feature attributes and feature geometries, which contains

* `sfc`, the list-column with the geometries for each feature (record), which is composed of

* `sfg`, the feature geometry of an individual simple feature.

There are currently two main approaches in R to handle geographic vector data.

First package to provide classes and methods for spatial data types in R is called [`sp`](https://cran.r-project.org/package=sp). Provides classes and methods to create _points_, _lines_, _polygons_, and _grids_ and to operate on them.

About 350 of the spatial analysis packages use the spatial data types that are implemented in `sp` i.e. they "depend" on the `sp` package and many more are indirectly dependent.

Many spatial R packages still depends on the **sp** package, therefore it is important to know how to convert **sp** to and from **sf** objects

Foundational structure for any spatial object in `sp` is the `Spatial` class. It has two "slots" ([new-style S4 class objects in R have pre-defined components called slots](http://stackoverflow.com/a/4714080)):

* a __bounding box__

* a __CRS class object__ to define the Coordinate Reference System

Basic structure extended depending on the characteristics of the spatial object (point, line, polygon).

* __Points__ most basic spatial data objects. Typically from a two-column matrix/dataframe with a column for latitude and one for longitude.

* __Lines__ are generated out of `Line` objects. A `Line` object is a collection of 2D coordinates and is generated out of a two-column matrix or a dataframe with a column for latitude and one for longitude. A `Lines` object is a __list__ of one or more `Line` objects.

* __Polygons__ are generated out of `Polygon` objects. A `Polygon` object is a collection of 2D coordinates with equal first and last coordinates and is generated out of a two-column matrix or a dataframe with a column for latitude and one for longitude. A `Polygons` object is a __list__ of one or more `Polygon` objects, for example islands belonging to the same country.

A very simple `Line` object:

A very simple `Lines` object:

Adds bounding box and the slot for the Coordinate Reference System or CRS. `SpatialPoints` can be directly generated out of the coordinates. `SpatialLines` and `SpatialPolygons` objects are generated using lists of `Lines` or `Polygons` objects respectively (more below).

See here for how to create a `SpatialLines` object:

Add attribute data, which will turn your `Spatial*` object into a `Spatial*DataFrame` object. ID fields are here required to match the data frame row names.

Create a simple `SpatialLinesDataframe`:

function | and what it does

------------ | ------------------------------------------------------

`bbox()` | returns the bounding box coordinates

`proj4string()` | sets or retrieves projection attributes using the CRS object.

`CRS()` | creates an object of class of coordinate reference system arguments

`spplot()` | plots all the attributes

`coordinates()` | set or retrieve the spatial coordinates.

`over(a, b)` | used to retrieve the polygon or grid indices on a set of points

`spsample()` | sampling of spatial points within the spatial extent of objects

[`sf`](https://cran.r-project.org/package=sf) implements a formal standard called ["Simple Features"](https://en.wikipedia.org/wiki/Simple_Features) that specifies a storage and access model of spatial geometries (point, line, polygon). A feature geometry is called simple when it consists of points connected by straight line pieces, and does not intersect itself.

This standard has been adopted widely, not only by spatial databases such as PostGIS, but also more recent standards such as GeoJSON.

If you work with PostGis or GeoJSON you may have come across the [WKT (well-known text)](https://en.wikipedia.org/wiki/Well-known_text) format, for example like these:

`POINT (30 10)

LINESTRING (30 10, 10 30, 40 40)

POLYGON ((30 10, 40 40, 20 40, 10 20, 30 10))`

`sf` implements this standard natively in R. Data are structured and conceptualized very differently from the `sp` approach.

In `sf` spatial objects are stored as a simple data frame with a special column that contains the information for the geographic coordinates. That special column is a list with the same length as the number of rows in the data frame. Each of the individual list elements then can be of any length needed to hold the coordinates that correspond to an individual feature.

Geometric objects (simple features) can be created from a numeric vector, matrix or a list with the coordinates. They are called `sfg` objects for Simple Feature Geometry.

See here for an example of how a LINESTRING `sfg` object is created:

In order to work our way towards a data frame for all features we create what is called an `sfc` object with all individual features, which stands for Simple Feature Collection. The `sfc` object also holds the bounding box and the projection information.

See here for an example of how a `sfc` object is created:

We now combine the dataframe with the attributes and the simple feature collection.

See here how its done.

There are many methods available in the `sf` package, to find out use `methods(class="sp")`

* provides **fast** I/O, particularly relevant for large files

* directly reads from and writes to spatial **databases** such as PostGIS

* compatibile with the *tidyverse*

* stay tuned for a new `ggplot` release that will be able to read and plot the `sf` format without the need of conversion to a data frame, like the `sp` format

`sp` and `sf` are _not_ only formats for spatial objects. Other spatial packages may use their own class definitions for spatial data (for example `spatstat`). Usuallly you can find functions that convert `sp` and increasingly `sf` objects to and from these formats.

- The structures in the **sp** packages are more complicated - `str(world_sf)` vs `str(world_sp)`

- Moreover, many of the **sp** functions are not "pipeable" (it's hard to combine them with the **tidyverse**)

```{r, eval = F}

`Error in UseMethod("filter_") :

no applicable method for 'filter_' applied to an object of class "c('SpatialPolygonsDataFrame', 'SpatialPolygons', 'Spatial', 'SpatialPolygonsNULL', 'SpatialVector')"`

Counterpart to `st_read()` is the `st_write` function, e.g. `st_write(world, 'data/new_world.gpkg')`. A full list of supported formats could be found using `sf::st_drivers()`.

```{r, eval = FALSE}

```{r, echo = FALSE}

```{r, eval=FALSE}

```{r, echo=FALSE}

```{r, eval=FALSE}

```{r, echo=FALSE}

```{r, warning=FALSE}

```{r, echo=FALSE}

The `st_set_geometry` function can be used to remove the geometry column:

It's a big topic which includes:

- Spatial subsetting

- Spatial joining/aggregation

- Topological relations

- Distances

- Spatial geometry modification

- Raster operations (map algebra)

See [Chapter 4](http://robinlovelace.net/geocompr/spatial-data-operations.html#spatial-operations-on-raster-data) of *Geocomputation with R*

Transform (warp) to a different projection:

```{r, warning = FALSE, message = FALSE, fig.height = 4}

- Basic maps of `sf` objects can be quickly created using the `plot()` function:

Syntax for creating plots is similar to that of ggplot2.

The add-on package `tmaptools` contains tool functions for reading and processing shape files.

```{r, fig.align='center', fig.height=4, message=FALSE}

```{r, eval=F}

```{r, echo=T, message=T,results='asis'}

<iframe src="04_assets/leaflet.html" width=100% height=400px></iframe>

Raster data is not yet closely connected to the **tidyverse**, however:

- Some functions from the **raster** package works well in `pipes`

- You can convert vector data to the `Spatial*` form using `as(my_vector, "Spatial")`before working on raster-vector interactions

- There are some initial efforts to bring raster data closer to the **tidyverse**, such as [tabularaster](https://github.com/hypertidy/tabularaster), [sfraster](https://github.com/mdsumner/sfraster), or [fasterize](https://github.com/ecohealthalliance/fasterize)

- The development of the [stars](https://github.com/r-spatial/stars), package focused on multidimensional, large datasets should start soon. It will allow pipe-based workflows

- Slides adapted from:

- "Robin Lovelace and Jakub Nowosad" draft book [_Geocomputation with R_ (to be published in 2018)](http://robinlovelace.net/geocompr/). Source code at https://github.com/robinlovelace/geocompr.

- [Claudia Engel's spatial analysis workshop](https://github.com/cengel/rspatial/blob/master/2_spDataTypes.Rmd)