Skip to content

Instantly share code, notes, and snippets.

@potterzot

potterzot/latlon2shape.R

Created December 19, 2017 21:07
Show Gist options
  • Star 0 You must be signed in to star a gist
  • Fork 0 You must be signed in to fork a gist
  • Save potterzot/811d2cb6ee4f2f3857f60158f1ddc7e0 to your computer and use it in GitHub Desktop.
Save potterzot/811d2cb6ee4f2f3857f60158f1ddc7e0 to your computer and use it in GitHub Desktop.
Download ZIP
Map a lat/lon point to a county or other geography
Raw
latlon2shape.R
#' Get county/state name from lat/long pair.
#'
#' Maps a lat/long pair to a geography.
#'
#' @importFrom maps map
#' @importFrom maptools map2SpatialPolygons
#' @importFrom sp over CRS SpatialPoints
#' @export
#'
#' @param lat list of numbers or number of the latitude in degrees.
#' @param lon list of numbers or number of the longitude in degrees.
#' @param projargs a CRS string. Default is "+proj=longlat +datum=WGS84".
#' @param shape either "county" or "state".
#' @return a shape name.
latlon2shape <- function(lat, lon,
projargs="+proj=longlat +datum=WGS84 +ellps=WGS84 +no_defs",
shape=c("county", "state")) {
shape = match.arg(shape)
# Prepare SpatialPolygons object with one SpatialPolygon
# per county
shapes <- map(shape, fill=TRUE, col="transparent", plot=FALSE)
IDs <- sapply(strsplit(shapes$names, ":"), function(x) x[1])
shapes_sp <- map2SpatialPolygons(shapes, IDs=IDs,
proj4string=CRS(projargs))
# Convert pointsDF to a SpatialPoints object
df_points = data.frame(x=lon, y=lat)
pointsSP <- SpatialPoints(df_points,
proj4string=CRS(projargs))
# Use 'over' to get _indices_ of the Polygons object containing each point
indices <- over(pointsSP, shapes_sp)
# Return the county names of the Polygons object containing each point
shapeNames <- sapply(shapes_sp@polygons, function(x) x@ID)
res = shapeNames[indices]
res = toupper(strsplit(res, ",")[[1]][c(2,1)])
res = list('county' = res[1], 'state' = res[2])
}
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment