[R-sig-Geo] creating point sample of polygons using spsample

[Tim Howard]

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

> Date: Thu, 13 Mar 2014 20:07:09 -0400
> From: Spencer Meyer <spencer.meyer at maine.edu>
> To: r-sig-geo at r-project.org 
> >Subject: [R-sig-Geo] creating point sample of polygons using spsample
> Message-ID: <ddc1888cd851e4e504906a4624fd9648 at mail.gmail.com>
> Content-Type: text/plain
> 
> Greetings - I am using R 2.15.2 on Windows 7 and am seeking help with the
> spsample function in R. I have a polygon dataset (n= 5733) and I'm trying
> to use the spsample function to sample the polygons to create a point
> dataset. I'm using sapply to create the sample based on code I found here:
> 
> http://casoilresource.lawr.ucdavis.edu/drupal/node/644 
> 
> 
> 
> Problem #1
> 
> The polygons are vastly different sizes (1-309,000 hectares) so I'd like to
> vary the sample size (n) based on the size of the polygon. I've tried
> passing a vector of the size of each polygon to the function, but I don't
> believe I'm linking the vector to the index (i) properly. I have also tried
> accessing the slot "area" of the x at polygons slot but I cannot seem to get a
> slot of a slot. I have also tried using the "stratified" type and adjusting
> the cellsize but that does not appear to change the sample size differently
> for different polygons.
> 
> 
> 
> Problem #2
> 
> I am having trouble joining the attributes from my original polygons back
> to my resulting points. I have used the sample from the website listed
> above to attach the polygon ID to each point but I cannot also get the
> attributes. The method for extracting the ID relies on accessing the
> polygons slot but then I have to access the data slot to get the attributes
> and I can't seem to do both.
> 
> 
> 
> Here is my code so far. As you can see I'm an R rookie, particularly when
> it comes to understanding how different data classes relate to each other.
> 
> Thank you for your help!
> 
> _____________________________________________________________________
> 
> 
> 
> #### 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]]