I have several hundred points that represent samples in a landscape:

enter image description here

As you can see, the points are unevenly spread and form large, dense clusters. To account for spatial autocorrelation, I want to subset or thin these points so that the points that remain are a minimum distance from each other. At the moment I've set up a pairwise distance matrix for them and am randomly removing points that are below the distance threshold until the only ones remaining meet my minimum distance criteria. However, this does not seem optimal from two different perspectives: first, different potential random subsets using this method can return different numbers of samples. Ideally, I'd retain as many points as possible. Second, different subsets may be distributed differently across the landscape. Ideally, I want a point subset with a distribution that more uniformly represents the landscape. Are there any recommended or general algorithms for optimally subsetting points (at least more optimal than my approach)?


There is a function pp.subsample in the spatialEco package that will subset your data using a spatial intensity function. This will, in effect, thin your data to a defined n, using the expected density as sample weights.

However, in looking at the spatial distribution of your data, unless you subsample to a handful of observations, this will not mitigate the spatial autocorrelation (pseudoreplication) issues in this particular data. You have a notable sample distribution problem.

Here is a worked example, from the functions help, that illustrates how this works.

First, lets add the required packages and some data to work with.


  trees <- as(bei, 'SpatialPoints')
  plot(trees, pch=20, col='black', cex=0.50)
    title("Observed population")

observed population

Now we can set a 10% subsample target, which is fairly arbitrary and one should define n based on their data and analysis goals. We then run the pp.subsample function to use the intensity estimate as sample weights in pulling the subsample.

n = round(length(trees) * 0.10, digits=0)  
trees.wrs <- pp.subsample(trees, n=n, window='hull')   

Here are the points plotted on the estimated intensity function. Please note that the sigma parameter is quite important. I am defaulting to the first-order Scott bandwidth estimate. This smooths the point process quite a bit. If you want more of a second-order process, either specify your own sigma value, or use a second-order sigma estimate such as diggle or stoyan.

intensity function

Plotting the results illustrates the expected spatially weighted sample distribution and a notable thinning of the data.

plot(trees, pch=20, col='black', cex=0.50)
  plot(trees.wrs, pch=20, col='red', cex=1.25, add=TRUE) 
    title('10% subsample (n=360) of population (n=3604)')
    legend('bottomright', legend=c('Original sample', 'Subsample'), 


| improve this answer | |
  • Thank you for the detailed answer. It's certainly better than the approach I had when I wrote the question, and better than anything else I was able to find online. I ended up figuring out my own algorithm that has the specific properties I'm looking for, and will also post it as an answer when I publish my work using it. – anjama Dec 20 '17 at 17:46

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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