Skip to content

Instantly share code, notes, and snippets.

# toddwschneider/state_analysis.R

Last active December 10, 2021 15:15
Show Gist options
• Save toddwschneider/6678635 to your computer and use it in GitHub Desktop.
This file contains 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
 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() }

### nghiaht commented Dec 27, 2013

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

to join this conversation on GitHub. Already have an account? Sign in to comment