Map projections

Projecting the surface of a three dimensional object onto a two dimensional surface produces distortions in the area, shape or distance between spatial features. Here I demonstrate how to visualize the impacts that projections can have on area, shape and distance using Tissot’s indicatrix.

R packages such as sf offer us the ability to specify projections by linking to the PROJ package. PROJ is a library that enable us to quickly and conveniently switch between a large number of different projections.

#Loading libraries
library(ggplot2)   # For generating the map visuals.
library(sf)        # For handling spatial objects, transformations etc.
library(dplyr)     # Data manipultaion and wrangling
library(rworldmap) # Data source for world map

First we need to import the map of the world

world<- getMap()%>%            #Importing map file
  st_as_sf()%>%                #Converting it into an 'sf' object
  st_transform(., crs = 4326)  #Transforming the coordinate reference system.

Tissot’s Indicatrices

Now we must set up our grid of Tissot’s indicatices. These are a grid of circles of equal size that are evenly spaced across longitude and latitude. They help us visualize and appreciate the distortion in shapes, areas, and distances that can occur through projection. This visualization technique was first presentd by the French mathematician Nicolas Auguste Tissot during the late 19th century.

lon <- seq(-135,135, 45) # Defining the longitiudes (vertical grid lines)
lat <- seq(-70,70,35)    # Defining the latitudes (horizontal grid lines)
         
#Generating coordinates for each time combindation of lat and long
latlon<- lapply(lat, data.frame, lon)%>%  
  bind_rows()
names(latlon) <- c("lat", "lon")
  
#Generating the circles to be used for our Tissot's indicatrices
tissots <-latlon%>%
  st_as_sf(., coords = c("lon", "lat"))%>% #Converting our latlon into `sf`
  st_set_crs(4326)%>%                      #Giving them a coordinate reference system
  st_buffer(., dist = units::set_units(500, "km")) #Generating circles with a 500km radius from our points.

head(tissots)
## Simple feature collection with 6 features and 0 fields
## Geometry type: POLYGON
## Dimension:     XY
## Bounding box:  xmin: -148.4014 ymin: -74.53223 xmax: 103.4179 ymax: -65.46998
## Geodetic CRS:  WGS 84
##                         geometry
## 1 POLYGON ((-148.2988 -70.441...
## 2 POLYGON ((-85.74652 -65.791...
## 3 POLYGON ((-47.35994 -65.573...
## 4 POLYGON ((-1.080087 -74.496...
## 5 POLYGON ((47.86902 -74.4485...
## 6 POLYGON ((80.50167 -67.0979...

Now we are ready to put the map and grid of circles together to visualise and explore the impacts projections have on area,shape and distance.

First lets start with the orthographic projection in which the coordinates are projected onto a sphere. Note the uniformity of the circles across latitude and longitude.

We use the coord_sf() function to set the projection we wish to plot.

ggplot()+
  geom_sf(data = world)+
  geom_sf(data = tissots)+
  coord_sf( crs= "+proj=ortho +lat_0=50 +lon_0 = 310")

Mercator

This is a commonly used map projection and one you are likely familiar with if you’ve zoomed out enough on google maps. Notice that the shape has been preserved, but the areas of the landmasses has dramatically expanded at higher latitudes (further north and south) and shrunk at lower latitudes (closer to the equator).

ggplot()+
  geom_sf(data = world%>%
            filter(continent != "Antarctica"))+ #Antarctica is removed as it results in a map that is too distorted.
  geom_sf(data = tissots)+
  coord_sf(crs = "+proj=merc")

Interrupted Goode Homolosine

This projection does well to preserve the shapes and areas overall, but in some areas results in severe distortion at certain locations, with the only large landmass notably impacted being greenland.

ggplot()+
  geom_sf(data = world)+
  geom_sf(data = tissots)+
  coord_sf(crs = "+proj=igh")

Azimuthal Equidistant

Area, size and distance are well preserved at 0 longitude but the distortion increases with increasing longitude.

ggplot()+
  geom_sf(data = world)+
  geom_sf(data = tissots)+
  coord_sf(crs = "+proj=aeqd")

Nicolosi Globular

This projection was invented over 1000 years ago by the Iranian polymath al-Biruni and gained popularity during the 17th century as it preserves the shape, area, and distances within Europe very well.

ggplot()+
  geom_sf(data = world)+
  geom_sf(data = tissots)+
  coord_sf(crs = "+proj=nicol")

These represent just a handful of the many hundreds of map projections available, each of which has its own merits and particular use cases. Next time you wish to display spatial data on a global or continental scale take a moment. Think about whether the area, shape, or apparent distances between the features you are mapping could be influenced by the projection you have chosen and whether this could influence how the interpretation.

Projections available can be viewed using the sf_proj_info() function available within the sf package.

For an interactive way to explore different projections: https://alanw.shinyapps.io/tissotsinidcatrix/

All code used to generate this document: https://github.com/AWard93/RPubs_spatial

For more information on PROJ: https://proj.org/en/9.3/about.html

Please feel free to take and use all code presented here. Please comment if you see any improvements that can be made to the code or have anything you would like to add or think I have missed.