county_dendist_map.Rmd

---
title: "County Density/Distance/Population Map"
author: "Jeff Erickson"
date: "Tuesday, June 17, 2014"
output: html_document
---

### LOAD REQUIRED LIBRARIES

```{r}
library(rgdal)
library(ggmap)
library(ggplot2)
library(rgeos)
library(knitr)
library(sqldf)
```

### IMPORT DATA

```{r}
#County Shapes from US Census
county.shapes <- readOGR(dsn="datasets", "gz_2010_us_050_00_500k")

#County population centroids from US Census
county.pop.center <- read.csv("datasets/CenPop2010_Mean_CO.txt", header=TRUE)

#County population density information from US Census
county.density <- read.csv("datasets/DEC_10_SF1_GCTPH1.CY07_with_ann.csv", header=TRUE, as.is=TRUE)

#Major city (primary statistical areas) from gen.maj.city.long.lat.Rmd
maj.city.long.lat <- read.csv("datasets/maj.city.long.lat.csv", header=TRUE)

#Google driving distances from https://github.com/jefferickson/county-city-driving-dist
county.city.driving.dist <- read.csv("datasets/county-city-driving-dist.csv", header=TRUE, as.is=TRUE)
```

### HELPER FUNCTIONS/DEFINITIONS

```{r}
missoula.fips <- "0500000US30063"
westchester.fips <- "0500000US36119"

deg.2.rad <- function(x) { #Convert vector from degrees to radians
  return(sapply(x, function(deg) {
    return(deg * (pi/180))
  }))
}

euclid.dist <- function(x) { #Simple Euclidean distance
  #arg: c(lat1, long1, lat2, long2)
  return(sqrt((x[1]-x[3])^2 + (x[2]-x[4])^2))
}

haversine.dist <- function(x) { 
  #Implementation of the haversine formula
  #http://en.wikipedia.org/wiki/Haversine_formula
  #http://stackoverflow.com/questions/27928/how-do-i-calculate-distance-between-two-latitude-longitude-points
  #arg: c(lat1, long1, lat2, long2)
  coords <- deg.2.rad(x)
  earth.r <- 6371 #earth radius in km
  d.long.rad <- coords[2] - coords[4]
  d.lat.rad <- coords[1] - coords[3]

  a <- sin(d.lat.rad/2)^2 + cos(coords[1]) * cos(coords[3]) * sin(d.long.rad/2)^2
  c <- 2 * atan2(sqrt(a), sqrt(1 - a))
  d <- earth.r * c
  return(d)
}

google.driving.dist <- function(x) {
  #Get Google driving distance
  #http://stackoverflow.com/questions/16863018/getting-driving-distance-between-two-points-lat-lon-using-r-and-google-map-ap
  #Has the following limitations:
  #100 elements per query.
  #100 elements per 10 seconds.
  #2500 elements per 24 hour period.
  #arg: c(lat1, long1, lat2, long2)
  library(XML)
  library(RCurl)
  origin <- paste(x[1], x[2], sep=",")
  dest <- paste(x[3], x[4], sep=",")
  url <- paste0('http://maps.googleapis.com/maps/api/distancematrix/xml?origins=', origin, '&destinations=', dest, '&mode=driving&sensor=false')
  xml.resp <- xmlParse(getURL(url))
  d.lab <- xmlValue(xmlChildren(xpathApply(xml.resp, "//distance")[[1]])$text)
  d <- as.numeric(sub(" km", "", d.lab))
  return(d)
}

z.score <- function(x) { #Calc z-scores
  #arg: a vector of values to z-score
  return((x - mean(x)) / sd(x))
}

link.func.unwgted.avg <- function(x) { #Unweighted average of each argument (row-by-row), then scaled to [0,1]
  #arg: data.frame of numerics to be averaged and scaled
  index.unscaled <- apply(x, 1, function(x) { return(mean(x)) })
  index.shifted <- index.unscaled - min(index.unscaled)
  index.scaled <- index.shifted / max(index.shifted)
  return(index.scaled)
}
```

### DATA TRANSFORMATIONS

```{r}
#We only want counties, not census tracts
county.density <- county.density[which(county.density$GEO.id2 == county.density$GCT_STUB.target.geo.id2), ]

#Create key that is used to merge
county.pop.center$GEO_ID <- paste("0500000US", formatC(county.pop.center$STATEFP, width=2, flag="0"), formatC(county.pop.center$COUNTYFP, width=3, flag="0"), sep="")

#For now, only the lower 48 states
county.density <- county.density[which(! (substr(county.density$GEO.id, 10, 11) %in% c("02", "52", "15", "72"))), ] #remove AK, HI, PR (2 codes)
county.pop.center <- county.pop.center[which(! (substr(county.pop.center$GEO_ID, 10, 11) %in% c("02", "52", "15", "72"))), ] #remove AK, HI, PR (2 codes)

#Convert field(s) to numeric. Remove (r[number]) pattern found in some HD01 values.
county.density$SUBHD0401 <- as.numeric(county.density$SUBHD0401)
county.density$HD01 <- as.numeric(gsub("\\(r[0-9]*\\)", "", county.density$HD01))

#clean up the driving distance data
county.city.driving.dist <- county.city.driving.dist[c("fips", "city_id", "driving_distance")]
```

### CALCULATE DISTANCES AND DEN-DIST-POP INDEX

```{r}
#Create a list of all unique combinations of counties and the major cities
county.city.list.long <- merge(county.pop.center$GEO_ID, maj.city.long.lat$id, all.x=TRUE, all.y=TRUE)
names(county.city.list.long) <- c("county.geo.id", "city.id")

#For each combination, calculate the distance
county.city.list.long$dist <- apply(county.city.list.long[c("county.geo.id", "city.id")], 1, function(x) {
  dist <- haversine.dist(as.numeric(
    c(
      county.pop.center[which(county.pop.center$GEO_ID == x[1]), c("LATITUDE", "LONGITUDE")],
      maj.city.long.lat[which(maj.city.long.lat$id == x[2]), c("lat", "long")]
    )
  ))
  return(dist)
})

#Merge on the county densities/populations
county.city.list.long <- merge(county.city.list.long, county.density[c("GEO.id", "SUBHD0401", "HD01")], by.x="county.geo.id", by.y="GEO.id", all.x=TRUE)

#Merge the driving distances
county.city.list.long <- merge(county.city.list.long, county.city.driving.dist, by.x=c("county.geo.id", "city.id"), by.y=c("fips", "city_id"))

#Find the shortest of the distances (driving and otherwise)
county.city.list.shortest <- sqldf("select `county.geo.id` as county_geo_id, SUBHD0401, HD01, min(dist) as dist, min(driving_distance) as driving_dist from `county.city.list.long` group by `county.geo.id`")

#transform density and population with log
county.city.list.shortest$den.log <- log(county.city.list.shortest$SUBHD0401)
county.city.list.shortest$pop.log <- log(county.city.list.shortest$HD01)

#let's "flip" the distances so that high values of both density and "distance" mean more urban than rural
county.city.list.shortest$dist.flip <- -(county.city.list.shortest$dist - (max(county.city.list.shortest$dist) + 1))
county.city.list.shortest$drive.dist.flip <- -(county.city.list.shortest$driving_dist - (max(county.city.list.shortest$driving_dist) + 1))

#z-score transforms of density and population and distance-flip
county.city.list.shortest$den.log.z <- z.score(county.city.list.shortest$den.log)
county.city.list.shortest$pop.log.z <- z.score(county.city.list.shortest$pop.log)
county.city.list.shortest$dist.flip.z <- z.score(county.city.list.shortest$dist.flip)
county.city.list.shortest$drive.dist.flip.z <- z.score(county.city.list.shortest$drive.dist.flip)

#create the den-dist-pop index
county.city.list.shortest$den.dist.pop.index <- link.func.unwgted.avg(county.city.list.shortest[c("den.log.z", "pop.log.z", "dist.flip.z")])
county.city.list.shortest$den.drivedist.pop.index <- link.func.unwgted.avg(county.city.list.shortest[c("den.log.z", "pop.log.z", "drive.dist.flip.z")])

#what is the difference between the two?
county.city.list.shortest$index.diff <- county.city.list.shortest$den.dist.pop.index - county.city.list.shortest$den.drivedist.pop.index

#create a {-1, 0, 1} flag for index value relative to Missoula County
county.city.list.shortest$rel.missoula <- apply(county.city.list.shortest, 1, function(x) {
  if (x["county_geo_id"] == missoula.fips) {
    return(0)
  } else if (x["den.drivedist.pop.index"] < county.city.list.shortest$den.drivedist.pop.index[which(county.city.list.shortest$county_geo_id == missoula.fips)]) {
    return(-1)
  } else if (x["den.drivedist.pop.index"] > county.city.list.shortest$den.drivedist.pop.index[which(county.city.list.shortest$county_geo_id == missoula.fips)]) {
    return(1)
  } else {
    return(NA)
  }
})
county.city.list.shortest$rel.missoula <- factor(county.city.list.shortest$rel.missoula, levels=c(-1, 0, 1), labels=c("More Rural than Missoula", "Missoula", "Less Rural than Missoula"))
```

### MERGE DEN-DIST-POP DATA WITH MAPPING INFO

```{r}
#Fortify the shapefile and use the FIPS codes as id
county.shapes.f <- fortify(county.shapes, region="GEO_ID")

#Merge on our demo data
county.shapes.f.dendistpop <- merge(county.shapes.f, county.city.list.shortest, by.x="id", by.y="county_geo_id")
```

### GRAPH THE DENSITY OF THE INDEX
```{r}
ggplot(county.city.list.shortest, aes(x=den.dist.pop.index)) + geom_density()

ggplot(county.city.list.shortest, aes(x=den.drivedist.pop.index)) + geom_density()
```

### MAP IT!

```{r}
ggplot(county.shapes.f.dendistpop[order(county.shapes.f.dendistpop$order), ], aes(x=long, y=lat, group=group)) + #careful: order matters here
  geom_polygon(aes(fill=den.dist.pop.index), colour="black") +
  coord_equal() +
  scale_fill_gradientn(values=c(0, median(county.shapes.f.dendistpop$den.dist.pop.index), 1), colours=c("red", "white", "green"), rescaler=function(x, ...) x, oob=identity) 

ggplot(county.shapes.f.dendistpop[order(county.shapes.f.dendistpop$order), ], aes(x=long, y=lat, group=group)) + #careful: order matters here
  geom_polygon(aes(fill=den.drivedist.pop.index), colour="black") +
  coord_equal() +
  scale_fill_gradientn("Den-Dist-Pop\nIndex", values=c(0, median(county.shapes.f.dendistpop$den.drivedist.pop.index), 1), colours=c("red", "white", "green"), rescaler=function(x, ...) x, oob=identity) +
  ggtitle(expression(atop("Density-Distance-Population Index", atop("By County, Compared to Median", "")))) +
  theme(axis.title=element_blank(), axis.text=element_blank(), axis.ticks=element_blank(), panel.grid=element_blank())

ggplot(county.shapes.f.dendistpop[order(county.shapes.f.dendistpop$order), ], aes(x=long, y=lat, group=group)) + #careful: order matters here
  geom_polygon(aes(fill=index.diff), colour="black") +
  coord_equal() + 
  scale_fill_gradientn("Difference", values=c(-0.07, 0, 0.07), colours=c("green", "white", "red"), rescaler=function(x, ...) x, oob=identity) +
  ggtitle(expression(atop("Driving Distance vs. \"As the Crow Flies\"", atop("Index Difference", "")))) +
  theme(axis.title=element_blank(), axis.text=element_blank(), axis.ticks=element_blank(), panel.grid=element_blank())
```

### CONCLUSIONS

```{r}
county.w.index <- merge(county.city.list.shortest[c("county_geo_id", "den.dist.pop.index", "den.drivedist.pop.index")], county.density[c("GEO.id", "GEO.display.label")], by="county_geo_id", by.y="GEO.id")
```

```{r}
ggplot(county.city.list.shortest, aes(x=den.drivedist.pop.index)) + 
  geom_density() + 
  geom_vline(xintercept=county.city.list.shortest$den.drivedist.pop.index[which(county.city.list.shortest$county_geo_id == missoula.fips)], linetype="longdash", colour="red") + 
  geom_text(x=county.city.list.shortest$den.drivedist.pop.index[which(county.city.list.shortest$county_geo_id == missoula.fips)] + 0.02, y=1, label="MISSOULA COUNTY, MT", colour="red", angle=90) +
  geom_vline(xintercept=county.city.list.shortest$den.drivedist.pop.index[which(county.city.list.shortest$county_geo_id == westchester.fips)], linetype="longdash", colour="red") +
  geom_text(x=county.city.list.shortest$den.drivedist.pop.index[which(county.city.list.shortest$county_geo_id == westchester.fips)] + 0.02, y=1.1, label="WESTCHESTER COUNTY, NY", colour="red", angle=90) +
  ggtitle(expression(atop("Density-Distance-Population Index", atop("Density Plot", "")))) +
  scale_x_continuous("Den-Dist-Pop Index") +
  scale_y_continuous("Density (%)")
```


```{r}
kable(head(county.w.index[order(county.w.index$den.dist.pop.index), ]), format="markdown", row.names=FALSE)
kable(head(county.w.index[order(county.w.index$den.drivedist.pop.index), ]), format="markdown", row.names=FALSE)
```

```{r}
kable(tail(county.w.index[order(county.w.index$den.dist.pop.index), ]), format="markdown", row.names=FALSE)
kable(tail(county.w.index[order(county.w.index$den.drivedist.pop.index), ]), format="markdown", row.names=FALSE)
```

```{r}
ggplot(county.shapes.f.dendistpop[order(county.shapes.f.dendistpop$order), ], aes(x=long, y=lat, group=group)) + #careful: order matters here
  geom_polygon(aes(fill=rel.missoula), colour="black") +
  coord_equal() +
  scale_fill_discrete("Relation to Missoula") +
  ggtitle(expression(atop("Density-Distance-Population Index", atop("By County, Relation to Missoula County", "")))) +
  theme(axis.title=element_blank(), axis.text=element_blank(), axis.ticks=element_blank(), panel.grid=element_blank())
```