title: "Spatial data"

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

Start thinking about:

* Question you want to answer

* Problem you want to solve

[<i class="fa fa-desktop"></i> If interested, you can download the R script associated with this presentation here](`r knitr::current_input()`).

Write out in a sentence what this code is doing. Make sure to catch the key points in your sentence

```{r eval=F}

* `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` gridded data (like satellite imagery)

* and a few others...

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

Features could include:

* a forest (polygon)

* a tree in the forest (point or polygon)

* a branch on the tree (line?)

* a complete image (multipoint, polygon, raster)

* a satellite image pixel of that forest (point or polgyon or raster)

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

* date of observation

* ...

"_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: coordinates in 2-, 3- or 4-dimensional space.

* **X** coordinate in X direction (typically longitude or similar)

* **Y** coordinate in Y direction (typically latitude or similar)

* **Z** coordinate denoting altitude

* **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.

The four possible cases then are:

1. 2D (XY): x and y, easting and northing, or longitude and latitude

2. 3D (XYZ): three-dimensional points

3. 3D (XYM): three-dimensional points where 3rd is some attribute space

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

| 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 |

Some formats only include these (e.g. [GeoJSON](https://tools.ietf.org/html/rfc7946))

10 more geometries 10 are rare:

* `CIRCULARSTRING`

* `COMPOUNDCURVE`

* `CURVEPOLYGON`

* `MULTICURVE`

* `MULTISURFACE`

* `CURVE`

* `SURFACE`

* `POLYHEDRALSURFACE`

* `TIN`

* `TRIANGLE`

SFs 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).

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

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

~350 of the spatial analysis packages use `sp` data types, so it is important to know how to convert **sp** to and from **sf** 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.

All spatial functions and methods in `sf` prefixed by `st_` (refering to _spatial and temporal_)

Simple features are implemented as R native data, using simple data structures (S3 classes, lists,

matrix, vector).

Stored as `data.frame` objects (or very similar `tbl_df`) with _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.

* `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.

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.

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.

Create a LINESTRING `sfg` object:

Create a `sfc` (Simple Feature Collection) object of individual features:

The `sfc` object also holds the bounding box and the projection information.

Add attributes (in a `data.frame`) to the `sfc` object to make a `sf` (Simple Features) object:

```{r, echo=F}

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

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

* compatibile with the *tidyverse*

* recent `ggplot` release can read and plot the `sf` format without conversion

`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`). Usually 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}

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=F}

```{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, message=F}

For more on `grid.arrange`, see [here](https://cran.r-project.org/web/packages/egg/vignettes/Ecosystem.html).

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

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

All the nice ggplot features are available

```{r, message=F}

[Leaflet](https://leafletjs.com/) leading open-source JavaScript library for mobile-friendly interactive maps. Weighing just about 38 KB of JS, it has all the mapping features most developers ever need.

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

<iframe id="test" style=" height:400px; width:100%;" scrolling="no" frameborder="0" src="world_leaflet.html"></iframe>

Page hosted on Github for free (except the domain name)...

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

- Some functions from `raster` work well in `pipes`

- Convert vector data to `Spatial*` form using `as(my_vector, "Spatial")` for raster-vector interactions

- Some early efforts to bring raster data into the **tidyverse**, including [tabularaster](https://github.com/hypertidy/tabularaster), [sfraster](https://github.com/mdsumner/sfraster), [fasterize](https://github.com/ecohealthalliance/fasterize), and [stars](https://github.com/r-spatial/stars) (multidimensional, large datasets).

- 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)