Quantifying spatial data (e.g. land cover) around points can be done in a variety of ways, some of which require considerable amounts of patience, clicking around, and/or cash for a license. Here’s a bit of code that I cobbled together to quickly extract land cover data from the National Land Cover Database for buffered regions around points (e.g. small ponds, point-count locations, etc.), in this case, U.S. capitals.

Here, extract_cover() is a function to do the extraction (with the help of the raster package), and extraction.R makes a parallel call to the function using the doMC package:

# Function to extract the data in R instead of Arc
# inputs: year, buffer size (in meters), point data file,
#   cover data directory, and output file directory
# output: csv file with site id, cover type, and % in buffer

extract_cover <- function(year, buffer,
                          point_d = "sites.csv",

  # load cover data
  filename <- paste(data_dir, "/nlcd_",
  NLCD <- raster(filename)

  # load site data
  sites <- read.csv(point_d, header=T)  
  coords <- sites[, c("longitude", "latitude")]

  #convert lat/lon to appropriate projection
  names(coords) <- c("x", "y")
  coordinates(coords) <- ~x + y
  proj4string(coords) <- CRS("+proj=longlat +ellps=WGS84 +datum=WGS84")
  crs_args <- NLCD@crs@projargs
  sites_transformed <- spTransform(coords, CRS(crs_args))

  #extract land cover data for each point, given buffer size
  Landcover <- extract(NLCD, sites_transformed, buffer=buffer)

  # summarize each site's data by proportion of each cover type
  summ <- lapply(Landcover, function(x){

  # generate land cover number to name conversions
  num.codes <- unique(unlist(Landcover))
  cover.names <- NLCD@data@attributes[[1]]$Land.Cover.Class[num.codes + 1]
  levels(cover.names)[1] <- NA # first level is ""
  conversions <- data.frame(num.codes, cover.names)
  conversions <- na.omit(conversions)
  conversions <- conversions[order(conversions$num.codes),]

  # convert to data frame
  mydf <- data.frame(id = rep(sites$id, lapply(summ, length)),
                     cover = names(unlist(summ)),
                     percent = unlist(summ)

  # create cover name column
  mydf$cover2 <- mydf$cover
  levels(mydf$cover2) <- conversions$cover.names

  # write output
  out_name <- paste(write_dir, "/",
                    year, "_cover_", buffer,
                    "_m_buffer.csv", sep="")
  write.csv(mydf, out_name)
# Script to actually extract the data at various buffer distances
# parallelized for speed (one core per year)

years <- c(2001, 2006, 2011)
nyears <- length(years)

# input vector of distances (in meters)
buffer_distances <- c(1000)

foreach (i=1:nyears) %dopar% {
  for (j in buffer_distances){
    extract_cover(year = years[i], buffer=j)