A07-Coordinate-Systems.Rmd

# Coordinate Systems in R 

```{r, include=FALSE}
# Code chunk used to add scrollbar to lengthy WKT output
options(width = 60)
```

```{r, echo = FALSE}
source("package_list.R")
get.pckg.info("A07-Coordinate-Systems.Rmd") 
```


## A note about the changes to the `PROJ` environment {-}

Newer versions of `sf` make use of the `PROJ 6.0` C library or greater. Note that the version of `PROJ` is not to be confused with the version of the `proj4` R package--the `proj4` and `sf` packages make use of the `PROJ` C library that is developed independent of R. You can learn more about the `PROJ` development at [proj.org](https://proj.org/). 

There has been a significant change in the `PROJ` library since the introduction of version `6.0`. This has had serious implications in the development of the R spatial ecosystem. As such, if you are using an older version of `sf` or `proj4` that was developed with a version of `PROJ` older than 6.0, some of the input/output presented in this appendix may differ from yours.

## Sample files for this exercise{-}

Data used in this exercise can be loaded into your current R session by running the following chunk of code.

```{r}
library(terra)

z <- gzcon(url("https://github.com/mgimond/Spatial/raw/main/Data/elev.RDS"))
elev.r <- unwrap(readRDS(z))
z <- gzcon(url("https://github.com/mgimond/Spatial/raw/main/Data/s_sf.RDS"))
s.sf <- readRDS(z)
```

We'll make use of two data layers in this exercise: a Maine counties polygon layer (`s.sf`) and an elevation raster layer (`elev.r`). The former is in an `sf` format and the latter is in a `SpatRaster` format. 

## Loading the `sf` package {-}

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

Note the versions of `GEOS`, `GDAL` and `PROJ` the package `sf` is linked to. Different versions of these libraries may result in different outcomes than those shown in this appendix. You can check the linked library versions as follows:

```{r}
sf_extSoftVersion()[1:3]
```


## Checking for a coordinate system {-}

To extract coordinate system (CS) information from an `sf` object use the `st_crs` function.


```{r attr.output = 'style="max-height: 100px;direction: ltr;"'}
st_crs(s.sf)
```



With the newer version of the `PROJ` C library, the coordinate system is defined using the *Well Known Text* (**WTK**/**WTK2**) format which consists of a series of `[...]` tags. The WKT format will usually start with a `PROJCRS[...]` tag for a projected coordinate system, or a `GEOGCRS[...]` tag for a geographic coordinate system.

The CRS output will also consist of a user defined CS definition which can be an **EPSG** code (as is the case in this example), or a string defining the datum and projection type.

You can also extract CS information from a `SpatRaster` object use the `st_crs` function. 

```{r attr.output = 'style="max-height: 100px;direction: ltr;"'}
st_crs(elev.r)
```


Up until recently, there has been two ways of defining a coordinate system: via the **EPSG** numeric [code](http://spatialreference.org/ref/epsg/) or via the **PROJ4** [formatted string](https://proj4.org/apps/proj.html). Both can be used with the `sf` and `SpatRast` objects.

With the newer version of the `PROJ` C library, you can also define an `sf` object's  coordinate system using the *Well Known Text* (**WTK**/**WTK2**) format. This format has a more elaborate syntax (as can be seen in the  previous outputs) and may not necessarily be the easiest way to manually define a CS. When possible, adopt an **EPSG** code which comes from a well established authority. However, if customizing a CS, it may be easiest to adopt a **PROJ4** syntax. 

## Understanding the Proj4 coordinate syntax {-}

The [PROJ4](http://proj4.org/) syntax consists of a list of parameters, each prefixed with the `+` character. For example, `elev.r`'s CS is in a UTM projection (`+proj=utm`) for zone 19 (`+zone=19`) and in an NAD 1983 datum (`+datum=NAD83`).

A list of a few of the PROJ4 parameters used in defining a coordinate system follows. Click [here](https://proj4.org/usage/projections.html) for a full list of parameters.

```
+a         Semimajor radius of the ellipsoid axis
+b         Semiminor radius of the ellipsoid axis
+datum     Datum name 
+ellps     Ellipsoid name 
+lat_0     Latitude of origin
+lat_1     Latitude of first standard parallel
+lat_2     Latitude of second standard parallel
+lat_ts    Latitude of true scale
+lon_0     Central meridian
+over      Allow longitude output outside -180 to 180 range, disables wrapping 
+proj      Projection name 
+south     Denotes southern hemisphere UTM zone
+units     meters, US survey feet, etc.
+x_0       False easting
+y_0       False northing
+zone      UTM zone
```

You can view the list of available projections `+proj=` [here](https://proj4.org/operations/projections/).

## Assigning a coordinate system {-}

A coordinate system definition can be passed to a spatial object. It can either fill a spatial object's empty CS definition or it can overwrite its existing CS definition (the latter should only be executed if there is good reason to believe that the original definition is erroneous). Note that this step **does not** change an object's underlying coordinate *values* (this process will be discussed in the next section).

We'll pretend that a CS definition was not assigned to `s.sf` and assign one manually using the `st_set_crs()` function. In the following example, we will define the CS using the proj4 syntax.

```{r}
s.sf <- st_set_crs(s.sf, "+proj=utm +zone=19 +ellps=GRS80 +datum=NAD83") 
```

Let's now check the object's CS. 

```{r attr.output = 'style="max-height: 100px;direction: ltr;"'}
st_crs(s.sf)
```

You'll note that the `User input:` field now shows the `proj4` string as defined in our call to the `st_set_crs()` function. But you'll also note that some of the parameters in the WKT string such as the `PROJCRS[...]` and `BASEGEOGCRS[...]` tags are not defined (`unknown`). This is not necessarily a problem given that key datum and projection information are present in that WKT string (make sure to scroll down in the output box to see the other WKT parameters). Nonetheless, it's not a bad idea to define the CS using EPSG code when one is available. We'll do this next.

The UTM NAD83 Zone 19N EPSG code equivalent is `26919`.

```{r}
s.sf <- st_set_crs(s.sf, 26919)
```

Let's now check the object's CS.

```{r attr.output = 'style="max-height: 100px;direction: ltr;"'}
st_crs(s.sf)  
```

Key projection parameters remain the same. But additional information is added to the WKT header.

You can use the PROJ4 string defined earlier for `s.sf` to define a raster's CRS using the `crs()` function as follows (here too we'll assume that the spatial object had a missing reference system or an incorrectly defined one).

```{r}
crs(elev.r) <- "+proj=utm +zone=19 +ellps=GRS80 +datum=NAD83"
```

Note that we do not need to define all of the parameters so long as we know that the default values for these undefined parameters are correct. Also note that we do not need to designate a hemisphere since the NAD83 datum applies only to North America.

Let's check the raster's CS:

```{r attr.output = 'style="max-height: 100px;direction: ltr;"'}
st_crs(elev.r)
```

To define a raster's CS using an EPSG code, use the following PROJ4 syntax:

```{r}
crs(elev.r) <- "+init=EPSG:26919"
```

```{r attr.output = 'style="max-height: 100px;direction: ltr;"'}
st_crs(elev.r) 
```

To recreate a CS defined in a software such as ArcGIS, it is best to extract the CS' **WKID/EPSG** code, then use that number to look up the PROJ4 syntax on [http://spatialreference.org/ref/](http://spatialreference.org/ref/). For example, in ArcGIS, the WKID number can be extracted from the coordinate system properties output. 

```{r fA2-arcgis, echo=FALSE, fig.cap = "An ArcGIS dataframe coordinate system properties window. Note the WKID/EPSG code of 26919 (highlighted in red) associated with the NAD 1983 UTM Zone 19 N CS.", fig.align='center', out.width=471}

knitr::include_graphics("img/ArcGIS_UTMNAD83.jpg")
```

That number can then be entered in the [http://spatialreference.org/ref/](http://spatialreference.org/ref/)'s search box to pull the Proj4 parameters (note that you must select **Proj4** from the list of syntax options).

```{r fA2-epsg, echo=FALSE, fig.cap = "Example of a search result for EPSG **26919** at [http://spatialreference.org/ref/](http://spatialreference.org/ref/). Note that after clicking the `EPSG:269191` link, you must then select the Proj4 syntax from a list of available syntaxes to view the projection parameters", fig.align='center', out.width=420}

knitr::include_graphics("img/EPSG_search.jpg")
```

Here are examples of a few common projections:

Projection                          WKID    Authority    Syntax
---------------------------------  ------- -----------   --------------------
UTM NAD 83 Zone 19N                 26919    EPSG        `+proj=utm +zone=19 +ellps=GRS80 +datum=NAD83 +units=m +no_defs`
USA Contiguous albers equal area    102003   ESRI        `+proj=aea +lat_1=29.5 +lat_2=45.5 +lat_0=37.5 +lon_0=-96 +x_0=0 +y_0=0 +ellps=GRS80 +datum=NAD83 +units=m +no_defs`
Alaska albers equal area            3338     EPSG        `+proj=aea +lat_1=55 +lat_2=65 +lat_0=50 +lon_0=-154 +x_0=0 +y_0=0 +ellps=GRS80 +datum=NAD83 +units=m +no_defs` 
World Robinson                      54030    ESRI        `+proj=robin +lon_0=0 +x_0=0 +y_0=0 +ellps=WGS84 +datum=WGS84 +units=m +no_defs`


## Transforming coordinate systems {-}

The last section showed you how to _define_ or _modify_ the coordinate system _definition_. This section shows you how to **transform** the coordinate values associated with the spatial object to a *different* coordinate system. This process calculates new coordinate values for the points or vertices defining the spatial object.

For example, to transform  the `s.sf`  vector object to a WGS 1984 geographic (long/lat) coordinate system, we'll use the `st_transform` function.

```{r attr.output = 'style="max-height: 100px;direction: ltr;"'}
s.sf.gcs <- st_transform(s.sf, "+proj=longlat +datum=WGS84")
st_crs(s.sf.gcs)
```

Using the EPSG code equivalent (`4326`) instead of the proj4 string yields:

```{r attr.output = 'style="max-height: 100px;direction: ltr;"'}
s.sf.gcs <- st_transform(s.sf, 4326)
st_crs(s.sf.gcs)
```

This approach may add a few more tags (These reflect changes in datum definitions in newer versions of the PROJ library) but, the coordinate values should be the same

To transform a raster object, use the `project()` function. 

```{r attr.output = 'style="max-height: 100px;direction: ltr;"'}
elev.r.gcs <- project(elev.r, y="+proj=longlat +datum=WGS84")
st_crs(elev.r.gcs)
```

If an EPSG code is to be used, adopt the `"+init=EPSG: ..."` syntax used earlier in this tutorial.

```{r}
elev.r.gcs <- project(elev.r, y="+init=EPSG:4326")
st_crs(elev.r.gcs)
```

A geographic coordinate system is often desired when overlapping a layer with a web based mapping service such as Google, Bing or OpenStreetMap (even though these web based services end up projecting to a projected coordinate system--most likely a Web Mercator projection). To check that `s.sf.gcs` was properly transformed, we'll overlay it on top of an OpenStreetMap using the `leaflet` package.

```{r fig.height=3}
library(leaflet)
leaflet(s.sf.gcs) %>% 
  addPolygons() %>% 
  addTiles()
```

<p>

Next, we'll explore other transformations using a `tmap` dataset of the world

```{r attr.output = 'style="max-height: 100px;direction: ltr;"'}
library(tmap)
data(World)  # The dataset is stored as an sf object
 
# Let's check its current coordinate system
st_crs(World)
```


The following chunk transforms the world map to a custom azimuthal equidistant projection centered on latitude `0` and longitude `0`. Here, we'll use the proj4 syntax.

```{r white-space: pre-wrap;  , }
World.ae <- st_transform(World, "+proj=aeqd +lat_0=0 +lon_0=0 +x_0=0 +y_0=0 +ellps=WGS84 +datum=WGS84 +units=m +no_defs")
```

Let's check the CRS of the newly created vector layer

```{r attr.output = 'style="max-height: 100px;direction: ltr;"'}
st_crs(World.ae) 
```

Here's the mapped output:

```{r fig.height=2.5, fig.width=2.5}
tm_shape(World.ae) + tm_fill() 
```

The following chunk transforms the world map to an Azimuthal equidistant projection centered on Maine, USA (69.8&deg; West, 44.5&deg; North) .
 
```{r fig.height=2.5, fig.width=2.5}
World.aemaine <- st_transform(World, "+proj=aeqd +lat_0=44.5 +lon_0=-69.8 +x_0=0 +y_0=0 +ellps=WGS84 +datum=WGS84 +units=m +no_defs")

tm_shape(World.aemaine) + tm_fill()  
```

The following chunk transforms the world map to a World Robinson projection.
 
```{r fig.height=2.5, fig.width=4.5}
World.robin <- st_transform(World,"+proj=robin +lon_0=0 +x_0=0 +y_0=0 +ellps=WGS84 +datum=WGS84 +units=m +no_defs")
tm_shape(World.robin) + tm_fill()  
```

The following chunk transforms the world map to a World sinusoidal projection.

```{r fig.height=2.5, fig.width=4.5}
World.sin <- st_transform(World,"+proj=sinu +lon_0=0 +x_0=0 +y_0=0 +ellps=WGS84 +datum=WGS84 +units=m +no_defs")
tm_shape(World.sin) + tm_fill()  
```

The following chunk transforms the world map to a World Mercator projection.

```{r fig.height=2.5, fig.width=2.5}
World.mercator <- st_transform(World,"+proj=merc +lon_0=0 +k=1 +x_0=0 +y_0=0 +ellps=WGS84 +datum=WGS84 +units=m +no_defs")
tm_shape(World.mercator) + tm_fill()  
```


### Reprojecting to a new meridian center {-}

An issue that can come up when transforming spatial data is when the location of the tangent line(s) or points in the CS definition forces polygon features to be split across the 180&deg; meridian. For example, re-centering the Mercator projection to -69&deg; will create the following output.

```{r fig.height=2.5, fig.width=3}
World.mercator2 <- st_transform(World, "+proj=merc +lon_0=-69 +k=1 +x_0=0 +y_0=0 +ellps=WGS84 +datum=WGS84 +units=m +no_defs")

tm_shape(World.mercator2) + tm_borders()
```


The polygons are split and R does not know how to piece them together.

One solution is to split the polygons at the new meridian using the `st_break_antimeridian` function before projecting to a new re-centered coordinate system.

```{r fig.height=2.5, fig.width=3}
# Define new meridian
meridian2 <- -69

# Split world at new meridian
wld.new <- st_break_antimeridian(World, lon_0 = meridian2)

# Now reproject to Mercator using new meridian center
wld.merc2 <- st_transform(wld.new, 
                          paste("+proj=merc +lon_0=", meridian2 ,
                                "+k=1 +x_0=0 +y_0=0 +ellps=WGS84 +datum=WGS84 +units=m +no_defs") )

tm_shape(wld.merc2) + tm_borders()

```

This technique can be applied to any other projections. Here's an example of a Robinson projection.

```{r fig.height = 2.5, fig.width=4}
wld.rob.sf <-  st_transform(wld.new, 
                            paste("+proj=robin +lon_0=", meridian2 ,
                                  "+k=1 +x_0=0 +y_0=0 +ellps=WGS84 +datum=WGS84 +units=m +no_defs") )
tm_shape(wld.rob.sf) + tm_borders()
```
 
## A note about containment {-}

While in theory, a point completely enclosed by a bounded area should always remain bounded by that area in any projection, this is not always the case in practice. This is because the transformation applies to the vertices that define the line segments and not the lines themselves. So if a point is inside of a polygon and very close to one of its boundaries in its native projection, it may find itself on the other side of that line segment in another projection hence outside of that polygon. In the following example, a polygon layer and point layer are created in a Miller coordinate system where the points are enclosed in the polygons.

```{r fig.height=2.5,fig.width=4}
# Define a few projections
miller <- "+proj=mill +lat_0=0 +lon_0=0 +x_0=0 +y_0=0 +ellps=WGS84 +datum=WGS84 +units=m +no_defs"
lambert <- "+proj=lcc +lat_1=20 +lat_2=60 +lat_0=40 +lon_0=-96 +x_0=0 +y_0=0 +ellps=GRS80 +datum=NAD83 +units=m +no_defs"

# Subset the World data layer and reproject to Miller
wld.mil <- subset(World, iso_a3 == "CAN" | iso_a3 == "USA") |>
  st_transform(miller)

# Create polygon and point layers in the Miller projection  
sf1 <- st_sfc( st_polygon(list(cbind(c(-13340256,-13340256,-6661069, -6661069, -13340256),
                               c(7713751, 5326023, 5326023,7713751, 7713751 )))), crs = miller) 

pt1 <- st_sfc( st_multipoint(rbind(c(-11688500,7633570), c(-11688500,5375780),
                                  c(-10018800,7633570), c(-10018800,5375780),
                                  c(-8348960,7633570), c(-8348960,5375780))),  crs = miller)
pt1 <- st_cast(pt1, "POINT") # Create single part points

# Plot the data layers in their native projection
tm_shape(wld.mil) +tm_fill(col="grey") + 
  tm_graticules(x = c(-60,-80,-100, -120, -140), 
                y = c(30,45, 60), 
                labels.col = "white", col="grey90") +
  tm_shape(sf1) + tm_polygons("red", alpha = 0.5, border.col = "yellow") +
  tm_shape(pt1) + tm_dots(size=0.2) 

```

The points are close to the boundaries, but they are inside of the polygon nonetheless. To confirm, we can run `st_contains` on the dataset:

```{r}
st_contains(sf1, pt1)
```

All six points are selected, as expected.

Now, let's reproject the data into a Lambert conformal projection.

```{r, fig.height=2.5, fig.width=4}
# Transform the data
wld.lam <- st_transform(wld.mil, lambert)
pt1.lam <- st_transform(pt1, lambert)
sf1.lam <- st_transform(sf1, lambert)

# Plot the data in the Lambert coordinate system
tm_shape(wld.lam) +tm_fill(col="grey") + 
  tm_graticules( x = c(-60,-80,-100, -120, -140), 
                 y = c(30,45, 60), 
                 labels.col = "white", col="grey90") +
  tm_shape(sf1.lam) + tm_polygons("red", alpha = 0.5, border.col = "yellow") +
  tm_shape(pt1.lam) + tm_dots(size=0.2)   

```

Only three of the points are contained. We can confirm this using the `st_contains` function:

```{r}
st_contains(sf1.lam, pt1.lam)
```

To resolve this problem, one should _densify_ the polygon by adding more vertices along the line segment. The vertices density will be dictated by the resolution needed to preserve the map's containment properties and is best determined experimentally.

We'll use the `st_segmentize` function to create vertices at 1 km (1000 m) intervals.

```{r fig.height=2.5, fig.width=4}
# Add vertices every 1000 meters along the polygon's line segments
sf2 <- st_segmentize(sf1, 1000)

# Transform the newly densified polygon layer
sf2.lam <- st_transform(sf2, lambert)

# Plot the data
tm_shape(wld.lam) + tm_fill(col="grey") + 
  tm_graticules( x = c(-60,-80,-100, -120, -140), 
                 y = c(30,45, 60), 
                 labels.col = "white", col="grey90") +
  tm_shape(sf2.lam) + tm_polygons("red", alpha = 0.5, border.col = "yellow") +
  tm_shape(pt1.lam) + tm_dots(size=0.2) 
```

Now all points remain contained by the polygon. We can check via:

```{r}
st_contains(sf2.lam, pt1.lam)
```

## Creating Tissot indicatrix circles {-}

Most projections will distort some aspect of a spatial property, especially area and shape. A nice way to visualize the distortion afforded by a projection is to create geodesic circles.

First, create a point layer that will define the circle centers in a lat/long coordinate system.

```{r}
tissot.pt <- st_sfc( st_multipoint(rbind(c(-60,30), c(-60,45), c(-60,60),
                                  c(-80,30), c(-80,45), c(-80,60),
                                  c(-100,30), c(-100,45), c(-100,60),
                                  c(-120,30), c(-120,45), c(-120,60) )),  crs = "+proj=longlat")
tissot.pt <- st_cast(tissot.pt, "POINT") # Create single part points
```

Next we'll construct geodesic circles from these points using the `geosphere` package.

```{r}
library(geosphere)

cr.pt <- list() # Create an empty list

# Loop through each point in tissot.pt and generate 360 vertices at 300 km
# from each point in all directions at 1 degree increment. These vertices
# will be used to approximate the Tissot circles
for (i in 1:length(tissot.pt)){
  cr.pt[[i]] <- list( destPoint( as(tissot.pt[i], "Spatial"), b=seq(0,360,1), d=300000) )
}

# Create a closed polygon from the previously generated vertices
tissot.sfc <- st_cast( st_sfc(st_multipolygon(cr.pt ),crs = "+proj=longlat"), "POLYGON" )
```

We'll check that these are indeed geodesic circles by computing the geodesic area of each polygon. We'll use the `st_area` function from `sf` which will revert to geodesic area calculation if a lat/long coordinate system is present.

```{r}
tissot.sf <- st_sf( geoArea =  st_area(tissot.sfc), tissot.sfc )
```

The true area of the circles should be $\pi * r^2$ or `r round(pi * 300000^2)` square meters in our example. Let's compute the error in the `tissot` output. The values will be reported as fractions.

```{r}
( (pi * 300000^2) -  as.vector(tissot.sf$geoArea) ) / (pi * 300000^2)
```

In all cases, the error is less than 0.1%. The error is primarily due to the discretization of the circle parameter. 

Let's now take a look at the distortions associated with a few popular coordinate systems. We'll start by exploring the Mercator projection.

```{r fig.height=3, fig.width=4}
# Transform geodesic circles and compute area error as a percentage
tissot.merc <- st_transform(tissot.sf, "+proj=merc +ellps=WGS84")
tissot.merc$area_err <- round((st_area(tissot.merc, tissot.merc$geoArea)) / 
                                tissot.merc$geoArea * 100 , 2)

# Plot the map
tm_shape(World, bbox = st_bbox(tissot.merc), projection = st_crs(tissot.merc)) + 
  tm_borders() + 
  tm_shape(tissot.merc) + tm_polygons(col="grey", border.col = "red", alpha = 0.3) + 
  tm_graticules(x = c(-60,-80,-100, -120, -140), 
                y = c(30,45, 60),
               labels.col = "white", col="grey80") +
  tm_text("area_err", size=.8, alpha=0.8, col="blue")

```

The mercator projection does a good job at preserving shape, but the area's distortion increases dramatically poleward.

Next, we'll explore the Lambert azimuthal equal area projection centered at 45 degrees north and 100 degrees west.

```{r fig.height=3, fig.width=4}
# Transform geodesic circles and compute area error as a percentage
tissot.laea <- st_transform(tissot.sf, "+proj=laea +lat_0=45 +lon_0=-100 +ellps=WGS84")
tissot.laea$area_err <- round( (st_area(tissot.laea ) - tissot.laea$geoArea) / 
                                 tissot.laea$geoArea * 100, 2)

# Plot the map
tm_shape(World, bbox = st_bbox(tissot.laea), projection = st_crs(tissot.laea)) + 
  tm_borders() + 
  tm_shape(tissot.laea) + tm_polygons(col="grey", border.col = "red", alpha = 0.3) + 
  tm_graticules(x=c(-60,-80,-100, -120, -140), 
           y = c(30,45, 60),
           labels.col = "white", col="grey80") +
  tm_text("area_err", size=.8, alpha=0.8, col="blue")
```

The area error across the 48 states is near 0. But note that the shape does become slightly distorted as we move away from the center of projection.

Next, we'll explore the Robinson projection.

```{r fig.height=2.5, fig.width=4}
# Transform geodesic circles and compute area error as a percentage
tissot.robin <- st_transform(tissot.sf, "+proj=robin  +ellps=WGS84")
tissot.robin$area_err <- round(  (st_area(tissot.robin ) - tissot.robin$geoArea) / 
                                   tissot.robin$geoArea * 100, 2)

# Plot the map
tm_shape(World, bbox = st_bbox(tissot.robin), projection = st_crs(tissot.robin)) + 
  tm_borders() + 
  tm_shape(tissot.robin) + tm_polygons(col="grey", border.col = "red", alpha = 0.3) + 
  tm_graticules(x=c(-60,-80,-100, -120, -140), 
                y = c(30,45, 60),
                labels.col = "white", col="grey80") +
  tm_text("area_err", size=.8, alpha=0.8, col="blue")
```

Both shape and area are measurably distorted for the north american continent.