public
Last active

  • Download Gist
state_analysis.R
R
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128
library(maptools)
library(geosphere)
 
# load USA state-level spatial data
# download from http://gadm.org
# click the 'download' tab
# select county = 'united states', file format = 'R', click ok
# download 'level 1' for state-level data
load("USA_adm1.RData")
 
# get coordinate matrix for a state
getCoordinatesByStateName = function(state) {
ix = which(gadm$NAME_1 == state)
sdf = gadm[ix,]
coords = do.call(rbind, lapply(sdf@polygons[[1]]@Polygons, function(p) { p@coords }))
return(coords)
}
 
# take a matrix of coordinates, remove points that are less than minDistance km apart
simplifyCoords = function(coords, minDistance = 10, maxToSkip = 1000) {
ncoords = nrow(coords)
i = 1
simplified = c()
while(i < nrow(coords)) {
upper = min(ncoords, i + maxToSkip)
coordsToConsider = matrix(coords[i:upper,], ncol = 2)
dists = spDistsN1(coordsToConsider, coords[i,], longlat = TRUE)
ixToSelect = which(dists > minDistance)[1]
if(is.na(ixToSelect)) { ixToSelect = length(dists) }
simplified = rbind(simplified, coordsToConsider[ixToSelect,])
i = i + ixToSelect - 1
}
return(simplified)
}
 
# plot a primary state and optionally some surrounding states
plotStates = function(primary, additional = NULL, pcol = '#cccccc', acol = '#eeeeee', ...) {
s1 = which(gadm$NAME_1 == primary)
s2 = which(gadm$NAME_1 %in% additional)
ix = c(s1, s2)
sdf = gadm[ix, ]
nstates = length(ix)
cols = c(pcol, rep(acol, nstates - 1))
plot(sdf, col=cols, lwd = 0.1, ...)
}
 
analyzeState = function(state, minDistance, animate = FALSE, gcPoints = 100) {
fullCoords = getCoordinatesByStateName(state)
simplifiedCoords = simplifyCoords(fullCoords, minDistance)
ncoords = nrow(simplifiedCoords)
maxStates = 0
maxp1 = NA
maxp2 = NA
maxStateNames = NA
for(p1Index in 1:(ncoords - 1)) {
print(p1Index)
p1 = simplifiedCoords[p1Index,]
for(p2Index in (p1Index + 1):ncoords) {
p2 = simplifiedCoords[p2Index,]
line = gcIntermediate(p1, p2, n = gcPoints, addStartEnd = FALSE)
points = SpatialPoints(data.frame(x = line[,1], y = line[,2]),
proj4string=CRS(" +proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs +towgs84=0,0,0"))
states = unique(over(points, gadm)$NAME_1)
states = states[!is.na(states)]
if(length(states) > maxStates) {
maxStates = length(states)
maxp1 = p1
maxp2 = p2
maxStateNames = as.character(states)
}
if(animate) {
plotStates(state, states[states != state])
lines(line)
}
}
}
return(list(n = maxStates, p1 = maxp1, p2 = maxp2, states = maxStateNames))
}
 
# states you want to calculate
stateNames = c('New York', 'Maryland', 'West Virginia')
 
# run the calculations and save results
for(state in stateNames) {
print(state)
result = analyzeState(state, minDistance = 10)
dump("result", file = paste("results/", state, ".R", sep=""))
}
 
# check the results and make some plots!
for(file in dir("results/", pattern = "*.R")) {
source(paste("results/", file, sep=""))
state = strsplit(file, ".", fixed = TRUE)[[1]][1]
otherStates = result$states[result$states != state]
line = gcIntermediate(result$p1, result$p2, n = 100, addStartEnd = TRUE)
png(filename=paste("results/", state, ".png", sep=""), width=480, heigh=480)
plotStates(state, otherStates)
lines(line)
points(line[c(1, nrow(line)),], pch=19, cex=0.5, col='red')
dev.off()
xlim = sort(c(result$p1[1], result$p2[1])) + c(-0.1, 0.1)
ylim = sort(c(result$p1[2], result$p2[2])) + c(-0.1, 0.1)
png(filename=paste("results/", state, "_zoom.png", sep=""), width=480, heigh=480)
plotStates(state, otherStates, xlim = xlim, ylim = ylim)
lines(line)
points(line[c(1, nrow(line)),], pch=19, cex=0.5, col='red')
dev.off()
}

Run the code, it took several hours but not finished yet. I don't know where the problems come from. Could you recheck it?

Please sign in to comment on this gist.

Something went wrong with that request. Please try again.