Created
August 5, 2016 03:33
-
-
Save dill/49b2ae223872ff736bf43687bda193fe to your computer and use it in GitHub Desktop.
Is Columus, OH the only place in the US where 1/2 the population is within 500 miles? Well, kinda
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
# is Columbus, OH really the only place in the US where you're within | |
# 500 miles of half of the US population? | |
# claim: https://448f59f74df57015bbb8-a9447b7dfa4ae38e337b359963d557c4.ssl.cf3.rackcdn.com/9987%20Brewdog%20EFP%20Prospectus%20USA%20A4%20v9.pdf | |
# "Columbus is within 500 miles of over half of the US population" | |
# get the US 2010 census data | |
# http://www2.census.gov/geo/docs/reference/cenpop2010/county/CenPop2010_Mean_CO.txt | |
cent <- read.csv("CenPop2010_Mean_CO.txt") | |
# make that data spatial | |
library(sp) | |
coords <- cbind(cent$LONGITUDE, cent$LATITUDE) | |
sp <- SpatialPoints(coords) | |
spdf <- SpatialPointsDataFrame(sp, cent) | |
proj4string(spdf) <- CRS("+proj=longlat +datum=WGS84") | |
# reproject | |
lcc_proj4 <- CRS("+proj=lcc +lat_1=20 +lat_2=60 +lat_0=40 +lon_0=-96 +x_0=0 + y_0=0 +ellps=GRS80 +datum=NAD83 +units=m +no_defs ") | |
# project using LCC | |
projd <- spTransform(spdf, CRSobj=lcc_proj4) | |
# distances between county centroids | |
distm <- as.matrix(dist(coordinates(projd), upper=TRUE, diag=TRUE)) | |
# iterate over the counties, finding the populations within 500 miles | |
winners <- c() | |
for(i in 1:nrow(distm)){ | |
ind <- which(distm[i,] < 804672) # 500 miles in meters | |
# 2010 US population was 309.3 million, there are no counties | |
# with > 150 million within 500 miles, so let's say 148 million | |
if(sum(cent$POPULATION[ind]) > 148e6){ | |
cat("HURRAH:", as.character(cent$COUNAME)[i], as.character(cent$STNAME)[i], "is a winner!\n") | |
winners <- c(winners, i) | |
} | |
} | |
par(mfrow=c(1,2)) | |
# plot where that town is | |
plot(projd, xlim=c(-2.1e6,2e6), ylim=c(-2e6,1.5e6)) | |
points(coordinates(projd)[winners,,drop=FALSE],pch=19,col="red") | |
# relax that a bit and look for town with >140 million within 500 miles | |
winners <- c() | |
for(i in 1:nrow(distm)){ | |
ind <- which(distm[i,] < 804672) # 500 miles in meters | |
# 2010 US population was 309.3 million, there are no counties | |
# with > 150 million within 500 miles, so let's say 148 million | |
if(sum(cent$POPULATION[ind]) > 140e6){ | |
cat("HURRAH:", as.character(cent$COUNAME)[i], as.character(cent$STNAME)[i], "is a winner!\n") | |
winners <- c(winners, i) | |
} | |
} | |
# plot where that town is | |
plot(projd, xlim=c(-2.1e6,2e6), ylim=c(-2e6,1.5e6)) | |
points(coordinates(projd)[winners,,drop=FALSE],pch=19,col="red") | |
#> unique(as.character(cent$STNAME)[winners]) | |
#[1] "Kentucky" "Maryland" "Ohio" "Pennsylvania" | |
#[5] "Virginia" "West Virginia" | |
Folks from the (albeit population small) non-continential portion of the U.S. might find problems with your analysis.
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment
plots...
