---
title: "crstools"
output: rmarkdown::html_vignette
vignette: >
  %\VignetteIndexEntry{crstools}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

```{r, include = FALSE}
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>"
)
```

`crstools` is a collection of tools to facilitate the choice of a Coordinate
Reference System (CRS) in R. CRSs are essential for mapping, as they define how
geographic data are projected onto a flat surface. The choice of CRS depends on
the type of map, the extent of the map, and the purpose of the map.
The `crstools` package provides 3 main functionalities:

1) `suggest_crs` provides suggestions for CRS for different types of maps, such
as equal-area maps, conformal maps, and equidistant maps, depending on the
extent and location of the area of interest.

2) `geom_tissot` visualises the distortion associated with a given crs, by
adding Tissot indicatrix to a `ggplot2` map (the distortion of circles on the map).

3) `add_crs` helps to add CRS to an image; it does so by defining Ground Control
Points (GCPs) and using GDAL to warp the image to a desired CRS.

## Installation

Currently, `crstools` is only available on GitHub. You can install it using the
`devtools` package:

```{r setup, eval=FALSE}
devtools::install_github("EvolEcolGroup/crstools")
```

## An overview of CRS codes

Coordinate Reference System (CRS) codes are standardized identifiers used to
define spatial reference systems, ensuring accurate geospatial data
representation. One common way to describe a CRS is using **Well-Known Text
(WKT)**, a human-readable format that specifies key parameters such as the
datum, projection, and coordinate units. Another widely used system is the
**EPSG code**, which is a numeric identifier assigned by the **European
Petroleum Survey Group (EPSG)** for commonly used CRS definitions, such as
**EPSG:4326** for WGS 84. Additionally, **Proj4** is a text-based format used by
the PROJ library, which provides concise parameter strings for defining
projections, transformations, and datums, making it particularly useful in GIS
software and programming environments. These different formats help ensure
interoperability between geospatial tools and datasets.

Let's take the **Web Mercator projection**, commonly used in web mapping
applications like Google Maps and OpenStreetMap. Below are its representations
in WKT, EPSG, and Proj4 formats:

1. **EPSG Code**  
   - **EPSG:3857** (also called "Pseudo-Mercator" or "Web Mercator")  

2. **Well-Known Text (WKT)**  
   ```wkt
   PROJCS["WGS 84 / Pseudo-Mercator",
      GEOGCS["WGS 84",
         DATUM["WGS_1984",
            SPHEROID["WGS 84",6378137,298.257223563,
               AUTHORITY["EPSG","7030"]],
            AUTHORITY["EPSG","6326"]],
         PRIMEM["Greenwich",0,
            AUTHORITY["EPSG","8901"]],
         UNIT["degree",0.0174532925199433,
            AUTHORITY["EPSG","9122"]],
         AUTHORITY["EPSG","4326"]],
      PROJECTION["Mercator_1SP"],
      PARAMETER["central_meridian",0],
      PARAMETER["scale_factor",1],
      PARAMETER["false_easting",0],
      PARAMETER["false_northing",0],
      UNIT["metre",1,
         AUTHORITY["EPSG","9001"]],
      AXIS["X",EAST],
      AXIS["Y",NORTH],
      AUTHORITY["EPSG","3857"]]
   ```

3. **Proj4 String**  
   ```
   +proj=merc +a=6378137 +b=6378137 +lat_ts=0 +lon_0=0 +x_0=0 +y_0=0 +k=1 +units=m +nadgrids=@null +wktext +no_defs
   ```

Each of these formats represents the same Web Mercator projection, allowing
different GIS tools and libraries to interpret and use the coordinate system
consistently.

In general, *WKT* is the most complete and flexible format, representing any
possible CRS; *EPSG* codes are standardized and widely recognized, but
only cover a limited number of possibilities; and *Proj4* string are concise and
easy to use in programming environments, but may not cover all possible
options. Since in `crstools` we attempt to provide a wide range of CRS options,
we use the *WKT* and *Proj4* format for the CRS suggestions.

## Spatial objects in R

In R, the `sf` (Simple Features) package is widely used for handling vector
geographic data, such as points, lines, and polygons. It provides an efficient
and standardized way to work with spatial data using the **simple features**
model, allowing seamless integration with the `tidyverse` and modern spatial
workflows. The `sf` package supports various spatial operations, including
transformations, intersections, and aggregations, while also providing
compatibility with the system libraries **GDAL, PROJ, and GEOS** for advanced spatial analysis. For
raster data, the `terra` package is the preferred choice, offering a powerful
and efficient framework for handling large raster datasets. `terra` is designed to
replace the older `raster` package, providing better performance and memory
management for processing multi-layer rasters, reprojecting, cropping, and
performing raster algebra. Together, `sf` and `terra` enable a comprehensive
geospatial analysis workflow in R, ensuring smooth interoperability between
vector and raster data.

`sf` and `terra` each have their own functions to query and set CRS: `sf::st_crs()`
and `terra::crs()`, respectively. These functions allow you to access and modify
the CRS of spatial objects, ensuring that your data is correctly projected and
aligned for spatial analysis. To avoid confusion between the two packages, in
this vignette we will use the convention of giving the package name before the
function name, such as `sf::st_crs()` and `terra::crs()`.

Note that vector data can also be represented in `terra` with the `SpatVector`
class, and raster data can be represented in `sf` with the `stars` class 
(via the `stars` package). These two options are less used, but they can each
be queried by `terra::crs()` and `sf::st_crs()` depending on the package they
belong to; `crstools` can handle these objects as well.

## Choosing a CRS

The `suggest_crs()` function provides suggestions for CRS for different types of
maps. The function takes the extent of the map and the type of distortion as
inputs and returns a CRS (or a list of CRSs if more than one option is available).
There are 3 types of distortion: *equal area*, *conformal*, and *equidistant*.
As suggested by the name, equal-area maps preserve the area of features, conformal
maps preserve angles, and equidistant maps preserve distances. A fourth option,
*compromise*, is also available, which tries to balance the trade-offs between
the other three types of distortion.

The choice of CRS depends on the extent of the map. Let us use South America
as an example. We will use `rnaturalearth` to get a map of that continent, as
an `sf` object (`sf` is generally used for ):
```{r get_s_america}
library(rnaturalearth)
library(sf)
s_america_sf <- ne_countries(continent = "South America", returnclass = "sf")
```

We use `ggplot2` to visualise it:
```{r plot_s_america}
library(ggplot2)
ggplot() +
  geom_sf(data = s_america_sf) +
  theme_minimal()
```


Let us check the CRS for this object:
```{r check_crs}
sf::st_crs(s_america_sf)
```
Note that, by default, `sf` gives us the CRS in WKT format. At the very bottom,
you can see that, as this is a standard CRS, the EPSG code is also given.
If we want to see
the Proj4 code, we can use the `sf::st_crs()` and extract the `$proj4string` slot:
```{r check_crs_epsg}
sf::st_crs(s_america_sf)$proj4string
```

Now, let's use `suggest_crs()` to get a suggestion for an equal-area CRS 
for this map:
```{r}
library(crstools)
s_am_equal_area <- suggest_crs(s_america_sf, distortion = "equal_area")
```

The function returns a list with two elements: `proj4` and `wkt`. Let's see:
```{r}
s_am_equal_area
```

We can now use this CRS to reproject the map. We can quickly inspect the projection
by plotting the map with the new CRS:
```{r plot_s_america_equidist}
ggplot() +
  geom_sf(data = s_america_sf) +
  coord_sf(crs = s_am_equal_area$proj4) +
  theme_minimal()
```

The Tissot indicatrix is a mathematical contrivance used in cartography to
characterize local distortions due to map projection. The `geom_tissot()` function
adds Tissot's indicatrix to a map, showing how circles are distorted by the
projection. Let's use it to add Tissot's indicatrix to the map of South America
with the equal-area projection:
```{r}
ggplot(data = s_america_sf) +
  geom_sf() +
  geom_tissot() +
  coord_sf(crs = s_am_equal_area$proj4) +
  theme_minimal()
```

Let us compare this to a simple Mercator projection:
```{r}
ggplot(data = s_america_sf) +
  geom_sf() +
  geom_tissot() +
  coord_sf(crs = "+proj=merc") +
  theme_minimal()
```

We can see how the area of the circles is distorted in the Mercator projection
as we move further away from the equator, while the equidistant projection
preserves the area of the circles better.

We could use the CRS to also project the `sf` object, and check that the new
CRS was indeed adopted:
```{r}
s_america_sf_equal_area <- sf::st_transform(s_america_sf, s_am_equal_area$proj4)
sf::st_crs(s_america_sf_equal_area)$proj4string
```

Now we can plot it, and the new CRS will be used automatically:
```{r plot_sa_equidist_proj}
ggplot() +
  geom_sf(data = s_america_sf_equal_area) +
  theme_minimal()
```

We can do the same with a raster object. We'll use a raster from the package
`pastclim`, which allows the easy retrieval of climate data. First, ensure
`pastclim` is installed and the data_path is set, using `set_data_path()`.

Let's get the annual mean temperature (variable "bio01") for Europe:
```{r get_raster_europe}
library(terra)
library(pastclim)

# get example dataset from cran
set_data_path(on_CRAN = TRUE)

europe_r <- region_slice(
  time_bp = 0,
  bio_variables = c("bio01"),
  dataset = "Example",
  ext = region_extent$Europe
)
```

A quick summary of the object gives the CRS, as "coord. ref.":
```{r}
europe_r
```

We can inspect in more detail the CRS for a terra object with:
```{r}
terra::crs(europe_r)
```
which gives a WKT.

We can ask for the "Proj4" string with:
```{r}
terra::crs(europe_r, proj = TRUE)
```

There is also a parameter "describe" which returns the EPSG code.

Let's plot it with `ggplot2`, using `tidyterra` to visualise the raster:

```{r plot_raster_s_america}
library(tidyterra)
ggplot() +
  geom_spatraster(data = europe_r) +
  scale_fill_viridis_c(na.value = "transparent") +
  theme_minimal()
```

We can now ask for a suggestion for an equal-area projection for this raster:
```{r}
europe_r_equal_area <- suggest_crs(europe_r, distortion = "equal_area")
europe_r_equal_area
```

Let's use it:
```{r plot_raster_europe_ea}
ggplot() +
  geom_spatraster(data = europe_r) +
  scale_fill_viridis_c(na.value = "transparent") +
  coord_sf(crs = europe_r_equal_area$proj4) +
  theme_minimal()
```

Let us use the Tissot indicatrix to assess how effective the equal area
projection is:
```{r}
ggplot() +
  geom_spatraster(data = europe_r) +
  scale_fill_viridis_c(na.value = "transparent") +
  geom_tissot(data = europe_r) +
  coord_sf(crs = europe_r_equal_area$proj4) +
  theme_minimal()
```

Comparing it to the Mercator projection:
```{r}
ggplot() +
  geom_spatraster(data = europe_r) +
  scale_fill_viridis_c(na.value = "transparent") +
  geom_tissot(data = europe_r) +
  coord_sf(crs = "+proj=merc") +
  theme_minimal()
```

We can also use the CRS to reproject the raster, and check that it has been
applied correctly:
```{r}
europe_r_equal_area <- terra::project(europe_r, europe_r_equal_area$proj4)
terra::crs(europe_r_equal_area, proj = TRUE)
```

If we now plot the raster, the CRS is automatically added:
```{r plot_raster_europe_ea_proj}
ggplot() +
  geom_spatraster(data = europe_r_equal_area) +
  scale_fill_viridis_c(na.value = "transparent") +
  theme_minimal()
```