I've completed a landcover classification on an Landsat Image and would like to start an accuracy assessment. I'd like a stratified random sample that has a sample number that varies based on each landcover values total image pixels. In other words a proportional stratified random sample.

My code so far is the following

spdf = sampleStratified(FinalClassImg, 200, sp = T)

This generates 200 points per landcover class but what I really need is a variable number of samples based on the overall landcover pixel totals.

I roughly calculated the percentages in the image and decided to use 5000 total points. I multiplied the percentages * 5000 and got the following numbers for the nine landcover classes

100, 439, 643, 172, 194, 45, 39, 27, 3341

These numbers represent the number of random samples I'd like selected within each landcover class.

I would like to feed this set of numbers into the above code in place of the 200

I tried this:

spdf = sampleStratified(FinalClassImg, c(100, 439, 643, 172, 194, 45, 39, 27, 3341), sp = T)

but I get a mess of warnings

1: In if (nrow(y) < size) { ... : the condition has length > 1 and only the first element will be used - fewer samples than requested found for strata: 2, 3, 4, 5, 9

and only the first number is used (100 samples per class).

Is there a way to do this or does the sampleStratified only work with a single number?

1 Answer 1


?sampleStratified says:

size: positive integer giving the number of items to choose

not a vector of integers.

You could do what you want by sampling over each category. This might not be very efficient and there's probably a function somewhere that does all this.

First a function that samples N from category C by stratified sampling of the r==C raster and keeping only the samples in the part where r==C:

> sampleNfromC = function(r,N,C){

So to get 2 samples from category 2 do:

> sampleNfromC(r,2,2)
  cell layer
3   37     2
4   70     2

Then loop this over your categories and counts.

> cats = c(2,3,4)
> counts=c(3,5,2)
> strat = do.call(rbind, lapply(1:3, function(i){sampleNfromC(r,counts[i],cats[i])}))
> strat
   cell layer
4    46     2
5    38     2
6     5     2
61   10     3
7     2     3
8     7     3
9    96     3
10   44     3
3    66     4
41   61     4
  • Very helpful reply! I really need the results to be in a spatial points dataframe so that I can pull the same info from the Accuracy Assessment image I have. The results of strat are a dataframe.
    – jport
    Jan 6, 2020 at 14:39
  • 1
    Use strat=cbind(strat,xyFromCell(r,strat$cell)) to get the X-Y coordinates for each cell and make a SpatialPointsDataFrame with coordinates(strat)=~x+y. Then set the CRS. Or modify my function a bit.
    – Spacedman
    Jan 6, 2020 at 16:02

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

Not the answer you're looking for? Browse other questions tagged or ask your own question.