Distance of Point Feature to Nearest Polygon in R

Distance of point feature to nearest polygon in R

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.

find the nearest polygon for a given point

If you're ok to convert to sf objects, you can find the nearest polygon to each point outside the polys using something on this lines:

set.seed(1)
library(raster)
library(rgdal)
library(rgeos)
library(sf)
library(mapview)
  
p <- shapefile(system.file("external/lux.shp", package="raster"))
p2 <- as(1.5*extent(p), "SpatialPolygons")
proj4string(p2) <- proj4string(p)
pts <- spsample(p2, n=10, type="random")
  
## Plot to visualize
plot(p, col=colorRampPalette(blues9)(12))
plot(pts, pch=16, cex=.5,col="red", add = TRUE)

# transform to sf objects
psf   <- sf::st_as_sf(pts) %>% 
    dplyr::mutate(ID_point = 1:dim(.)[1])
polsf <- sf::st_as_sf(p)

# remove points inside polygons
in_points  <- lengths(sf::st_within(psf,polsf))
out_points <- psf[in_points == 0, ]

# find nearest poly
nearest <- polsf[sf::st_nearest_feature(out_points, polsf) ,]  %>% 
    dplyr::mutate(id_point = out_points$ID)
nearest

> Simple feature collection with 6 features and 6 fields
> geometry type:  POLYGON
> dimension:      XY
> bbox:           xmin: 5.810482 ymin: 49.44781 xmax: 6.528252 ymax: 50.18162
> geographic CRS: WGS 84
>   ID_1       NAME_1 ID_2           NAME_2 AREA                       geometry id_point
> 1    2 Grevenmacher    6       Echternach  188 POLYGON ((6.385532 49.83703...        1
> 2    1     Diekirch    1         Clervaux  312 POLYGON ((6.026519 50.17767...        2
> 3    3   Luxembourg    9 Esch-sur-Alzette  251 POLYGON ((6.039474 49.44826...        5
> 4    2 Grevenmacher    7           Remich  129 POLYGON ((6.316665 49.62337...        6
> 5    3   Luxembourg    9 Esch-sur-Alzette  251 POLYGON ((6.039474 49.44826...        7
> 6    2 Grevenmacher    6       Echternach  188 POLYGON ((6.385532 49.83703...        9
> 

#visualize to check
mapview::mapview(polsf["NAME_2"]) + mapview::mapview(out_points)

HTH !

How to do spatial join of polygons to nearest point using polygon boundry?

Consider this piece of code; what it does is:

iterates over your polygons object, creating 1 row in results dataframe for each iteration
for a given polygon finds indices of nearest points; this will be, depending on the value of no_matches column, vector of length one to three
for each iteration find distance from your polygon to three points in the nearest vector; should the nearest vector have length less than 3 (as specified by no_matches column) NA will be returned.

Given the need to find more than 1 nearest object (which I kind of missed in my previous answer, now deleted) you will be better off with nngeo::st_nn() than sf::st_nearest_feature(). Using sf::st_distance() for calculating distance remains.

It is necessary to iterate (via a for cycle or some kind of apply) because the k argument of nngeo::st_nn() is not vectorized.

library(sf)
library(dplyr)

points <- st_read("./example/point.shp")
polygons <- st_read("./example/polygon.shp")

result <- data.frame() # initiate empty result set

for (i  in polygons$id) { # iterate over polygons
  
  # ids of nearest points as vector of lenght one to three
  nearest <- nngeo::st_nn(polygons[polygons$id == i, ], 
                          points, 
                          k = polygons[polygons$id == i, ]$no_matches,
                          progress = F)[[1]]
  
  # calculate result for current polygon
  current <- data.frame(id = i,
                        # id + distance to first nearest point
                        point_id_1 = points[nearest[1], ]$point_id,
                        distance_1 = st_distance(polygons[polygons$id == i, ],
                                                 points[nearest[1], ]),
                        # id + distance to second nearest point (or NA)
                        point_id_2 = points[nearest[2], ]$point_id,
                        distance_2 = st_distance(polygons[polygons$id == i, ],
                                                 points[nearest[2], ]),
                        # id + distance to third nearest point (or NA)
                        point_id_3 = points[nearest[3], ]$point_id,
                        distance_3 = st_distance(polygons[polygons$id == i, ],
                                                 points[nearest[3], ])
  )
  
  
  # add current result to "global" results data frame
  result <- result %>% 
    bind_rows(current)
  
}

# check results
result

Finding the nearest distances between a vector of points and a polygon in R

I think you get only one distance value because you overwrite g in every iteration of your for-loop. I do however not not know if this is the only problem because I cannot reproduce your issue without suitable data.
Try changing the last loop to this:

g = rep(NA, dim(spdf)[1])
for(i in 1:dim(spdf)[1]){
  g[i] = gDistance(spdf[i,],land_poly)
}

How do I find the polygon nearest to a point in R?

Here's an answer that uses an approach based on that described by mdsumner in this excellent answer from a few years back.

One important note (added as an EDIT on 2/8/2015): rgeos, which is here used to compute distances, expects that the geometries on which it operates will be projected in planar coordinates. For these example data, that means that they should be first transformed into UTM coordinates (or some other planar projection). If you make the mistake of leaving the data in their original lat-long coordinates, the computed distances will be incorrect, as they will have treated degrees of latitude and longitude as having equal lengths.

library(rgeos)

##  First project data into a planar coordinate system (here UTM zone 32)
utmStr <- "+proj=utm +zone=%d +datum=NAD83 +units=m +no_defs +ellps=GRS80"
crs <- CRS(sprintf(utmStr, 32))
pUTM <- spTransform(p, crs)
ptsUTM <- spTransform(pts, crs)

## Set up containers for results
n <- length(ptsUTM)
nearestCanton <- character(n)
distToNearestCanton <- numeric(n)

## For each point, find name of nearest polygon (in this case, Belgian cantons)
for (i in seq_along(nearestCanton)) {
    gDists <- gDistance(ptsUTM[i,], pUTM, byid=TRUE)
    nearestCanton[i] <- pUTM$NAME_2[which.min(gDists)]
    distToNearestCanton[i] <- min(gDists)
}

## Check that it worked
data.frame(nearestCanton, distToNearestCanton)
#       nearestCanton distToNearestCanton
# 1             Wiltz           15342.222
# 2        Echternach            7470.728
# 3            Remich           20520.800
# 4          Clervaux            6658.167
# 5        Echternach           22177.771
# 6          Clervaux           26388.388
# 7           Redange            8135.764
# 8            Remich            2199.394
# 9  Esch-sur-Alzette           11776.534
# 10           Remich           14998.204

plot(pts, pch=16, col="red")
text(pts, 1:10, pos=3)
plot(p, add=TRUE)
text(p, p$NAME_2, cex=0.7)

Sample Image

Find nearest distance from spatial point with direction specified

Thankyou to @patL and @Wimpel, I've used your suggestions to come up with a solution to this problem.

First I create spatial lines of set distance and bearings from an origin point using destPoint::geosphere. I then use gIntersection::rgeos to obtain the spatial points where each transect intersects the coastline. Finally I calculate the distance from the origin point to all intersect points for each transect line respectively using gDistance::rgeos and subset the minimum value i.e. the nearest intersect.

load packages

pkgs=c("sp","rgeos","geosphere","rgdal") # list packages
lapply(pkgs,require,character.only=T) # load packages

create data

coastline

coords<-structure(list(lon =c(-6.1468506,-3.7628174,-3.24646, 
-3.9605713,-4.4549561,-4.7955322,-4.553833,-5.9710693,-6.1468506), 
lat=c(53.884916,54.807017,53.46189,53.363665,53.507651,53.363665,53.126998,53.298056,53.884916)), class = "data.frame", row.names = c(NA,-9L))
l=Line(coords)
sl=SpatialLines(list(Lines(list(l),ID="a")),proj4string=CRS("+init=epsg:4326"))

point

sp=SpatialPoints(coords[5,]+0.02,proj4string=CRS("+init=epsg:4326"))
p=coordinates(sp) # needed for destPoint::geosphere

create transect lines

b=seq(0,315,45) # list bearings

tr=list() # container for transect lines

for(i in 1:length(b)){
    tr[[i]]<-SpatialLines(list(Lines(list(Line(list(rbind(p,destPoint(p=p,b=b[i],d=200000))))),ID="a")),proj4string=CRS("+init=epsg:4326")) # create spatial lines 200km to bearing i from origin
    }

calculate distances

minDistance=list() # container for distances

for(j in 1:length(tr)){ # for transect i
    intersects=gIntersection(sl,tr[[j]]) # intersect with coastline
    minDistance[[j]]=min(distGeo(sp,intersects)) # calculate distances and use minimum
}

do.call(rbind,minDistance)

In reality the origin point is a spatial point data frame and this process is looped multiple times for a number of sites. There are also a number of NULL results when carry out the intersect so the loop includes an if statement.

calculate distance from all parts of polygon to closest points

Here is my solution. I am using sf whenever possible. From my experience sf is not completely compatible with the raster functions yet, so there are a couple of workarounds here that aren't too ugly.

I am using different base data than what you provided.

Base Data

library(sf)
library(raster)
library(magrittr)

set.seed(1)

## We will create your polygons from points using a voronoi diagram
x <- runif(10, 640000, 641000)
y <- runif(10, 5200000, 5201000)

myPolyPoints <- data.frame(id = seq(x), x = x, y = y) %>%
    st_as_sf(coords = c("x", "y"))

## Creating the polygons here
myPolygons <- myPolyPoints$geometry %>%
    st_union %>%
    st_voronoi %>%
    st_collection_extract

myPolygons <- st_sf(data.frame(id = seq(x), geometry = myPolygons)) %>%
    st_intersection(y = st_convex_hull(st_union(myPolyPoints)))

## Creating points to query with buffers then calculate distances to
polygonExt <- extent(myPolygons)
x <- runif(50, polygonExt@xmin, polygonExt@xmax)
y <- runif(50, polygonExt@ymin, polygonExt@ymax)

myPoints <- data.frame(id = seq(x), x = x, y = y) %>%
    st_as_sf(coords = c("x", "y"))

## Set projection info
st_crs(myPoints) <- 26910
st_crs(myPolygons) <- 26910

## View base data
plot(myPolygons$geometry)
plot(myPoints$geometry, add = T, col = 'blue')

## write out data
saveRDS(list(myPolygons = myPolygons,
             myPoints = myPoints),
        "./basedata.rds")

The base data I generated looks like so:

View of base data

Distance Processing

library(sf)
library(raster)
library(magrittr)

## read in basedata
dat <- readRDS("./basedata.rds")

## makeing a grid of points at a resolution using the myPolygons extent
rast <- raster(extent(dat$myPolygons), resolution = 1, vals = 0, crs = st_crs(dat$myPoints))

## define a function that masks out the raster with each polygon, then
## generate a distance grid to each point with the masked raster
rastPolyInterDist <- function(maskPolygon, buffDist){
    maskPolygon <- st_sf(st_sfc(maskPolygon), crs = st_crs(dat$myPoints))
    mRas <- mask(rast, maskPolygon)
    cent <- st_centroid(maskPolygon)
    buff <- st_buffer(cent, buffDist)
    pSel <- st_intersection(dat$myPoints$geometry, buff)

    if(length(pSel) > 0){
        dRas <- distanceFromPoints(mRas, as(pSel, "Spatial"))
        return(dRas + mRas)
    }
    return(mRas)
}

dat$distRasts <- lapply(dat$myPolygons$geometry,
                        rastPolyInterDist,
                        buffDist = 100)

## merge all rasters back into a single raster
outRast <- dat$distRasts[[1]]

mergeFun <- function(mRast){
     outRast <<- merge(outRast, mRast)
}

lapply(dat$distRasts[2:length(dat$distRasts)], mergeFun)

## view output
plot(outRast)
plot(dat$myPoints$geometry, add = T)
dat$myPolygons$geometry %>%
    st_centroid %>%
    st_buffer(dist = 100) %>%
    plot(add = T)

The results can be seen below. You can see that there is a condition handled when the buffered centroid does not intersect any locations found in its polygon.

View of results

Using your base data, I have made the following edits to how your data is read and processed in R.

OP Base Data

library(raster)
library(sf)
library(magrittr)

url <- "https://www.dropbox.com/s/25n9c5avd92b0zu/example.zip?raw=1"
download.file(url, "example.zip")
unzip("example.zip")

myPolygons <- st_read("myPolygon.shp") %>%
    st_transform(st_crs("+proj=robin +datum=WGS84"))

myPoints <- st_read("myPoints.shp") %>%
    st_transform(st_crs("+proj=robin +datum=WGS84"))

centroids <- st_centroid(myPolygons)
buffer <- st_buffer(centroids, 5000)

plot(myPolygons, col="green")
plot(buffer, col="blue", add = T)
plot(centroids, pch = 20, col = "white", add = T)
plot(myPoints, pch = 20, col = "red", add = T)

saveRDS(list(myPoints = myPoints, myPolygons = myPolygons), "op_basedata.rds")

Distance Processing Using OP Data

To use the calculation routine I have suggested, you just need to modify the resolution of the starting raster, and the buffer distance input. Othewise it should behave the same once you read your data into R as I have outlined above.

library(sf)
library(raster)
library(magrittr)

## read in basedata
dat <- readRDS("./op_basedata.rds")

## makeing a grid of points at a resolution using the myPolygons extent
rast <- raster(extent(dat$myPolygons), resolution = 100, vals = 0, crs = st_crs(dat$myPoints))

## define a function that masks out the raster with each polygon, then
## generate a distance grid to each point with the masked raster
rastPolyInterDist <- function(maskPolygon, buffDist){
    maskPolygon <- st_sf(st_sfc(maskPolygon), crs = st_crs(dat$myPoints))
    mRas <- mask(rast, maskPolygon)
    cent <- st_centroid(maskPolygon)
    buff <- st_buffer(cent, buffDist)
    pSel <- st_intersection(dat$myPoints$geometry, buff)

    if(length(pSel) > 0){
        dRas <- distanceFromPoints(mRas, as(pSel, "Spatial"))
        return(dRas + mRas)
    }
    return(mRas)
}

dat$distRasts <- lapply(dat$myPolygons$geometry,
                        rastPolyInterDist,
                        buffDist = 5000)

## merge all rasters back into a single raster
outRast <- dat$distRasts[[1]]

mergeFun <- function(mRast){
     outRast <<- merge(outRast, mRast)
}

lapply(dat$distRasts[2:length(dat$distRasts)], mergeFun)

## view output
plot(outRast)
plot(dat$myPoints$geometry, add = T)
dat$myPolygons$geometry %>%
    st_centroid %>%
    st_buffer(dist = 5000) %>%
    plot(add = T)

View of OP results

Find nearest features using sf in R

It looks like the solution to my question was already posted -- https://gis.stackexchange.com/questions/243994/how-to-calculate-distance-from-point-to-linestring-in-r-using-sf-library-and-g -- this approach gets just what I need given an sf point feature 'a' and sf polygon feature 'b':

closest <- list()
for(i in seq_len(nrow(a))){
    closest[[i]] <- b[which.min(
    st_distance(b, a[i,])),]
}

Distance of Point Feature to Nearest Polygon in R