########################################################################## # PROGRAM: BreedingDensityAreas (BDA.R) # USE: CALCULATES BREEEDING DENSITY AREAS BASED ON POPULATION POINT COUNTS # REQUIRES: SP SpatialPointsDataFrame OBJECT WITH POP FIELD # sp, rgeos, spatstat packages # Version >= R 2.13.0 # # ARGUMENTS: # x SpatialPointsDataFrame OBJECT # pop POPULATION COUNT FIELD # dist BUFFER DISTANCE ( p < 0.75 dist = 6400m, p >= 0.75 dist = 8500m) # k.bw KDE BANDWIDTH DISTANCE # p TARGET PERCENT OF POPULATION # # REFERENCES: # Doherty, K.E., J.D. Tack, J.S. Evans, D.E. Naugle (2010) Mapping breeding densities of # greater sage-grouse: A tool for range-wide conservation planning. Bureau of Land # Management. Report Number L10PG00911. # # EXAMPLES: # require(rgdal) # require(sp) # setwd("D:/ANALYSIS/SageGrouse") # sdata = readOGR(dsn=getwd(), layer="leks") # bda75 <- BDA(sdata, pop="CountUSE", p=0.75, dist=8500) # plot(bda75) # points(sdata, pch=18) # writeOGR(bda75, dsn=getwd(), layer="BDA75", driver="ESRI Shapefile") # # UPDATE: # 01/20/2012 - CHANGED CODE TO DISSOLVE MULTIPART BUFFER RESULTS TO SINGLEPART AND CONVERT # TO SpatialPolygonsDataFrame. THIS ALLOWS FOR WRITING TO SHAPEFILE. # # CONTACT: # Jeffrey S. Evans, Ph.D. # Senior Landscape Ecologist # The Nature Conservancy # Central Science/DbyD # Affiliate Assistant Professor # Environment and Natural Resources # University of Wyoming # Laramie, WY 82070 # jeffrey_evans@tnc.org # (970) 672-6766 ########################################################################## BDA <- function(x, pop, p=0.75, dist=8500, k.bw=1000) { if (!require(sp)) stop("sp PACKAGE MISSING") if (!require(rgeos)) stop("rgeos PACKAGE MISSING") if (!require(spatstat)) stop("spatstat PACKAGE MISSING") if (!inherits(x, "SpatialPointsDataFrame")) stop("MUST BE SP SpatialPointsDataFrame OBJECT") pn=(1:length(names(x@data))) [names(x@data)==pop] pop.n <- sum(x@data[,pn]) * p sdata.ppp <- ppp(coordinates(x)[,1], coordinates(x)[,2], marks=x@data[,pn], window=ripras(coordinates(x))) win <- sdata.ppp$window kde <- density(sdata.ppp, weights=sdata.ppp$marks, sigma=k.bw, at="points") pop.den <- as.data.frame(cbind(X=sdata.ppp$x, Y=sdata.ppp$y, COUNT=sdata.ppp$marks, KDE=(round(kde*100000,digits=10)) )) pop.den <- pop.den[order(-pop.den$KDE),] i=0; j=0 results <- as.data.frame(array(0, dim=c( 0, 4 ))) while(i <= pop.n) { j=j+1 i = i + pop.den[,"COUNT"][j] results <- rbind(results, pop.den[j,]) } pts <- SpatialPointsDataFrame(coords=results[c("X","Y")], results) pop.buff <- gBuffer(pts, byid=FALSE, id=NULL, width=dist, joinStyle="ROUND", quadsegs=10) polys <- pop.buff@polygons[[1]]@Polygons pl <- vector("list", length(polys)) for (i in 1:length(polys)) { pl[i] <- Polygons(list(polys[[i]]), i) } pop.buff <- SpatialPolygons(pl) row.ids=sapply(slot(pop.buff, "polygons"), function(i) slot(i, "ID")) pop.buff <- SpatialPolygonsDataFrame(pop.buff, data.frame(FID=as.numeric(row.ids)) ) return( pop.buff ) }