Geographic Coordinate Systems
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.
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.
Projection in 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 fieldsproj4string(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 unprojected 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 fieldsproj4string(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 European Petroleum Survey Group (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,0The 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"