In arc2r/book: All things R Spatial or A guide to spatial R for Recovering ArcGIS Users

Projections and Transformations Toolset {#projections}

For this chapter, you will need the following R Packages:

library(arc2r)
library(sf)
library(raster)

If there is anything special about spatial, it's the Coordinate Reference System. We will not get into the topic of CRS and Projections here, but it is advisable to get up to speed on this topic if you are dealing with geodata. From this toolset, we will cover the following tools:

• Define Projection
• Project
• Project Raster
• Rescale (Resizes a raster by the specified x and y scale factors)

Define Projection {#projections-defineprojection}

Overwrites the coordinate system information (map projection and datum) stored with a dataset. This tool is intended for datasets that have an unknown or incorrect coordinate system defined.

If your dataset comes without an assigned CRS, things can get hairy. Usually, you can make an informed guess as to which CRS the dataset is associated with, but in the end only the data producer can tell you. Let's take the dataset bezirke from the package arc2r. As you can see, the CRS field in the data's header is NA.

data("bezirke")

Calling the function st_crs() on bezirke confirms this.

st_crs(bezirke)

In ArcGIS, we would use the Tool "Define projection". In R, we use the function st_crs(dataset) <- value, where value is the CRS of our dataset. It can be either

1. A string accepted by GDAL (type character)
2. A valid EPSG value (type integer) or
3. An object of class crs

We find option two to be the simplest approach. Since we know the dataset bezirke is in the new swiss coordinate system (CH1903+ / LV95), we can find out the EPSG code (e.g. by consulting epsg.io). We will not get into options 1 and 3 and only

st_crs(bezirke) <- 2056

If we look at our header again, we see that we have a new field projected CRS which is in fact displaying the correct name based on the EPSG code we provided.

bezirke

```{block2, type = "rmdnote", echo = TRUE, purl = FALSE} <!--

TODO revise this

--> st_set_crs function does not reproject the coordinates of the given dataset. In other words, it does not affect the actual geometry column of the sf object. st_tranform on the other hand indeed does indeed reproject the dataset to another coordinate

### Project {#projections-project}

> Projects spatial data from one coordinate system to another

Now that `bezirke` has an assigned `CRS` (see \@ref(projections-defineprojection), we can transform it into a new coordinate system. To emphasize this, will will visualize the dataset using `base::plot()` with the option `axes = TRUE`, to visualize the current coordinate system.  

```r
plot(bezirke["area_km2"], axes = TRUE)

bezirke_wgs84 <- st_transform(bezirke, 4326)

plot(bezirke_wgs84["area_km2"], axes = TRUE)

st_crs(bezirke)

bezirke_swiss <- st_transform(bezirke, 2056)

# retrieve the coordinate system
st_crs(bezirke_swiss)

Project Raster {#projections-projectraster}

Transforms a raster dataset from one coordinate system to another

Working with Raster datasets in GIS of operations is of equal importance, as working with vector ones. One of the spatial properties of raster datasets is the the Coordinate Reference System (CRS). CRS is the specific system that “associates” the raster coordinates (which are just pairs of x/y values) to geographic locations. In ArcGIS pro the tool for projecting a raster dataset is called Project Raster (Data Management). Let's see how we can perform the same operation with R.

# Dataset derived from the spatial interpolation of all the available "recycling points"
# in the city of Wädenwil
data("recycling_raster")

# Dataset representing the public transport quality in the city of Wädenswil
data("public_transport_waedi")# CRS -> WGS84 

# Plot the raster dataset - World Geodetic System 1984 
plot(public_transport_waedi,las=1,
     main = "Quality of public transport in the city of Wädenwil - CRS: WGS84",
     cex.main=1,font.main=4)

We can use the projectRaster() function to reproject a raster into a new CRS. The first argument of the aforementioned function is the raster dataset we want to reproject, while the second one is the dataset to whose projection we are targeting to. So, in our case, we are targeting to the coordinate system of the raster_recycling dataset.
It is important to remember that raster reprojection only works when the raster object has already a defined CRS.

# Transform the coordinate system of the raster dataset publicTransport_CH
# into the Swiss Coordinate system - CH1903+LV95
publicTransport_CH = projectRaster(public_transport_waedi, recycling_raster)

# Plot the raster dataset - Swiss Coordinate System CH1903+LV95 
plot(publicTransport_CH,las=1,
     main = "Quality of public transport in the city of Wädenwil - CRS: CH1903+LV95",
     cex.main=1,font.main=4)

arc2r/book documentation built on March 5, 2021, 2:10 p.m.