# How to associate point to nearest state boundary outside of the state the point is in

I have a point within a state and I want to find 1) the distance to the nearest state border and 2) the bordering state's name. I am able to do 1). For 2). I am able to get the state name corresponding to the nearest boundary, but sometimes this state is the state the point is in, and not actually the bordering state's name.

Modifying and extending the code from an excellent answer of this post (Calculate distance between points and nearest polygon in R), I do:

library(sp)
library(spData)
library(tidyverse)

states <- as(us_states,"Spatial")

#create some points
pts <- data.frame(x1=c(-100.5, -98.6, -98), x2=c(35, 41, 44))

# compute the shortest distance between points and polygons
# (from ?dist2Line): "returns matrix with distance and lon/lat of the nearest point" &
# "the ID (index) of (one of) the nearest objects"; distance is in meters (default)
dist.mat <- geosphere::dist2Line(p = pts, line = states)

# bind results with original points
pts.wit.dist <- cbind(pts, dist.mat)
pts.wit.dist[1:3,]

pts.sp <- sp::SpatialPoints(coords      = pts[,c("x1","x2")], # order matters
proj4string = states@proj4string)
plot(pts.sp, col="red")

# plot arrows to indicate the direction of the great-circle-distance
for (i in 1:nrow(pts.wit.dist)) {
arrows(x0 = pts.wit.dist[i,1],
y0 = pts.wit.dist[i,2],
x1 = pts.wit.dist[i,4],
y1 = pts.wit.dist[i,5],
length = 0.1,
col = "green")
}

a <- pts.wit.dist
a <- a %>% mutate(miles = distance*0.000621)

#create crosswalk between ID (from dist2Line) and state name (from state shapefile)
state_names <- data.frame(us_states\$NAME) %>% mutate(ID=row_number())
names(state_names)[1] <- "state"

b <- left_join(pts.wit.dist, state_names, by="ID")

b

This produces:

> b
x1 x2  distance        lon      lat ID        state
1 -100.5 35  45558.45 -100.00039 35.00162 23        Texas
2  -98.6 41 110992.87  -98.59800 40.00294 37     Nebraska
3  -98.0 44 118342.74  -98.49855 42.99856 22 South Dakota

However, the state is incorrect in all cases. It should be:

> b
x1 x2  distance        lon      lat ID        state
1 -100.5 35  45558.45 -100.00039 35.00162 23        Oklahoma
2  -98.6 41 110992.87  -98.59800 40.00294 37     Kansas
3  -98.0 44 118342.74  -98.49855 42.99856 22 Nebraska

It's tricky because the borders obviously are borders between two states, and I am not sure how to make it give me the state that is not the one in which the point is contained.

How can I do this?

• I think you have to drop the state the point is in from the search set - imagine a point near the california coast - without dropping california from the search set the shortest distance will be to the coast, and not (say) to a point on the CA-NV border so you could identify it as having NV as the nearest other state. Commented Aug 13, 2022 at 7:03

From my point of view, your problem can be seen as a nearest neighbour search using points and polygons, considering the condition to exclude the state your current point is located. My take making use of sf and nngeo:

library(sf)
library(nngeo)

# create sf objects with same crs -> WGS 84
states <- us_states |> st_transform(st_crs(4326))

pts <- data.frame(x = c(-100.5, -98.6, -98),
y = c(35, 41, 44)) |> st_as_sf(coords = c("x", "y"), crs = 4326)

# init empty list
results <- list(states = list(),
distances = list())

for (i in 1:dim(pts)[1]) {

# select current point
pt <- pts[i, ]
state_pt <- states[pt, ]

# which neighbours to consider? exlude state containing pt first
n_states <- dplyr::filter(states, NAME != state_pt[["NAME"]])

# determine nn
result <- nngeo::st_nn(pt, n_states, k = 1, returnDist = TRUE)

# return nn state and corresponding distance to the border
results[["states"]][[i]] <- n_states[result[["nn"]] |> unlist(), ]
results[["distances"]][[i]] <- result[["dist"]] |> unlist()
}

results
#> \$states
#> \$states[[1]]
#> Simple feature collection with 1 feature and 6 fields
#> Geometry type: MULTIPOLYGON
#> Dimension:     XY
#> Bounding box:  xmin: -103.0024 ymin: 33.63787 xmax: -94.43122 ymax: 37.0001
#> Geodetic CRS:  WGS 84
#>    GEOID     NAME REGION            AREA total_pop_10 total_pop_15
#> 19    40 Oklahoma  South 180971.4 [km^2]      3675339      3849733
#>                          geometry
#> 19 MULTIPOLYGON (((-103.0022 3...
#>
#> \$states[[2]]
#> Simple feature collection with 1 feature and 6 fields
#> Geometry type: MULTIPOLYGON
#> Dimension:     XY
#> Bounding box:  xmin: -102.0517 ymin: 36.99308 xmax: -94.60703 ymax: 40.00308
#> Geodetic CRS:  WGS 84
#>   GEOID   NAME  REGION            AREA total_pop_10 total_pop_15
#> 9    20 Kansas Midwest 213037.1 [km^2]      2809329      2892987
#>                         geometry
#> 9 MULTIPOLYGON (((-102.0517 4...
#>
#> \$states[[3]]
#> Simple feature collection with 1 feature and 6 fields
#> Geometry type: MULTIPOLYGON
#> Dimension:     XY
#> Bounding box:  xmin: -104.0532 ymin: 40 xmax: -95.30829 ymax: 43.00078
#> Geodetic CRS:  WGS 84
#>    GEOID     NAME  REGION            AREA total_pop_10 total_pop_15
#> 36    31 Nebraska Midwest 200272.3 [km^2]      1799125      1869365
#>                          geometry
#> 36 MULTIPOLYGON (((-104.0531 4...
#>
#>
#> \$distances
#> \$distances[[1]]
#> [1] 45507.54
#>
#> \$distances[[2]]
#> [1] 110868.8
#>
#> \$distances[[3]]
#> [1] 118392.9

Et voilà - Oklahoma, Kansas, Nebraska with distances attached should be your desired output.

• Thanks @falk-env! Is there a way to also return the distance between the point and the nearest neighbor state border? Commented Aug 13, 2022 at 15:42
• Yeah, there actually is making use of nngeo::st_nn(..., returnDist = TRUE). I edited my answer to consider this. By the way, distances are given in meters. Commented Aug 13, 2022 at 16:50
• Thank you! This is great. Commented Aug 13, 2022 at 16:51