Overview of Coordinate Reference Systems (CRS) in R

Coordinate reference systems CRS provide a standardized way of describing locations. Many different CRS are used to describe geographic data. The CRS that is chosen depends on when the data was collected, the geographic extent of the data, the purpose of the data, etc.
In R, when data with different CRS are combined it is important to transform them to a common CRS so they align with one another. This is similar to making sure that units are the same when measuring volume or distances. EPSG codes for commonly used CRS (in the U.S.) Latitude/Longitude WGS84 (EPSG: 4326) ## Commonly used by organizations that provide GIS data for the entire globe or many countries. CRS used by Google Earth NAD83 (EPSG:4269) ##Most commonly used by U.S. federal agencies. NAD27 (EPSG: 4267) ##Old version of NAD83 Projected (Easting/Northing) UTM, Zone 10 (EPSG: 32610) ## Zone 10 is used in the Pacific Northwest Mercator (EPSG: 3857) ## Tiles from Google Maps, Open Street Maps, Stamen Maps CRS in R for sp classes: Some spatial data files have associated projection data, such as ESRI shapefiles. When readOGR is used to import these data this information is automatically linked to the R spatial object. To retrieve the CRS for a spatial object: proj4string(x) To assign a known CRS to spatial data: proj4string(x) <- CRS("+init=epsg:28992") or, type in the PROJ.4 attributes: proj4string(x) <-CRS("+proj=utm +zone=10 +datum=WGS84") To transform from one CRS to another: newData <- spTransform(x, CRS("+init=epsg:4238")) or, reference the CRS of another spatial object: newData <- spTransform(x, proj4string(OtherData)) for rasters: To retrieve or describe the CRS: projection(x) projection(x) <- CRS(+init=epsg:28992) To transform from one CRS to another: newData <- projectRaster(x, proj4string(OtherData)) Finding CRS data for your data - metadata - the .prj file of shape files (automatically detected by R when data is loaded with readOGR) - data source (Google Earth = WGS84long/lat) - can be challenging!

Package sp and rgdal is used to assign and transform CRS in R: library(rgdal) library(sp) In R, the notation used to describe the CRS is proj4string from the PROJ.4 library. It looks like this: +init=epsg:4121 +proj=longlat +ellps=GRS80 +datum=GGRS87 +no_defs +towgs84=-199.87,74.79,246.62
There are various attributes of the CRS, such as the projection, datum, and ellipsoid. Some of the options for each variable can be obtained in R with projInfo: 1. Projection: projInfo(type = "proj") 2. Datum: projInfo(type = "datum") 3. Ellipsoid: projInfo(type = "ellps")

EPSG codes A particular CRS can be referenced by its EPSG code (i.e., epsg:4121). The EPSG is a structured dataset of CRS and Coordinate Transformations. It was originally compiled by the, now defunct, European Petroleum Survey Group. Here are some websites: http://www.epsg-registry.org/ http://spatialreference.org/ (although I find these kind of confusing). In R, the details of a particular EPSG code can be obtained: CRS("+init=epsg:4326") Which returns: +init=epsg:4326 +proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs +towgs84=0,0,0 A data frame of all the CRSs can be created: EPSG <- make_EPSG() which returns a data frame with columns: code: EPSG code note: notes as included in the file prj4character: PROJ.4 attributes for the projection The EPSG dataframe can be searched: EPSG[grep("4269", EPSG$code),] EPSG[grep("HARN", EPSG$note),] tmp <- EPSG[(grep("longlat", EPSG$prj4)), ]

Page 1, Melanie Frazier

ellipses and datums

An overview of ellipses, datums, and projections
The Ellipse: Describes the generalized shape of the Earth. All mapping and coordinate systems begin with this description. There are lots of ways to do each step, resulting in lots of coordinate reference systems.

Datums
Provides the information needed to anchor the abstract coordinates to the Earth. The datum defines an origin point of the coordinate axes and defines the direction of the axes. I think of this as the base information needed to draw the imaginary coordinate lines on a globe or map. The datum always specifies the ellipsoid that is used, but the ellipsoid does not specify the datum! Datums are based on specific ellipsoids and sometimes have the same name as the ellipsoid. Available options in R projInfo(type = "datum"): name ellipse definition WGS84 WGS84 towgs84=0,0,0 GGRS87 GRS80 towgs84=199.87,74.79,246.62 NAD83 GRS80 towgs84=0,0,0 NAD27 clrk66 nadgrids=@conus, @alaska, @ntv2_0.gsb, @ntv1_can.dat potsdam bessel towgs84=606.0,23.0,413.0 carthage clark80 towgs84=-263.0,6.0,431.0 hermannskogel bessel towgs84=653.0,-212.0,449.0 ire65 mod_airy towgs84=482.530, -130.596, 564.557, -1.042, -0.214, -0.631, 8.15 nzgd49 intl towgs84=59.47, -5.04, 187.44, 0.47, -0.1, 1.024, -4.5993 OSGB36 airy towgs84=446.448, 125.157, 542.060, 0.1502, 0.2470, 0.8421, -20.4894

The Datum: Defines origin and orientation of the coordinate axes (as well the size/shape of Earth) A Map A 2D representation of the 3D Earth with Easting/Northing coordinates

A Globe A 3D ellipse with Lat/Long coordinates

The Projection: Project the globe onto a 2D surface

Ellipses
Determining the shape of the earth is the first step in developing a CRS. An ellipse is a simple model describing the basic shape of the Earth. All mapping and coordinate systems are based on this shape. The Earth is almost spherical, however there is a tiny bulge at the equator that makes it ~0.33% larger than at the poles. The ellipsoid is an approximation and does not fit the Earth perfectly. There are different ellipsoids in use, some are designed to fit the whole Earth (WGS84, GRS80) and some are designed to fit a local region (NAD27). Local ellipses can be more accurate for the area they were designed for, but are not useful in other parts of the world. The modern trend is to use a global ellipsoid for compatibility, such as WGS84. The local-best fitting ellipsoid is now considered an old-fashioned concept, but many maps are based on these ellipsoids.

Some common ellipsoid models:

Page 2, Melanie Frazier

Projected
vs.

Unprojected

There are two general options: (1) unprojected (a.k.a. Geographic): Latitude/Longitude for referencing location on the ellipsoid Earth, and (2) projected: Easting/Northing for referencing location on 2D representations of Earth (the creation of maps) To see the options in R: projInfo(type = "proj")

Unprojected/Geographic: Lat/Long
Locations on Earths three-dimensional spherical surface are referenced using Latitude and Longitude. The Latitude and Longitude coordinates for a particular location will differ depending on the CRS and when the measurement was taken. The 3 most common in U.S.: WGS84 (EPSG: 4326) +init=epsg:4326 +proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs +towgs84=0,0,0 ## CRS used by Google Earth and the U.S. Department of Defense for all their mapping. Tends to be used for global reference systems. GPS satellites broadcast the predicted WGS84 orbits. NAD83 (EPSG:4269) +init=epsg:4269 +proj=longlat +ellps=GRS80 +datum=NAD83 +no_defs +towgs84=0,0,0 ##Most commonly used by U.S. federal agencies. Aligned with WGS84 at creation, but has since drifted. Although WGS84 and NAD83 are not equivalent, for most applications they are considered equivalent. NAD27 (EPSG: 4267) +init=epsg:4267 +proj=longlat +ellps=clrk66 +datum=NAD27 +no_defs +nadgrids=@conus,@alaska,@ntv2_0.gsb,@ntv1_can.dat ##Has been replaced by NAD83, but is still encountered!

http://en.wikipedia.org/wiki/North_American_Datum NAD83 vs. WGS84 The initial definition of NAD83(1986) was equivalent to WGS84. Over time, the two systems have diverged, primarily due to refinements (i.e. realizations) made to both systems. The most significant change occurred in 1994 when WGS84 was aligned to a different reference frame. This correction resulted in coordinate differences of up to 1-1.5m with respect to NAD83. Some divergence has also been due to tectonic activity. A fundamental difference between the two reference systems is that NAD83 is intended to track the movements of the North American plate and therefore location coordinates should remain essentially constant over time for this region. WGS84 is referenced to The International Reference Meridian. WGS84 lat/long coordinates are stationary with respect to the average of all global tectonic motions, but this means they are in motion relative to any particular region or country. WGS84 latitude and longitude coordinates are valid for only a specific date (epoch). NAD83 moves at about 2.5 cm/y relative to the ITRF/WGS84 systems. They also use different ellipsoid models (NAD83: GRS80 vs. WGS84: WGS84) but the ellipsoids are very similar and this doesnt cause any meaningful differences. Transforming between NAD83 and WGS84 typically isnt recommended for most applications because standard transformations can introduce error that is large relative to their difference. Complicating the matter, the difference between NAD83 and WGS84 varies with time and location. Both systems have frequent new realizations due to more data and improved techniques. NAD27 vs. NAD83 NAD27 is based on the Clarke Ellipsoid of 1866, whereas NAD83 is based on the GRS 80 ellipsoid. The Clarke Ellipsoid is much less accurate.

Datum shift between NAD27 and NAD83

A point with a given latitude and longitude in NAD27 may be many tens of meters from another point having the identical latitude and longitude in NAD83.

Page 3, Melanie Frazier

Projected data
Projected: Easting/Northing

This would all be a lot simpler if we only used globes. But, alas, they can be inconvenient, so we need maps that can be viewed on paper or computer screens.

Universal Transverse Mercator (UTM)

The UTM projection is commonly used in research because it tends to be more locally accurate, and furthermore, it has attributes that make the estimating distance easy and accurate. Positions are described using Easting and Northing coordinates. The mercator projection preserves angles and direction, but distorts distance. To minimize this distortion, the UTM divides the Earth into sixty zones, and uses a secant transverse Mercator projection in each zone. Over time, and depending on location, UTM coordinates have been based on different ellipsoid models. For areas within the contiguous United States the Clarke 1866 ellipsoid was originally used. The WGS84 ellipsoid is now often used. The difference between UTM coordinates for the Clarke 1986 and current ellipsoids can be over 200 meters. When describing UTM projections, two pieces of information must be provided: (1) the zone and, (2) the ellipsoid model. This information is represented in R as: proj4string(x) <-CRS("+proj=utm +zone=10 +datum=WGS84") According to Wikipedia: Distortion increases in each UTM zone as the boundaries are approached. However, it is often convenient or necessary to use a single grid when locations are located in two adjacent zones. Ideally, the coordinates of each location should be measured on the grid for the zone in which they are located, but it is possible to overlap measurements into an adjoining zone for some distance when necessary. For example, standard convention is to use Zone 10 for all of Oregon, even though Oregon is split between zone 10 and 11.

The elliptical Earth can be projected onto a flat surface (i.e., a paper map). Map coordinates of a point are computed from its ellipsoidal latitude and longitude by a standard formula known as a map projection. It is impossible to flatten a round object without distortion, and this results in trade-offs between area, direction, shape, and distance. For example, there is a trade-off between distance and direction because both features can not be simultaneously preserved. There is no "best" projection, but some projections are better suited to different applications. These websites provide good descriptions of various projections: http://www.radicalcartography.net/?projectionref http://egsc.usgs.gov/isb/pubs/MapProjections/projections.html http://gothos.info/2011/04/common-map-projectiondefinitions/ Projection examples library(maps); library(ggplot2); library(mapproj) states <- map_data("state") usamap <- ggplot(states, aes(x=long, y=lat, group=group)) + geom_polygon(fill="white", colour="black") usamap + coord_map("mercator")

Mercator preserves direction and is useful for navigation. But distances and areas are distorted, especially near the polar regions usamap + coord_map("azequalarea")

Azimuthal Equal Area preserves area, but not direction.