creating point sample of polygons using spsample - R-SIG-Geo

Fri, Mar 14, 2014 5:14 AM #

Spencer,
If you have the energy to try a different approach, I've had good success using the spsurvey package for this kind of thing. It generates spatially-balanced random points and can accommodate what you need. The setup is kind of complex, but there are some pretty good vignettes, including this one:
 
http://cran.at.r-project.org/web/packages/spsurvey/vignettes/Area_Design.pdf 
 
You want the stratified, equal probability, GRTS survey design (begins on page 9). 

Cheers, 
Tim

_____________________________________________________________________



#### Create Sample Points in R
------------------------------------------------------

library(rgdal)



## check polygons for holes

pls <- slot(poly.proj, "polygons")   #extract polygons from
spatialpolygondataframe

pls_new <- lapply(pls, checkPolygonsHoles)  # check for holes

poly.new <- SpatialPolygonsDataFrame(SpatialPolygons(pls_new,

            proj4string=CRS(proj4string(poly.proj))), data=as(poly.proj,
"data.frame"))  # fix holes



## create point sample within polygons

## followed http://casoilresource.lawr.ucdavis.edu/drupal/node/644 



# try three different approaches

polyPts = sapply(slot(poly.new, 'polygons'), function(i) spsample(i, n=100,
type='hexagonal', offset=c(0,0)))  # works but doesn't vary sample by
polygon size



polyPts2 = sapply(slot(poly.new, 'polygons'), function(i) spsample(i,
n=as.integer((i at area/10000)), type='hexagonal', offset=c(0,0)))  # works
but creates weird data structure that can't be used in merge below



polyPts3 = sapply(slot(poly.new, 'polygons'), function(i) spsample(i,
n=myn[i], type='hexagonal', offset=c(0,0)))  # should work in theory but
gets S4 index error.



# we now have a list of SpatialPoints objects, one entry per polygon of
original data

# stack into a single SpatialPoints object

plot(poly.new)

polyPts.merged <- do.call('rbind', polyPts) # this works when n above is
constant but not when varied

polyPts2.merged <- do.call('rbind', polyPts2)  # can't rbind the data. gets
error about proj4string being null even though poly.new has it set

polyPts2.merged2 = ldply(polyPts2, rbind) # doesn't work - mixing object
classes, I think.

polyPts2.merged3 = ldply(polyPts2, data.frame) # doesn't work either
because the result should be spatial*dataframes



length(polyPts.merged)  # check to see if the sample points actually varied
between the two methods

length(polyPts.merged2)



points(polyPts.merged, col='red', pch=3, cex=.25)  # plot results

points(polyPts.merged2, col='blue', pch=3, cex=.25)



# add an id, based on source polygon:

ids <- sapply(slot(poly.new, 'polygons'), function(i) slot(i, 'ID')) #
extract the original IDs

# determine the number of ACTUAL sample points generated for each poly

npts <- sapply(polyPts, function(i) nrow(i at coords))

# generate a reconstituted vector of point IDs

pt_id <- rep(ids, npts)

# promote to SpatialPointsDataFrame

polyPts.final <- SpatialPointsDataFrame(polyPts.merged,
data=data.frame(poly_id=pt_id))

# check:

plot(poly.new)

points(polyPts.final, col=polyPts.final$poly_id, pch=3, cex=0.5)

polyPts.final at proj4string <- poly.new at proj4string  # reassing projection



# get original attributes

test = over(polyPts.final, poly.new)

test2 = cbind(poly.new, test)

test3 = aggregate(polyPts.final, poly.new)

test2 = spCbind(polyPts.final,poly.new)



# another partial approach to recoving attribute data - doesn't work.

o <- match(polyPts.final$pt_id, poly.new$poly_id)

xtra1 <- poly.new[o,]

row.names(xtra1) <- xx$FIPSNO

xx1 <- spCbind(xx, xtra1)

names(xx1)

identical(xx1$CNTY_ID, xx1$CNTY_ID.1)

_____________________________________________________________________





Spencer R. Meyer
Center for Research on Sustainable Forests

University of Maine
spencer.meyer at maine.edu 

[[alternative HTML version deleted]]