Distance of point feature to nearest polygon in R

user1738753 picture user1738753 · May 8, 2013 · Viewed 8.5k times · Source

I working on a project at the moment, where I have a point feature -- the point feature includes a 142 points -- and multiple polygon (around 10). I want to calculate the distance between every single point and the nearest polygon feature in R.

My current approach is tedious and a bit long winded. I am currently planning to calculate the distance between every single point and every single polygon. For example, I would calculate the distance between the 142 points and Polygon A, the distance between the 142 points and Polygon B, the distance between 142 points and Polygon C, etc. Here is a sample code of one of these distance calculations:

dist_cen_polya <- dist2Line(centroids_coor, polygonA_shp)

After doing these calculations, I would write a code to choose the minimum/nearest distance between every single point and the closest polygon. The issue is that this procedure is tedious.

Does anyone know a package/code which would minimize the effort/computational time of the calculation? I would really like to use a package that compare a single to point to the nearest polygon feature or calculates the distance between a point and all polygons of interest?

Thank you.

Answer

Jeffrey Evans picture Jeffrey Evans · May 9, 2013

Here I am using the gDistance function in the rgeos topology library. I am using a brute force double loop but it is surprisingly fast. It takes less than 2 seconds for 142 points and 10 polygons. I am sure that there is a more elgant way to perform the looping.

   require(rgeos)

    # CREATE SOME DATA USING meuse DATASET
    data(meuse)
      coordinates(meuse) <- ~x+y
        pts <- meuse[sample(1:dim(meuse)[1],142),]  
    data(meuse.grid) 
      coordinates(meuse.grid) = c("x", "y") 
        gridded(meuse.grid) <- TRUE 
          meuse.grid[["idist"]] = 1 - meuse.grid[["dist"]]    
        polys <- as(meuse.grid, "SpatialPolygonsDataFrame")
          polys <- polys[sample(1:dim(polys)[1],10),]   
    plot(polys)
      plot(pts,pch=19,cex=1.25,add=TRUE)      

    # LOOP USING gDistance, DISTANCES STORED IN LIST OBJECT
    Fdist <- list()
      for(i in 1:dim(pts)[1]) {
        pDist <- vector()
          for(j in 1:dim(polys)[1]) { 
            pDist <- append(pDist, gDistance(pts[i,],polys[j,])) 
          }
        Fdist[[i]] <- pDist
      } 

    # RETURN POLYGON (NUMBER) WITH THE SMALLEST DISTANCE FOR EACH POINT  
    ( min.dist <- unlist(lapply(Fdist, FUN=function(x) which(x == min(x))[1])) ) 

    # RETURN DISTANCE TO NEAREST POLYGON
    ( PolyDist <- unlist(lapply(Fdist, FUN=function(x) min(x)[1])) ) 

    # CREATE POLYGON-ID AND MINIMUM DISTANCE COLUMNS IN POINT FEATURE CLASS
    pts@data <- data.frame(pts@data, PolyID=min.dist, PDist=PolyDist)

    # PLOT RESULTS
    require(classInt)
    ( cuts <- classIntervals(pts@data$PDist, 10, style="quantile") )
       plotclr <- colorRampPalette(c("cyan", "yellow", "red"))( 20 )
         colcode <- findColours(cuts, plotclr)
    plot(polys,col="black")
      plot(pts, col=colcode, pch=19, add=TRUE)

The min.dist vector represents the row number of the polygon. For instance you could subset the nearest polygons by using this vector as such.

near.polys <- polys[unique(min.dist),]

The PolyDist vector contain the actual Cartesian minimum distances in the projection units of the features.