Calculating the minimal distance to the nearest polygon for multipolygons in different crs

Question

I'm currently working on 31 different multipolygon shapefiles for which I need to generate a minimal distance to the nearest polygon "ANPs" (which is a multipolygon itself). I haven't been able to do it easily. I have had to modify the crs of the "ANPs" shapefile and disband the multipolygons in order to obtain individual polygons and then finding the closest neighbor and calculate minimal distance for it.

In the picture, you can see some of the 31 polygon shapefiles and the centroids of the "ANPs" shapefile. For each district, I want to generate the minimal distance for the closest blue dot.

I haven't been able to come up with an iterative/loop function to do this code automatically, but here is what I have came up with for:

`st_crs(ANPs)    #returns, "MEXICO_ITRF_2008"
 st_crs(DISTAgs) #returns, "WGS 84 / UTM zone 13N"

 ANPs_84=st_transform(ANPs, 32613)  #transforms to WGS 84 to UTM zone 13 N.

 ANPs= st_cast(ANPs_84,  "POLYGON") #disbanding to polygons. 
 DISTAgs=st_cast(DISTAgs, "POLYGON")

 ##disbanding ANPs results in the following error message.
 ##In st_cast.sf(ANPs, "POLYGON") :
 ## repeating attributes for all sub-geometries for which they may not be constant

 ANPsc=st_centroid(ANPs)            #calculating centroids.
 DISTAgsc=st_centroid(DISTAgs)

 st_nearest_feature(DISTAgsc, ANPsc) #nearest feature`

And that's as far as I got.

So, Should I worry about the multiple polygons generated by st_cast? (I think not because they are labeled). Is there an interative way to do this for all the other DISTx shapefiles? Will this be to innacurate?

DISCLAIMER: I'm new to cartography, GIS and more or less to RStudio.

Answer 1

Here is an example of calculating the distance to nearest, from one set of polygons to another set of polygons (or multi-polygons). The st_nn function uses st_distance internally, but it has parallel processing enabled so that it can be faster depending on the number of cores used (five, in this example). For 250/250 circular polygons, this takes 17 seconds.

library(nngeo)

# Sample data
set.seed(1)
n = 500
x = data.frame(lon = (runif(n) * 2 - 1) * 70, lat = (runif(n) * 2 - 1) * 70)
x = st_as_sf(x, coords = c("lon", "lat"), crs = 4326)
x = st_transform(x, 32630)
x = st_buffer(x, 100000)
x1 = x[c(TRUE, FALSE), ]
x2 = x[c(FALSE, TRUE), ]

# Plot
plot(x1, border = "blue")
plot(x2, border = "red", add = TRUE)

# Find distance to nearest - Parallel processing
start = Sys.time()
result = st_nn(x1, x2, k = 1, parallel = 5, returnDist = TRUE)
end = Sys.time()
end - start

## Time difference of 17.31085 secs

# Result
x1$dist_to_nearest = sapply(result$dist, "[", 1)
x1

## Simple feature collection with 250 features and 1 field
## geometry type:  POLYGON
## dimension:      XY
## bbox:           xmin: -7614809 ymin: -9040221 xmax: 12380790 ymax: 9045934
## projected CRS:  WGS 84 / UTM zone 30N
## First 10 features:
##                          geometry dist_to_nearest
## 1  POLYGON ((-2844767 965638.8...       610530.09
## 2  POLYGON ((1969966 2507177, ...       759565.50
## 3  POLYGON ((-4074219 -551688....       324383.07
## 4  POLYGON ((9454452 -3374190,...            0.00
## 5  POLYGON ((2514100 4188392, ...       445759.94
## 6  POLYGON ((-1934989 6518013,...            0.00
## 7  POLYGON ((1736155 -7916389,...       369668.55
## 8  POLYGON ((4674295 4175327, ...        43115.42
## 9  POLYGON ((4248319 2666428, ...       170073.76
## 10 POLYGON ((-822747.6 -261759...       604320.54

# Plot
plot(x1, reset = FALSE)
plot(x2, border = "red", add = TRUE)
plot(st_connect(x1, x2, result$nn), add = TRUE)

If you can please provide further details on your data (how many features? what is the area of interest? a map?) and/or share a subset of your data, then maybe there will be more specific ideas about why in your case the calculation takes >11 hours.

Answer 2

In R there are warnings, which do execute the code but give you feedback on something you should be considering, and errors, which halt the functions; repeating attributes for all sub-geometries for which they may not be constant is a warning, so you get an output.

Two remarks: 1) If you set the CRS to 4326 st_set_crs(4326) you will get ellipsoidal distances in meters, which may even be better than distances on two UTM zones and 2) Your data, if it is the same of the picture, is not in 32613, but in 32614 and 32615.

library(sf)
library(dplyr)

p = st_sfc(st_point(c(-97.4, 17.4))) %>% st_set_crs(4326)

mat = matrix(c(-96, 17.4, -96,16.5),2,byrow = T)
buff = st_sfc(st_multipoint(mat)) %>%
        st_set_crs(4326)  

a = data.frame(a_value = 1500, geometry = buff) %>% st_as_sf() %>%
    st_transform(32614) %>%   #just to make a buffer polygon
    st_buffer(10000) %>% st_transform(4326)

This is where the warning comes, observe how Units is meters, and crs is 4326:

a %>% st_cast("POLYGON") %>%
      st_distance(p)

Units: [m]
         [,1]
[1,] 169330.8
[2,] 138768.9
Warning message:
In st_cast.sf(., "POLYGON") :
  repeating attributes for all sub-geometries for which they may not be constant

This is considering them a single multigeometry:

a %>% st_distance(p)
Units: [m]
     [,1]
[1,] 138768.9

plot(a, graticule = T, axes = T)
plot(p, add = T)