1
votes

I am trying to script a random allocation routine. The sampling design has a region divided into numerous polygons or strata. A set, but different number of samples are to be allocated randomly to each strata (a minimum of 2 samples but as many as 7 in some strata). As such, I have a shapefile of the strata and in its' attribute table, the strata names and number of samples desired in each.

  • STRATA; SAMPLES
  • 440; 4
  • 441; 2
  • 5Z3; 4
  • 5Z1; 7
  • 560; 2

I found some good documentation on these types of sampling designs (http://casoilresource.lawr.ucdavis.edu/drupal/book/export/html/519) although I'm having some issues implementing the routine for my own needs. Following this link, I have been working with both rgdal and maptools packages. My working script is, as follows:

# read in strata boundary shapefile 
strataboundaries <- readOGR('strataboundaries.shp', layer='strataboundaries')
#sample allocation to strata
allocation <- sapply(slot(strataboundaries, 'polygons'), function(i) spsample(i, n=4, type='random'))
allocation.merge <- do.call('rbind', allocation)
stratumID <- sapply(slot(strataboundaries, 'polygons'), function(i) slot(i, 'ID'))
sample <- sapply(allocation, function(i) nrow(i@coords))
sampleID <- rep(stratumID, sample)
allocation.final <- SpatialPointsDataFrame(allocation.merge, 
data=data.frame(poly_id=sampleID))
plot(strataboundaries, col="lightcyan", bborder="black", axes=TRUE, bg="lightsteelblue1")
points(allocation.final, col="red", pch=3, cex=0.8)
#write out shapefile containing sampling locations
allocation.final@proj4string <- strataboundaries@proj4string
writeOGR(allocation.final, ".", "allocation", driver='ESRI Shapefile')

However, the sampling intensity per stratum is not static (where I have n=4). I need this to reflect the column in the attributes table, which indicates the respective number of samples required for a given strata. I would also like to assign the strata names back to the sampling locations that have been allocated.

Ideally, the routine would iterate though each polygon and assign n number of samples randomly within (as indicated in the attribute table) and be written as a shapefile. Any direction would be greatly appreciated. I'm relatively new to the program so I apologize if I've missed something obvious.

[1] "SpatialPolygonsDataFrame"
attr(,"package")
[1] "sp"

Formal class 'SpatialPolygonsDataFrame' [package "sp"] with 5 slots
  ..@ data       :'data.frame': 66 obs. of  3 variables:
  .. ..$ OBJECTID_1: int [1:66] 1 2 3 4 5 6 7 8 9 10 ...
  .. ..$ Stratum1  : Factor w/ 66 levels "440","441","442",..: 65 64 63 62 61 60 12 11 7 49 ...
  .. ..$ Primary   : int [1:66] 2 2 4 5 2 7 2 2 2 2 ...
  ..@ polygons   :List of 66
  .. ..$ :Formal class 'Polygons' [package "sp"] with 5 slots
  .. .. .. ..@ Polygons :List of 1
  .. .. .. .. ..$ :Formal class 'Polygon' [package "sp"] with 5 slots
  .. .. .. .. .. .. ..@ labpt  : num [1:2] -68.3 40.4
  .. .. .. .. .. .. ..@ area   : num 0.769
  .. .. .. .. .. .. ..@ hole   : logi FALSE
  .. .. .. .. .. .. ..@ ringDir: int 1
  .. .. .. .. .. .. ..@ coords : num [1:654, 1:2] -66.3 -66.3 -66.4 -66.4 -66.4 ...
  .. .. .. ..@ plotOrder: int 1
  .. .. .. ..@ labpt    : num [1:2] -68.3 40.4
  .. .. .. ..@ ID       : chr "0"
  .. .. .. ..@ area     : num 0.769
1

1 Answers

0
votes

Unfortunately I can't help you to carry out an appropriate sampling without opening your shapefile. But let me give you some advices, which maybe can be helpful:

  1. If the number of points vary and is not static like you mentioned above so it will be cleverly to use a for-loop, which will help to distingush the number of points each time is to be distributed over particular polygonal area.
  2. SpatialPolygons - Objects are so constracted that your first polygon will be found under "x@polygons[[i]]@Polygons[[1]]" here i=1 is first polygon, i=2 is the second polygon and so on at the end i=n is your n-th polygon, x is your SpatialPolygonsDataFrame.

Hopefully I could help you to make a step further.