Screen Points by Proximity

#This code removes points from a SpatialPointsDataFrame that do not meet a user-specified proximity
#rule, i.e., tosses out points that are within some specified distance of another point in the same
#SpatialPointsDataFrame.
 
#Note that this code is not necessarily optimized, meaning it may be possible to remove points in a
#different order that would result in fewer points having to be removed.
 
#Load required packages
require(sp)       #for handling spatial data types
require(rgeos)    #for the gDistance function
 
#Load example data for use as InputPoints
data(meuse)                     #Part of sp package
coordinates(meuse) <- ~x+y      #Make SpatialPointsDataframe
InputPoints <- meuse            #Save as InputPoints
head(InputPoints@data)          #View first few lines of data slot
plot(InputPoints, pch=1)        #Plot to visualize

#Specify the minimum distance points can be from eachother
MinDist <- 150 #unit is meters for projected coordinate systems
 
#Create a temp.id column used later in data.frame manipulation
InputPoints@data$temp.id<- 1:nrow(InputPoints)
 
#Create a complete distance matrix
d<- gDistance(InputPoints, byid=TRUE) #uses rgeos package
#Remove the 0s on the diagonal (distance of point to itself)
d[row(d) == col(d)] <- NA
 
#Loop through each point and determine which points are within MinDist
outDF<- data.frame()
for(i in 1:nrow(d)){
    d.match<- unique(which(d[i,]<MinDist))
    if(length(d.match)>0){
        outDF<- rbind(outDF, data.frame(inPt=i, closePt=paste(d.match, collapse=",")))
    }#close if statement
}#close for loop
 
#Loop through each close point and remove points as needed to produce OutputPoints
outDF$closePt<- as.character(outDF$closePt)     #Adjust data type in outDF
OutputPoints<- InputPoints                      #Copy InputPoints to OutputPoints
 
for(i in 1:nrow(outDF)){
    #First check to see if current outDF point still exists in the OutputPoints dataset
    if(!is.na(match(outDF$inPt[i], OutputPoints$temp.id))){
        #Point still exists in the dataset, so remove the other points that are too close
        kill<- which(OutputPoints@data$temp.id %in% unlist(strsplit(outDF[i,2], ",")))
        if(length(kill)>0){
            OutputPoints<- OutputPoints[-kill,]
        }#close inner if statement 
    }#close outer if statement
}#close for loop
 
#Plot to visualize
plot(InputPoints, col="dark gray", pch=1)                   #started with 155 points in example
plot(OutputPoints, col="green", pch=20, add=TRUE)           #83 retained in example
legend("topleft", cex = 1.2, bty = "n", legend = c("InputPoints", "OutputPoints"),
       col = c("gray", "green"), pch = c(1,20))

#Print summary of the fraction of InputPoints retained in OutputPoints
print(paste(nrow(OutputPoints), " (or", round((nrow(OutputPoints)/nrow(InputPoints))*100, 1),
            "%) of your original points were retained.", sep=""),quote=FALSE)
 
#Test to make sure no points are now too close in OutputPoints
g<- gDistance(OutputPoints, byid=TRUE)  #Uses rgeos package
g[row(g) == col(g)] <- NA               #Remove the 0s on the diagonal
min(g, na.rm=TRUE)                      #Print the minimum distance between points
 
#Remove temporary data created in the process
rm(d,g,outDF,MinDist,d.match,i,kill)