4

I need to randomly generate ~ 40K points throughout my region of interest, but I need to weight the point selection by a set of probability weights, stored as a raster - such that higher probability areas have more of the random points and lower probability areas have fewer points. Additionally, I need to indicate a minimum distance between the points. I found that I could do this in ArcMap using the "create spatially balanced points" tool in the GeoStatistical Analyst package - but I do not have a license for that package. Is there a way to do this in QGIS or R? I found the sp package in R to randomly select points - but I do not see how I can weight that selection given the "type" options.

Any suggestions on how to solve this problem? Other R packages that I should check out?

Improve this question
1
  • 1
    A minimum distance constraint can interfere so strongly with the probability weights that sometimes a solution is impossible or else it is quite non-random. Could you explain the purpose of this random point generation exercise?
    – whuber
    Commented Jan 13, 2017 at 18:04

1 Answer 1

Reset to default
4

I don't understand the "minimum distance" thing you mentioned.

Here's how to generate 1000 points uniformly within cells but with the number in each cell weighted by the cell value:

Make a test 3x4 raster with some positive random numbers:

> set.seed(12)
> r = raster(matrix(runif(12),3,4))

Get the cell half-width for later:

> hs = res(r)/2

Now work out which cell each of our 1000 points is going in by sampling from the number of cells (12) with replacement, weighted by the value in the cells:

> ptscell = sample(1:12, 1000, prob=r[], replace=TRUE)

Now find the centre of those 1000 cell numbers:

> centres = xyFromCell(r,ptscell)

And generate random uniform points in the cell by using the centre and the half-width/height from earlier:

> pts = cbind(runif(nrow(centres),centres[,1]-hs[1],centres[,1]+hs[1]),runif(nrow(centres),centres[,2]-hs[2],centres[,2]+hs[2]))

Voila:

> plot(r)
> points(pts)

points weighted by raster value

Improve this answer
6
  • Thank you! 2 questions: 1. I am running into an error with this line (when I use my own raster): ptscell = sample(1:417, 1000, prob=r[], replace=TRUE) Error in sample.int(length(x), size, replace, prob) : incorrect number of probabilities Is the "length(x) the number of of grid cells in the ratser? I used zonal stats in qGIS to get this count and and I still get this number, but I still get the error1. 2. Concerning the minimum distance - I would like each point to be distance by at least 100 m.
    – user65148
    Commented Jan 13, 2017 at 19:43
  • That "417" has to be the number of cells in the raster. So its the number of rows times the number of columns. 417 can only be 3x139, which seems unlikely. Use 1:ncell(r).
    – Spacedman
    Commented Jan 13, 2017 at 19:47
  • Generating N points in a cell such that no two points are closer than X is possible via certain inhibition process functions in the spatstat package. But note its impossible if NpiX^2 > the cell area. Also, I'm not sure how you'd keep points that distance away at cell boundaries, unless you generate the points within a reduced cell so there's an empty guard area at the cell edge.
    – Spacedman
    Commented Jan 13, 2017 at 19:51
  • Thanks for the tip of the number of cells. I also just realized that you can get this # from R "dimensions : 2160, 4320, 9331200 (nrow, ncol, ncell)", so I changed my code, and got the same error: ptscell = sample(1:9331200, 1000, prob=r[], replace=TRUE) That solved that problem.
    – user65148
    Commented Jan 13, 2017 at 20:28
  • But now I am getting a new error, due to that I indicated that there were NA's using: r[r==-3.402823e+38] <- NA ptscell = sample(1:9331200, 1000, prob=r[], replace=TRUE) Error in sample.int(length(x), size, replace, prob) : NA in probability vector Any thoughts?
    – user65148
    Commented Jan 13, 2017 at 20:31

Your Answer

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

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