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

Question

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(states, add=TRUE)

# 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?

Answer 1

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.