我正在尝试编写一个随机分配例程的脚本。采样设计将一个区域划分为许多多边形或地层。一组但数量不同的样本将随机分配给每个地层(至少2个样本，但在某些地层中多达7个)。因此，我有一个地层的形状文件，在它的"属性表"中，每个地层的名称和所需样本的数量。
- 地层；样品
- 440；4
- 441；2
- 5Z3；4
- 5Z1；7
- 560；2
我发现了一些关于这些类型的采样设计的好文档(http://casoilresource.lawr.ucdavis.edu/drupal/book/export/html/519)尽管我在为自己的需要执行例程时遇到了一些问题。在这个链接之后，我一直在使用rgdal和maptools包。我的工作脚本如下：
# 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')
然而，每个地层的采样强度并不是静态的(其中n=4)。我需要它来反映属性表中的列，该列指示给定地层所需的相应样本数。我还想将地层名称重新分配给已分配的采样位置。
理想情况下，例程将遍历每个多边形，并在其中随机分配n个样本(如属性表中所示)，并作为shapefile编写。任何指示都将不胜感激。我对这个节目还比较陌生，所以如果我错过了一些明显的东西，我很抱歉。
[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
很遗憾，如果不打开shapefile，我无法帮助您进行适当的采样。但让我给你一些建议，也许会有所帮助：
- 如果点数不同，并且不是像您提到的那样是静态的因此，使用for循环将是明智的，这将有助于区分每次的点的数量是分布在特定的多边形区域上
- SpatialPolygons-对象如此受约束，以至于您的第一个多边形将在此处的">x@polygons[[i]]@polygons[[1]]"下找到i=1是第一个多边形，i=2为第二个多边形，依此类推i=n是您的第n个多边形，x则是您的SpatialPolygonsDataFrame
希望我能帮助你更进一步。