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.meyer@maine.edu

	[[alternative HTML version deleted]]

