Map projections. Everyone’s favorite subject. Ok, ok, you may have detected some sarcasm. In this mini-tutorial, I am going to try to convince you that *yes*, in fact, map projections are very cool and very important. Check out this video for a reminder of the importance of projection. Let’s go!

A GCS uses a sphere to locate you on Earth’s surface.1 YES, I know, the Earth is not a sphere, but more of a potato… we’ll get to that..

2 Visualization taken from this amazing reference

Have you ever wondered why longitude and latutide are measured in degrees? This is why! Your location on the globe is described using the angle at which we can find your location relative to midway between the poles (that’s the equator y’all) and the prime meridian.3 The prime meridian was established in 1884 at the International Meridian Conference in DC, one of the coolest conferences ever held. The intersection of the equator and the prime merdian form coordinates `(0,0)`

in our GCS. Latitude ranges from -90 degrees at the South Pole to 90 degrees at the North Pole. This should make sense if you look at the figure above. Longitude runs from 180 degrees when heading east of the prime meridian and -180 degrees when traveling west. This leaves parts of Fiji and the Chukotka Autonomous Okrug with some confusing coordinates.

There are many different GCS, each defined by the way we think of Earth’s shape. Some GCS project latitude and longitude onto a sphere, others a spheroid. They are also defined by the way we think of the `(0,0)`

point in the middle of the earth. A *datum* defines the position of our abstract Earth in sphere or spheroid form relative to the center of the Earth. The following is the best visualization I’ve found of this, also taken from this amazing reference:

What’s important to understand here is that while latitude and longitude make it easy to locate exact positions on our planet’s surface, they are NOT UNIFORM. Look at the globe… as latidude moves north, the regions defined by each grid cell become smaller. This means we shouldn’t measure distances or areas using unprojected latitude and/or longitude. To do that, we need projected coordinate systems.

PCS, unline GCS, are defined on a two dimensional surface, i.e. lengths, angles, and areas are constant across these two dimensions. PCS are generally tailored to the pecularities of specific regions. This is a way of addressing Eath’s lumpiness and distortions. To translate a 3D shape onto a 2D surface, we use *map projections*. Imagine a light shining through Earth onto a 2D surface, as shown below. Yet again, the visualization is provided by this amazing reference. One of the most important things you can remember about map projections is that any time you transform a 3D object into a 2D object, you will distort shape, area, and representation of the 3D object. Map projections, therefore, are not neutral. The projection you choose often is not neutral. For example, conformal projections preserve local shape at the expense of distorting larger regions. Equal area projections preserve the area of shapes displayed at the expense of angle and scale. Equidistant projections preserve distance between points but distorts scale.

`R`

You can determine the projection of any `sp`

object by calling the `proj4string()`

function. This function reports the projection string in the `proj4string`

slot of the `sp`

object.

```
library(rgdal)
states <- readOGR("./data/states.shp")
```

```
## OGR data source with driver: ESRI Shapefile
## Source: "./data/states.shp", layer: "states"
## with 51 features
## It has 5 fields
```

`proj4string(states)`

`## [1] "+proj=longlat +datum=NAD83 +no_defs +ellps=GRS80 +towgs84=0,0,0"`

There are many parameters you can add to your `Spatial`

object’s projection. The full list can be found here. Our `states`

object contains `proj`

, `datum`

, `ellps`

, and `towgs84`

. Let’s unpack this. `proj`

refers to the projection of the object, in our case, the data is simply in a Geographic Coordinate System (i.e. using latitude and longitude to describe location), so `+proj=longlat`

. Really, we can think of this dataset as *un*projected since it’s being displayed using a GCS. Our datum is `NAD83`

, or the North American Datum of 1983, a commonly used reference for the US. This is indicated in the `proj4string`

by `+datum=NAD83`

. We have no additional projection info since we’re working with a GCS (`+no_defs`

). We’re using the Geodetic Reference System 1980 as our model of Earth’s shape and indicating this by `+ellps=GRS80`

. Finally, the `+towgs84`

is an additional that allows to transform our datum to one of the most commonly used datums, WGS84. In our case, we’re not making the transform, so we have `+towgs84=0,0,0`

. In this class, you will be mostly working with data that already stores projection information, so you just need to know how to read the `proj4string()`

output to identify your projection. At other times, you may read in spatial data and realize it has no projection information (`proj4string(shp) = NA`

). If you know the projection information, you can write your own `proj4string`

and assign it using `shp <- CRS("+proj=longlat, etc, etc")`

.

Let’s break apart another `proj4string`

output for practice. First, let’s look at a shapefile of conservation easements4 What’s a conservation easement? *“In a conservation easement, a landowner voluntarily agrees to sell or donate certain rights associated with his or her property - often the right to subdivide or develop - and a private organization or public agency agrees to hold the right to enforce the landowner’s promise not to exercise those rights. In essence, the rights are forfeited and no longer exist.”* More info can be found here. in Utah.

`CE <- readOGR("./data/ConservationEasements.shp")`

```
## OGR data source with driver: ESRI Shapefile
## Source: "./data/ConservationEasements.shp", layer: "ConservationEasements"
## with 159 features
## It has 15 fields
```

`proj4string(CE)`

`## [1] "+proj=utm +zone=12 +datum=NAD83 +units=m +no_defs +ellps=GRS80 +towgs84=0,0,0"`

Much of the GIS data managed by the state government of Utah will come in this flavor of projection. This data is projected using the Universal Transverse Mercator projection, in Zone 12, with the 1983 North American Datum. Nice! Notice, importantly, that since this is a *projected* CRS (not a geographical CRS in latitude and longitude), since we have projected our data onto a 2D surface. We can see an additional argument in our `proj4string`

that tells us the units of our projection, `units=m`

or meters!

Another common `proj4string`

notation involves a code called the `epsg`

. This comes from the **E**uropean **P**etroleum **S**urvey **G**roup (now Oil and Gas Producers), who, starting in the 1980s, began collecting and standardizing geodetic parameters to help normalize projection. Yes, big gas and oil even influences projections. You can use the `EPSG`

code to build your `proj4string`

. The full list of `EPSG`

entries can be found here. www.spatialreference.org is another great reference. If we know, for example, that the `EPSG`

for `WGS84 UTM Zone 12N`

, the UTM zone for Utah, is `32612`

(just Google `WGS84 ZOne 12N EPSG`

and you should find the `EPSG`

), we can build the `proj4string`

as follows:

`CRS("+init=epsg:32612")`

```
## CRS arguments:
## +init=epsg:32612 +proj=utm +zone=12 +datum=WGS84 +units=m +no_defs
## +ellps=WGS84 +towgs84=0,0,0
```

The `CRS()`

function helps us build `proj4string`

strings. To learn more, read through `?CRS()`

.

Finally, what if we have data in one projection and we want to *reproject* that data into another projection? We’ll cover the basics with `spTransform()`

here.5 For advanced projection questions, check out *Applied Spatial Data Analysis with R* by Bivand et al. Pages 84-88 are paricularly helpful. Imagine you want to transform your conservation easement shapefile into the UTM Zone 12N projection. How would we do this in

`R`

?```
CEp <- spTransform(CE, CRS("+init=epsg:32616"))
proj4string(CEp)
```

`## [1] "+init=epsg:32616 +proj=utm +zone=16 +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0"`

`spTransform()`

takes first an argument pointing to the object to be transformed, in our case `CE`

, and an object of class `CRS`

that describes the projection. Remember, `CRS`

refers to coordinate reference system. In `R`

, to describe a `Spatial`

object’s projection, we must format the projection information into a `CRS`

object. We can build a `CRS`

object by feeding the `CRS()`

function the `EPSG`

information for our projection of interest. If we happen to know the entire `proj4string`

string, the we can also feed the `CRS()`

function the entire string, i.e.:

```
CEp <- spTransform(CE, CRS("+init=epsg:32616 +proj=utm +zone=16 +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0"))
proj4string(CEp)
```

`## [1] "+init=epsg:32616 +proj=utm +zone=16 +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0"`