Module 52 Geographic Information Systems

Learning goals

  • Understand the difference between vector and raster data, and when to use which.
  • Learn how to work with vector and raster data in base maps in R.
  • Learn how to work with vector and raster data in leaflet maps.

 

For some good examples of bad maps – the kind of thing we are trying to help you avoid – check out this Twitter feed.

We don’t often read scatter plots in life, but we read maps all the time. Learning how to work with spatial data and create beautiful maps are essential skills as a data scientist. To get there, you need to become familiar with Geographical Information Systems, or GIS.

In a GIS, spatial data are tied together in a database that makes it (1) easy to map the data and (2) easy to access all sorts of information underlying each spatial feature.

GIS matters because most data, it turns out, have an important geographic component. There are nearly always important geographic questions you can ask about your data.

Geographic data structures

In your work with spatial data, you are going to encounter two main types of data:

Vectors

Vectors are discrete shapes built up by individual points.

There are three common types of vectors:

  1. A point is a single location (zero-dimensional).

  2. A line is a set of points connected by straight lines (1-dimensional).

  3. A polygon is a shape formed by a line whose beginning and end are the same location: it is a line that begins and ends at the same spot and thus encloses an area. Polygons have an inside and an outside.

Vectors are best used to represent discrete objects or areas: the location of an observation, a river, a road, a building, a country’s boundary, etc.

Most vectors have two types of data: geographic information, such as coordinates, and tabular data with attributes for the vector. For this reason, vectors are usually saved in shapefiles, which typically require a folder of several different files to use. That’s why most shapefiles are downloaded as zipped folders.

Rasters

Rasters are grids. They have rows and columns, just like a dataframe, and each grid cell represents a single value. Rasters are best used for continuous phenomena – in other words, things that occur everywhere: temperature, elevation, humidity, etc.

We interact with rasters all the time: think of smartphone screens and digital photos. Such images are simply rasters of pixels. Each pixel has a value representing color and brightness.

Getting to work

To use R as a GIS, first download a bunch of libraries that are likely to be useful:

Some learning exercises with Tennessee zip codes

We’re going to do some exercises with Tennessee zip codes.

  1. What is a zip code?

  2. What kind of data format is zip code data likely to be in?

  3. Load up some tabular data on zip codes by running the below:

  1. Inspect the data. What is the unit of observation.

  2. Make a histogram of population.

  3. Which zip code is the most populous?

  4. Which city has the most zip codes?

  5. What is the total population of Tennessee?

  6. Make a chart of the 5 most populous cities in Tennessee, showing their population.

  7. Load up some geographic data on Tennessee zip codes by running the below:

  1. Plot the map using the plot function

  2. Run the below in order to convert the format to a ggplot-friendly dataframe.

  1. Make a ggplot of the tng object.

  2. In tnz, make a new variable named percent_water which shows the percentage of each zip code that is water (in area).

  3. Bring data from tnz into tng by running the below:

  1. Modify the above ggplot so that the filled color of each polygon reflects the new variable, percent_water.

  2. Make a leaflet map with no polygons on it.

  3. Add the Tennessee zip codes to the leaflet map.

  4. Make the width of the lines thinner and make the background totally transparent.

  5. Bring the pop column of tn into tnz using left_join.

  6. We want to make a map showing population density. We’ll start by making a color palette:

  1. Make a choropleth of population in leaflet.

  1. Make a choropleth of percentage water in leaflet. Make it blue.

  2. Create a dataframe of the “centroid” of each zip code by running the below:

  1. Define the location, in longitude, latitude, of Sewanee
  1. Use distm from the geosphere package to get the distance from Sewanee to each zip code’s centroid.
  1. Deposit the distances object as a column in the tnz object
  1. Make a leaflet choropleth of distance to Sewanee.

Working with rasters

To download some raster data, user the function getData() from the raster package:

This function is great for getting country-level data. Use ccodes() to get the codes for each country:

                   NAME ISO3 ISO2       NAME_ISO       NAME_FAO
1           Afghanistan  AFG   AF    AFGHANISTAN    Afghanistan
2 Akrotiri and Dhekelia  XAD <NA>           <NA>           <NA>
3                 Åland  ALA   AX  ÅLAND ISLANDS           <NA>
4               Albania  ALB   AL        ALBANIA        Albania
5               Algeria  DZA   DZ        ALGERIA        Algeria
6        American Samoa  ASM   AS AMERICAN SAMOA American Samoa
             NAME_LOCAL      SOVEREIGN       UNREGION1 UNREGION2 continent
1           Afghanestan    Afghanistan   Southern Asia      Asia      Asia
2 Akrotiri and Dhekelia United Kingdom    Western Asia      Asia      Asia
3                 Åland        Finland Northern Europe    Europe    Europe
4             Shqiperia        Albania Southern Europe    Europe    Europe
5            Al Jaza'ir        Algeria Northern Africa    Africa    Africa
6        American Samoa  United States       Polynesia   Oceania   Oceania

Since the United States consists of several distant territories (Hawaii, Alaska, Samoa, etc.), the object usa is actually a list containing an object for each territory:

Let’s plot the continental United States:

How do you interpret this map? What type of data does this raster object contain?

Now plot a separate territory:

Combining rasters & vectors

Now let’s try to combine raster data, such as elevation (stored in the object usa[[1]]), with vector data, such as state boundaries (stored in the object tn).

To do so, let’s crop the elevation raster to only the data contained with the Tennessee boundary:

Dealing with projections

Let’s load in some data specifically for the area of Sewanee, TN:

This zipped folder has data files for both rasters and vectors.

Read in the raster data for elevation:

Now let’s try to map Sewanee’s raster data onto the tn vector of boundaries:

Oh no! That didn’t work. Why not? Have a look at coordinates:

It seems like these two objects are using different coordinate systems. Let’s check out what projection these objects are using:

The tn coordinates are in classic lat/long, but the elevation coordinates appear to be in UTM (Easting/Northing coordinates).

Let’s convert elevation’s coordinate system to than of tn:

Great – now let’s try out plot again:

Phew – worked this time (even though it still ain’t pretty).

Now let’s bring in some vectors pertaining to the Sewanee Domain: boundaries, roads, and structures:

This time let’s make a zoomed-in plot of Sewanee’s land, with elevation and the Domain boundary:

Uh oh. Same problem as before. We need to “reproject” boundary to latitude longitude.

This time, since we are now transforming a vector, we need to use the function spTransform() instead of projectRaster().

Okay, let’s retry our plot:

Let’s use crop and mask to get just the elevation for the domain.

If we wanted more information about this domain_elev raster, we could use some special exploration functions for rasters:

Now let’s add roads to the plot. This time, we will preemptively re-project the roads data.

What is we only want to use the roads that cross the Domain boundary? To do so, we will user the over() function:

Now let’s re-do the plot:

Combining rasters, vectors, & leaflet

Let’s place Sewanee’s elevation raster onto a leaflet map.

First, let’s make a color palette (i.e., a chloropleth) to represent elevation:

We add the raster to the leaflet map using the addRasterImage() function:

Excercises

Adding to the Sewanee leaflet map

1. Add structures to your leaflet map of elevation_ll from above. Hint: you’ll need to use addPolygons and you’ll need to reproject structures as structures_ll…

2. Add popups to your structures.

3. Make your structures a different color and remove the border *Hint:8 you’ll need to use the stroke argument.

4. Get the elevation of each structure by running:

4. Use the structure_elevation object to add a new column to the structures_ll object.

5. Make a histogram of the elevation of Sewanee buildings.

6. How low is the lowest building on the domain? How low is it?

7. How high is the highest building on the domain? Which building is it?

 

A global map of world health

Let’s read in a new dataset of polygons, this time the national boundaries of the world’s nations:

OGR data source with driver: ESRI Shapefile 
Source: "/tmp/world_small", layer: "TM_WORLD_BORDERS_SIMPL-0.3"
with 246 features
It has 11 fields
Integer64 fields read as strings:  POP2005

Now let’s read in some indicator data of global health:

8. Check out this data:

9. Let’s focus in on the rate of inpatient care use among adults:

10. Now let’s join this tabular data to the shapefile data according to country:

11. OK, let’s get a map started. First, let’s establish our base map:

12. Next, let’s create a color scale for the value of interest:

13. And use this palette to color-code each country:

14. Now add a legend:

15. Let’s stage some fancy text:

 

16. Now make a choropleth map of BMI for men, where the darker the shade of red, the higher the BMI for each country.

17. Remove the borders from the map

18. Add a legend on the top right of the map

19. Make the NA color blue

20. Make the hover label a combination of the country and BMI value

21. Make the title of the legend “BMI”

22. Create a function that takes an indicator name as an input and creates a map.