Joins and Geospatial Data
Intro to Data Science for Biology
March 31, 2016
Map maker, map maper, make me a map!
Today we’re going to look at how we can use joins and geospatial data more explicitly to make maps. Maps are among the first data visualizations that ever occured. And some of the most powerful. They’re also one of the places where joins become incredibly important, as to put data on a map we have to join our data with a geospatial description of the map we want.
Death from Heart Disease
Today’s data set that we’ll be using is a data set of heart disease mortality from the CDC
#load the data and prep
#for some data manipulation
library(readxl)
library(dplyr)
heart_disease <- read_excel("./join_maps/hd_all.xlsx",
na="Insufficient Data")
head(heart_disease)## Source: local data frame [6 x 4]
##
## State County Death_Rate FIPS_Code
## (chr) (chr) (dbl) (dbl)
## 1 Alabama Autauga 463.0 1001
## 2 Alabama Baldwin 391.4 1003
## 3 Alabama Barbour 533.1 1005
## 4 Alabama Bibb 511.1 1007
## 5 Alabama Blount 425.6 1009
## 6 Alabama Bullock 483.2 1011OK, we see that we have state, county, and information on death. FIPS codes, FYI, are standardized county codes. We’ll be ignoring them.
Introducing: Maps
There are a LOT of ways to get map data into R. We’re going to begin with the simplest - grabbing it from an R package. ggplot2 works in tandem with the maps package to provide a few standardized sets of maps for easy plotting. Let’s take a look at one of counties in the U.S. lower 48.
#install the mapdata library if you
#don't have it
library(mapdata)## Loading required package: maps## Warning: package 'maps' was built under R version 3.2.3##
## # ATTENTION: maps v3.0 has an updated 'world' map. #
## # Many country borders and names have changed since 1990. #
## # Type '?world' or 'news(package="maps")'. See README_v3. #library(ggplot2)
#map_data gets us one of the select maps
map_df <- map_data("county")
head(map_df)## long lat group order region subregion
## 1 -86.50517 32.34920 1 1 alabama autauga
## 2 -86.53382 32.35493 1 2 alabama autauga
## 3 -86.54527 32.36639 1 3 alabama autauga
## 4 -86.55673 32.37785 1 4 alabama autauga
## 5 -86.57966 32.38357 1 5 alabama autauga
## 6 -86.59111 32.37785 1 6 alabama autaugaOK, we have latitude and longitude of county borders, a group (each county has one group), and both a region and subregion. Note we don’t have states and counties - this map is a bit broader than that. It includes cities and US Territories. Also note that capitalization is wonky - it’s all lower case.
To show you how we would use this data
ggplot(data=map_df, mapping = aes(x = long, y = lat, group = group)) +
geom_polygon()
geom_polygon is our newest friend, and it creates polygons from the lat/long paths above. There are other ways to use this map, but it’s the clearest.
Getting ready to join maps and data
OK, so, we want to see how death from heart disease varies by county! To do this, we’ll need to join the two data frames. Note, the number of rows in map_df is YUGE compared to those in our data set, as they define county borders. So we’re going to be merging one data point to many rows in the map_df data frame. But, we need to do some prep work first.
First, we need to fix up the columns in our heart_disease data set to not be capitalized, and we need to rename those columns in our map data frame.
map_df <- map_df %>%
rename(State = region,
County = subregion)
heart_disease <- heart_disease %>%
mutate(State = tolower(State),
County = tolower(County))Now let’s test out the join! Before we do the big join, let’s look at mismatch. This is the perfect job for anti_join
bad_match <- anti_join(heart_disease, map_df)## Joining by: c("State", "County")nrow(bad_match)## [1] 112head(bad_match)## Source: local data frame [6 x 4]
##
## State County Death_Rate
## (chr) (chr) (dbl)
## 1 virgin islands of the u.s. saint thomas (county equivalent) NA
## 2 wisconsin st. croix 247.8
## 3 virginia winchester 336.0
## 4 virginia williamsburg 294.2
## 5 virginia waynesboro 364.3
## 6 virginia staunton 339.2
## Variables not shown: FIPS_Code (dbl)Uh oh! 112 rows of bad matches! Why? Well, we can see one of the reasons quickly - the US Virgin Islands. They are not in the map data set. Second, though, we can see st. croix - indeed, the map data frame has no . or ’ characters in it, so, we should filter those out of heart_disease and then re-check.
heart_disease <- heart_disease %>%
mutate(County = gsub("\\.", "", County)) %>%
mutate(County = gsub("\\'", "", County))
bad_match2 <- anti_join(heart_disease, map_df)## Joining by: c("State", "County")nrow(bad_match2)## [1] 82bad_match2## Source: local data frame [82 x 4]
##
## State County Death_Rate
## (chr) (chr) (dbl)
## 1 virgin islands of the u.s. saint thomas (county equivalent) NA
## 2 virgin islands of the u.s. saint john (county equivalent) NA
## 3 virgin islands of the u.s. saint croix (county equivalent) NA
## 4 virginia winchester 336.0
## 5 virginia williamsburg 294.2
## 6 virginia waynesboro 364.3
## 7 virginia staunton 339.2
## 8 virginia salem 388.1
## 9 virginia roanoke city 429.8
## 10 virginia richmond city 378.4
## .. ... ... ...
## Variables not shown: FIPS_Code (dbl)OK, much better - down to 82. And looking at what’s left, those are cities, not counties, so, should not affect our map making.
Are we good the other way?
bad_match3 <- anti_join(map_df, heart_disease)## Joining by: c("State", "County")#because we get a data frame return
length(unique(bad_match3$County))## [1] 8unique(bad_match3$County)## [1] "yellowstone national" "de witt" "de soto"
## [4] "la porte" "du page" "washington"
## [7] "de kalb" "la moure"Only 8. Were we trying to make this perfect, we’d try and figure out why, but, for now, let’s soldier on. 8 is acceptable.
Joining maps and data
So we know that there are some territories and cities left in the heart disease data, and we don’t want to worry about them for the moment. Given that we’ve got some mismatch on both sides, we want to use an inner_join - although we could also merge so that the entire map is retained, and thus see what counties we have missing data for. Try it both ways.
heart_disease_map_data <- inner_join(heart_disease, map_df)## Joining by: c("State", "County")And now - let’s plot it!
heart_map <- ggplot(data=heart_disease_map_data,
mapping = aes(x = long, y = lat, group = group,
fill = Death_Rate)) +
geom_polygon()
heart_map
Making pretty choropleths
So, this is a choropleth map. There are some issues - the scale is hard to resolve, and everything is hard to see. There are a lot of ways we could rectify it. Here are a few tricks.
First, a better scale. Maybe a gradient?
heart_map +
scale_fill_gradient(low = "white", high = "red")
OK, not bad! Not great, but much better. Are there gradients you’d prefer?
However, the bigger problem is that we have a huge range to cover. One common way to make choropleths is to bin data into categories and go from there. Binning in combination with a nice color scale - say via Color Brewer - cab be great. Remember Color Brewer? http://colorbrewer2.org/ actually uses maps as it’s demo! And there’s a scale_fill_brewer function in ggplot2.
So let’s make bins and plot that instead!
heart_disease_map_data <- heart_disease_map_data %>%
mutate(death_rate_bins = cut_interval(Death_Rate, 5))
heart_map_binned <- ggplot(data=heart_disease_map_data,
mapping = aes(x = long, y = lat, group = group,
fill = death_rate_bins)) +
geom_polygon()
heart_map_binned +
scale_fill_brewer(palette="Reds")
MUCH nicer. Now we can really start to see some hot spots.
Is all map data available via packages?
Most map data comes from separate files from R. You’ll see shapefiles most commonly as files that define borders. Shapefiles are funky things. Let’s take a look at the shapefile for a common classification of marine ecoregions of the world. These are coastal regions binned by biogeographic boundaries. To load them, we’ll need the sp and rgdal libraries. These might take some doing to install properly, but, it’s worth it.
library(sp)
library(rgdal)
ecoregions <- readOGR("./join_maps/MEOW-TNC/meow_ecos.shp", layer="meow_ecos")## OGR data source with driver: ESRI Shapefile
## Source: "./join_maps/MEOW-TNC/meow_ecos.shp", layer: "meow_ecos"
## with 232 features
## It has 9 fieldsplot(ecoregions)
Neat - these kinds of files can be plotted on their own! But, we’re going to take this a bit further.
Looking under the hood of a shapefile
So, if we dig into this object, what do we see.
class(ecoregions)## [1] "SpatialPolygonsDataFrame"
## attr(,"package")
## [1] "sp"#It's an S4 object
slotNames(ecoregions)## [1] "data" "polygons" "plotOrder" "bbox" "proj4string"So, it’s a SpatialPolygonsDataFrame. Weird. And it has all of these pieces to it - called slots - as it is what we call an S4 object. More on that another time, but, suffice to say, instead of the $ notation we use an @ notation to access pieces of it. Such as
head(ecoregions@data)## ECO_CODE ECOREGION PROV_CODE
## 0 20192 Agulhas Bank 51
## 1 20053 Aleutian Islands 10
## 2 20072 Amazonia 13
## 3 20194 Amsterdam-St Paul 52
## 4 20228 Amundsen/Bellingshausen Sea 61
## 5 20109 Andaman and Nicobar Islands 24
## PROVINCE RLM_CODE REALM
## 0 Agulhas 10 Temperate Southern Africa
## 1 Cold Temperate Northeast Pacific 3 Temperate Northern Pacific
## 2 North Brazil Shelf 4 Tropical Atlantic
## 3 Amsterdam-St Paul 10 Temperate Southern Africa
## 4 Continental High Antarctic 12 Southern Ocean
## 5 Andaman 5 Western Indo-Pacific
## ALT_CODE ECO_CODE_X Lat_Zone
## 0 189 192 Temperate
## 1 50 53 Temperate
## 2 64 72 Tropical
## 3 191 194 Temperate
## 4 209 228 Polar
## 5 105 109 TropicalHere we see the ecoregion names and some other properties of them. Other slots are filled with information that defines polygons, etc.
Using joins to make ecoregion data frame for plotting
So, to make this into a data frame for mapping, we have to work a littl ebit of magic. Step 1, we need to give the data slot an ID. This is necessary for a join, as the polygons all have an ID as well
slotNames(ecoregions@polygons[[1]])## [1] "Polygons" "plotOrder" "labpt" "ID" "area"ecoregions@data$id = rownames(ecoregions@data)Great! Step 2, we need to turn those polygons into a data frame of points. There’s a function for this in ggplot2, called fortify
ecoregions_points = fortify(ecoregions, ECOREGION="id")## Regions defined for each Polygonshead(ecoregions_points)## long lat order hole piece id group
## 1 28.35993 -36.64435 1 FALSE 1 0 0.1
## 2 28.28350 -36.67695 2 FALSE 1 0 0.1
## 3 28.22727 -36.69682 3 FALSE 1 0 0.1
## 4 28.13197 -36.74429 4 FALSE 1 0 0.1
## 5 28.00760 -36.77946 5 FALSE 1 0 0.1
## 6 27.86775 -36.84124 6 FALSE 1 0 0.1And now the last step - a join! We need to join the original data defining names of ecoregions with the data defining their borders.
ecoregions_df = inner_join(ecoregions_points, ecoregions@data, by="id")Voila! We have ecoregions in a data frame format!
ggplot(data=ecoregions_df,
mapping = aes(x = long, y = lat, group = group)) +
geom_polygon()
Enter Regional Temperature Information
OK, we have some polygons - now what are we going to do with them? One use is to visualize climate change over time. I’ve been working on a project to look at climate change in Ecoregions that contain kelp. So, for each Ecoregion, I’ve gotten all of the temperature data from 1950 - the present from http://www.metoffice.gov.uk/hadobs/hadisst/data/download.html - a great source for SST data.
sst_data <- read.csv("./join_maps/hadsst_regions.csv",
stringsAsFactors=FALSE)
head(sst_data)## ECOREGION DateName tempC Year
## 1 Agulhas Bank 1950.01.16 21.137959 1950
## 2 Aleutian Islands 1950.01.16 3.909404 1950
## 3 Bassian 1950.01.16 14.309411 1950
## 4 Beaufort Sea - continental coast and shelf 1950.01.16 -1.800000 1950
## 5 Cape Howe 1950.01.16 18.455363 1950
## 6 Celtic Seas 1950.01.16 10.399248 1950How can we use this data to see climate change on a map? There are a lot of methods. One simple one is to visualize change by decade. First, we have to use a quick trick to calculate decades - divide by ten, round to a whole number, and then multiple by ten.
sst_decadal <- sst_data %>%
mutate(Decade = round(Year/10)*10) %>%
group_by(ECOREGION, Decade) %>%
summarise(tempC = mean(tempC, na.rm=T)) %>%
ungroup()Now, we can look at temperature by decade, but, we can’t just look at raw temperature. The color palatte would be too spread out, and we couldn’t compare, say, change in the tropics to change in the Arctic. One way that scientists have gotten around this is to look at temperature anomoly - that is, deviation from a long-term mean. So let’s calculate the decadal deviation from the long-term average for each ecoregion.
sst_decadal <- sst_decadal %>%
group_by(ECOREGION) %>%
mutate(tempC_anomoly = tempC - mean(tempC)) %>%
ungroup()OK, now we’ve got something worth plotting!
Visualizing Climate Change
There are a lot of ways we can visualize this. We can look at long-term trends, of course, and even look at each region as a single data point to get average long-term trends
ggplot(data=sst_decadal, mapping=aes(x=Decade, y=tempC_anomoly)) +
geom_line(mapping=aes(group=ECOREGION), colour="grey") +
theme_bw() +
stat_smooth(method="lm", color="red", fill=NA) +
ylab("Decadal Temperature Anomoly")
But how do we put this information on a map so that we can truly see it?
As we did with states, we can join the temperature data to the marine ecoregion data.
sst_map_data <- inner_join(sst_decadal, ecoregions_df)## Joining by: "ECOREGION"## Warning in inner_join_impl(x, y, by$x, by$y): joining character vector and
## factor, coercing into character vectorOnce that is done, we can use this decadal anomoly as a fill, with facets to show different decades.
ggplot(data=sst_map_data,
mapping = aes(x = long, y = lat,
group = group, fill = tempC_anomoly)) +
borders("world", lwd=0.1, colour="black") +
geom_polygon(alpha=0.9) +
facet_wrap(~Decade) +
scale_fill_gradient(low = "blue", high = "red") +
theme_bw(base_size=14) +
coord_equal()
It’s a different way of looking at the same data. But what do you see?