[R-sig-Geo] creating point sample of polygons using spsample
Tim Howard
tghoward at gw.dec.state.ny.us
Fri Mar 14 13:14:51 CET 2014
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]]
More information about the R-sig-Geo
mailing list