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. Visualization

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.

 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.

  • 1
    You don't have to cast from MULTIPOLYGON to POLYGON and also don't have to calculate the centroids. Function st_distance gives the shortest distance to the nearest boundary point between (MULTI)POLYGONs. So st_distance(ANPs, DISTAgs) should do, returning a distance matrix of shortest distances between each pair of features in the two layers. For more specific example, please share a sample of the data you are working with. Nov 8, 2020 at 6:54
  • I've re coded everything in order to generate GEODIST a 300 row data frame that encapsulates the geometries to analyze. I've also have created an index to order the ANPs geometries to match the ANPs nearest observation for each GEODIST observation. With that in mind, I generated st_distance(GEODIST, ANPs, by_element=TRUE) but it has taken around 11 hours without giving a result? I based this operation on this question: stackoverflow.com/questions/53854803/….
    – socapal
    Nov 12, 2020 at 3:10
  • 1
    Thanks! I've posted an answer, please see below. Nov 12, 2020 at 9:23

2 Answers 2


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.


# Sample data
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)
## 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.

  • I created a Gitub repository in order to share some of the files with you. I tried running st_nn but it doesn't seem to work correctly. I'm very grateful that you have helped me with all of this!
    – socapal
    Nov 13, 2020 at 8:48
  • Thanks, I briefly went over the code but don't really understand what are you trying to do, and what the problem is. What do you mean by "doesn't seem to work correctly"? Can you please elaborate, perhaps in a new question? Nov 13, 2020 at 17:35
  • 1
    Hello! A bit dumb on my part, but I didn't even think about generating centroids in order to run them using st_nn(). I now have, it worked magnifically! Only one minute for all the computations which had taken 11 hours before. Great package :)
    – socapal
    Nov 20, 2020 at 6:26
  • Thanks! Keep in mind that "nearest centroid" does not necessarily imply "nearest feature" (in terms of minimal distance between polygon boundaries). The amount of bias when dealing with centroids depends on how large the distances are, how much are the polygons circular rather than complex shape, etc. But it surely is much faster :-) Nov 20, 2020 at 17:39

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.


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)) %>%

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") %>%

Units: [m]
[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,] 138768.9

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

enter image description here

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

Not the answer you're looking for? Browse other questions tagged or ask your own question.