06-transform.Rmd

# Reprojections and Transformations {#transform}

## Prerequisites {-}

- This chapter requires the packages **tidyverse**, **sf**, and **spData**:

```{r, message=FALSE}
library(tidyverse)
library(sf)
library(spData)
```

## Introduction

<!-- random notes -->
<!-- CRS also represent spatial relationship between datasets. -->
<!-- Therefore, spatial operations could only be correctly performed on data with the same CRS. -->

<!-- While the concept of CRS applies to both data types, conversion between coordinate rerefence systems differs between vector and raster. -->
<!-- Tranformation from one CRS to another in vector data changes only coordinates of vertices, keeping the values intact. -->
<!-- (for most of the case is better to reproject vector than raster) -->
<!-- rasters: transformation means change of the coordinates of (special case of resampling) -->
<!-- changes in dimensions, resolution, extent -->
<!-- change shape and attributes) -->
<!-- different methods of computing values after transformation, such as ngb or bilinear  -->

<!-- - crs - a heart of spatial data -->
<!-- - short history  -->
<!-- - types of crs (geographic vs cartesian; local vs regional vs global) -->
<!-- - objectives - 1/ to combine different datasets, 2/ area calculations, 3/ distance mesasurement, 4/ navigation, 5/ spatial data representations -->
<!-- - proj -->
<!-- - proj4 + epsg -->
<!-- - the most popular epsg -->
<!--in most of the cases reproject vector, not raster-->

As stated in Chapter \@ref(crs-intro), it is important to understand which CRS you are working in when undertaking spatial operations.
Many spatial operations assume that you are using a *projected* CRS (on a Euclidean grid with units of meters rather than a geographic 'lat/lon' grid with units of degrees).
The GEOS engine underlying most spatial operations in **sf**, for example, assume your data is in a projected CRS.
For this reason **sf** contains a function for checking if geometries have a geographic or projected CRS.
This is illustrated below using the example of the *Greenwich point* which came to define 0 degrees longitude:

```{r}
greenwich = st_sf(geometry = st_sfc(st_point(c(0, 51.5))))
st_is_longlat(greenwich)
```

The results show that when geographic data is created from scratch, or is loaded from a source that has no CRS metadata, the CRS is unspecified by default.
Spatial operations on objects without a CRS run on the implicit assumption that they are projected, even when in reality they are not.
This can be seen by creating a huge buffer of 10 degrees around the `greenwich` point:

```{r}
greenwich_buff = st_buffer(greenwich, dist = 10)
```

Brief consideration of what has happened should set alarm bells ringing:
the buffer will be highly distorted when projected onto the surface of the Earth.
For this reason, when **sf** (and other spatial packages) know that geometries are in lat/long coordinates they emit a warning.
This is illustrated in the code below which sets the CRS of `greenwich` to a lat/lon CRS (the commonly used EPSG 4326 in this case), checks to ensure that R thinks it's a geographic CRS, and then re-applies the buffer:

```{r}
greenwich_latlon = st_set_crs(greenwich, 4326)
st_is_longlat(greenwich_latlon)
greenwich_buff_latlon = st_buffer(greenwich_latlon, 10)
```

The results show that, as expected, `4326` is a geographic (lat/lon) CRS.
As a result a warning message is emitted to warn the user that the operation may not work correctly and that, if the operation was intended, the distance should be in degrees (not meters or some other Euclidean distance measurement).
The seemingly small difference in setting the CRS may seem inconsequential but it can have a huge impact.
This is illustrated in Figure \@ref(fig:crs-buf), which shows how the two buffers are plotted, with the (correctly) defined CRS of the latter object being dramatically elongated in the north-south direction due to the thinning of the vertical lines of longitude towards the Earth's poles.  

```{r, echo=FALSE}
par(mfrow = c(1, 2))
```

```{r crs-buf, fig.cap="Buffers on geographic data with undefined (top) and defined (bottom) CRSs."}
plot(greenwich_buff, graticule = st_crs(4326))
plot(greenwich_buff_latlon, graticule = st_crs(4326))
```

## CRS transformation

While CRSs can be set manually, it is more common in real world applications to *transform* a known CRS into another.
A typical example is when geometry data is provided in a geographic CRS but you want to do spatial operations, which require it to be in a projected CRS.
Let's use real-world examples to illustrate this.

### Vector data

The dataset `cycle_hire_osm` represents all cycle hire locations across London, take from OpenStreetMap (OSM).
It is automatically loaded by the **spData** package, meaning we do not have to load it, and its CRS can be queried as follows:

```{r}
st_crs(cycle_hire_osm)
```

Let's create a new version of it in a projected CRS, using the 'magic number' (a value to be explained subsequently) of 27700:

```{r}
cycle_hire_projected = st_transform(cycle_hire_osm, 27700)
st_crs(cycle_hire_projected)
```

Note that the result shows that the `epsg` has been updated and that `proj4string` element of the CRS now contains, among other things `+proj=tmerc` (meaning it is a projected CRS using the [tranverse Mercator](https://en.wikipedia.org/wiki/Transverse_Mercator_projection) projection) and `+units=m` (meaning the units of the coordinates are meters).
Another function, from the **rgdal** library, provides a note about the CRS (its name):

```{r}
crs_codes = rgdal::make_EPSG()[1:2]
dplyr::filter(crs_codes, code == 27700)
```

The result shows that the EPSG code 27700 represents the British National Grid, a result that could have been found by searching online for "[CRS 27700](https://www.google.com/search?q=CRS+27700)".
This projection is clearly inappropriate for the data: the coordinates represent degrees of longitude and latitude, and this can also be seen by plotting it over a basemap, e.g. with the **mapview** package: `mapview::mapview(sf_points)`.

The formula that converts a geographic point into a point on the surface of the Earth is provided by the `proj4string` element of the `crs` (see [proj4.org](http://proj4.org/) for further details):

```{r}
st_crs(27700)$proj4string
```

```{block2 type='rmdnote'}
The EPSG code can be found inside the `crs` attribute of the object's geometry.
It is hidden from view for most of the time except when the object is printed but can be can identified and set using the `st_crs` function, for example `st_crs(cycle_hire_osm)$epsg`.
```

<!--
- st_as_sf(x, coords = c("x","y"))
- st_crs(x)
- st_transform(x, crs)
- ==
- !st_is_longlat(x)
- st_set_crs(x, crs)
- st_proj_info("proj");st_proj_info("ellps");st_proj_info("datum");st_proj_info("units")
- st_bbox
- st_wrap_dateline
-->

### Raster data

<!--
- data? one numerical and one categorical
- projectRaster
- an issue of resampling (comparision of old and new values)
-->

<!-- ## Affine transformations -->

<!-- ### Translating -->

<!-- ### Scaling -->

<!-- ### Rotating -->

<!-- ### Reflecting -->

<!-- ### Shearing -->

<!-- Todo: add content on simplifying using mapshaper and other packages (e.g. sf) -->