32 Map Projections
A spatial projection is a mathematical representation of a coordinate space used to identify geospatial objects. Because the earth is both non-flat and non-spheroid, we must use mathematical approaches to describe the shape of the earth in a coordinate space. We do this using an ellipsoid—a simplified model of the shape of the earth. Common ellipsoids include:
- NAD27 (North American Datum of 1927) based upon land surveys
- NAD83 based upon satellite data measuring the distance of the surface of the earth to the center of the plant. This is also internationally known as GRS80 (Geodetic Reference System 1980) internationally.
- WGS84 (World Geodetic System 1984) is a refinement of GRS80 done by the US military that was used in the development of GPS systems (and subsequently for all of us).
Onto this ellipsoid, we must define a set of reference locations (in 3-space) called datum that help describe the precise shape of the surface.
32.1 Projections
A projection onto an ellipsoid is a way of converting the spherical coordinates, such as longitude and latitude, into 2-dimensional coordinates we can use. There are three main types of approaches that have been used to develop various projections. (see wikipedia for some example imagery of different projections).
These include:
- Azimuthal: An approach in which each region of the earth is projected onto a plane tangential to the surface, typically at the pole or equator. Cylindrical: This approach projects the surface of the earth onto a cylinder, which is ‘unrolled’ like a large map. This approach ‘stretches’ distances in a east-west fashion, which is why Greenland looks so large…
- Conic: Another ‘unrolling’ approach, though this time instead of a cylinder, it is projected onto a cone.
All projections produce bias in area, distance, or shape (some do so in more than one), so there is no ‘optimal’ projection. To give you an idea of the consequences of these projections, I’ll use the United States map as an example and we can visualize how it is projected onto a 2-dimensional space using different projections.
32.1.1 Equatorial Projections
These are projections centered on the Prime Meridian (Longitude=0)
Mercator Projection
library(maps)
map("state",proj="mercator")
MollWeide Projection
map("state",proj="mollweide")
Gilbert Projection
map("state",proj="gilbert")
Cylequalarea Projection
Some projections require additional parameters, this one is based upon equally spaced and straight meridians, equal area, and true centered on a particular Latitude. I used the centroid of the US.
map("state",proj="cylequalarea",par=39.83)
32.1.2 Azimuth Projections
These projections are centered on the North Pole with parallels making concentric circles. Meridians are equally spaced radial lines.
Orthographic Projection
map("state",proj="orthographic")
Stereographic Projection
map("state",proj="orthographic")
Perspective Projection
Here the parameter is the distance (in earth radii) the observer is looking.
map("state",proj="perspective", param=8)
Gnomonic Projection
map("state",proj="gnomonic")
32.1.3 Polar Conic Projections
Here projections are symmetric around the Prime Meridian with parallel as segments of concentric circles with meridians being equally spaced.
map("state",proj="conic",par=39.83)
map("state",proj="lambert",par=c(30,40))
32.1.4 Miscellaneous Projections
Square Projection
map("state",proj="square")
Hexagon Projection
map("state",proj="hex")
Bicentric Projection
map("state",proj="bicentric", par=-98)
Guyou Projection
map("state",proj="guyou")
There are a lot of ways to project a 3-dimensional surface onto a 2-dimensional representation. Be aware of what you are using and how you are using it when plotting materials.
32.1.5 Reprojecting Rasters
When working with rasters, we can reproject these onto other projections rather easily. Here is an example from the worldclim elevation tile we used previously (see 31.1).
library(raster)
alt <- raster("./spatial_data/alt.tif")
e <- extent( c(-115,-109,22,30) )
baja_california <- crop( alt, e )
baja_california## class : RasterLayer
## dimensions : 960, 720, 691200 (nrow, ncol, ncell)
## resolution : 0.008333333, 0.008333333 (x, y)
## extent : -115, -109, 22, 30 (xmin, xmax, ymin, ymax)
## coord. ref. : +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0
## data source : in memory
## names : alt
## values : -202, 2263 (min, max)We can now project it to another projection, lets say Lambert Conic Conformal.
library(rgdal)
projection(baja_california)## [1] "+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0"baja_lcc <- projectRaster( baja_california, crs="+proj=lcc +lat_1=48 +lat_2=33 +lon_0=-100 +ellps=WGS84")
baja_lcc## class : RasterLayer
## dimensions : 1060, 878, 930680 (nrow, ncol, ncell)
## resolution : 845, 935 (x, y)
## extent : -1610637, -868727, 2790826, 3781926 (xmin, xmax, ymin, ymax)
## coord. ref. : +proj=lcc +lat_1=48 +lat_2=33 +lon_0=-100 +ellps=WGS84
## data source : in memory
## names : alt
## values : -202, 2205.154 (min, max)These two projections influence the region as shown below.
32.1.6 Projecting In GGPlot
As usual, there is probably a way to plot these values in ggplot to make the output just a little bit more awesome. Projections of data in ggplot displays can be manipulated by appending a coord_* object to the plot. Here are two examples using a mercator and azimuth equal area projection of the state maps.
library(ggplot2)
states <- map_data("state")
map <- ggplot( states, aes(x=long,y=lat,group=group))
map <- map + geom_polygon( fill="white",color="black")
map <- map + xlab("Longitude") + ylab("Latitude")
map + coord_map("mercator") 
Conversely, we can plot it using the equal area Azimuth projection
map + coord_map("azequalarea")
or fisheye
map + coord_map("fisheye",par=3)
or any other projection available listed in the mapproject() function.
32.2 Coordinate Systems
In R, we typically are dealing with a combination of data that we’ve collected and that we’ve attained from some other provider. In most GIS applications, the coordinate systems we encounter are either:
- UTM (Universal Transverse Mercator) measuring the distance from the prime meridian for the x-coordinate and the distance from the equator (often called northing in the northern hemisphere) for the y-coordinate. These distances are in meters and the globe is divided into 60 zones, each of which is 6 degrees in width. Geographic coordinate systems use longitude and latitude. For historical purposes these are unfortunately reported in degrees, minutes, seconds, a temporal abstraction that is both annoying and a waste of time (IMHO).
- Decimal degrees, while less easy to remember, are easier to work with in R.
- State Planar coordinate systems are coordinate systems that each US State has defined for their own purposes. They are based upon some arbitrarily defined points of reference and another pain to use (IMHO). Given the differential in state area, some states are also divided into different zones. Maps you get from municipal agencies may be in this coordinate system. If your study straddles different zones or even state lines, you have some work ahead of you…
It is best to use a system that is designed for your kind of work. Do not, for example, use a state plane system outside of that state as you have bias associated with the distance away from the origin. That said, Longitude/Latitude (decimal degrees) and UTM systems are probably the easiest to work with in R.
In R, we use rgdal to project points. Here I load in the coordinates of the populations in the Arapatus attenuatus data set and make a SpatialPoints object out of it. Setting the proj4string() here does not project the data, I am just specifying that the data are already in the lat/long WGS84 format.
library(sp)
library(gstudio)
data(arapat)
coords <- strata_coordinates( arapat )
pts <- SpatialPoints( coords[,2:3] )
proj4string(pts) <- CRS("+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0")
pts## class : SpatialPoints
## features : 39
## extent : -114.2935, -109.1263, 23.0757, 29.32541 (xmin, xmax, ymin, ymax)
## coord. ref. : +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0The CRS() function holds the definition of the projection and interfaces between the PROJ.4 and RGDAL libraries. To project a set of data points into a new coordinate systems, we use spTransform() and pass it the definition of the new system to use.
32.2.1 Changing Datum
We will first begin by looking at differences in the actual datum used to record the loation of plots. Here I compare the decimal Longitude/Latitude we’ve used thus far with that from the Universal Transverse Mercator (UTM).
pts.utm <- spTransform(pts, CRS("+proj=utm +zone=12 +datum=WGS84"))
summary( pts.utm )## Object of class SpatialPoints
## Coordinates:
## min max
## Longitude 180128 686925.2
## Latitude 2552540 3248545.0
## Is projected: TRUE
## proj4string :
## [+proj=utm +zone=12 +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0]
## Number of points: 39You can see the tranformations in the coordinate system by comparing the plots below. The relative position of each point is the same.
