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@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@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@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@proj4string <- poly.new@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.me...@maine.edu [[alternative HTML version deleted]] _______________________________________________ R-sig-Geo mailing list R-sig-Geo@r-project.org https://stat.ethz.ch/mailman/listinfo/r-sig-geo
