I want to group together polygons from a shapefile a way that all the clusters formed contained approximately the same number of polygons. The problem is that the source polygons have vastly different sizes. Therefore, I cannot simply make a grid and lay it over my polygons.

The final goal of this work is to do some kind of "tiling" of the polygons to call computer heavy analysis on a sub-sample of my data.

The clustering has to be done strictly from a spatial contiguous polygons point of view. I don't care about any attributes in the shape.

There are some similar questions out there (for example here and here), however my question differs as I don't use attribute clustering, have vastly different size polygons and can't use ArcGIS.

I'm open for a solution is R, QGIS, SAGA or GRASS.

The figure below shows an example of what I mean by vastly different size:

enter image description here

with the red circle being:

enter image description here

My data are openly available at : http://www12.statcan.gc.ca/census-recensement/2011/geo/bound-limit/files-fichiers/2016/lda_000b16a_e.zip

I was thinking of scripting the solution so I can change the size of the cluster, but my ideal cluster size would be around 1000 polygons per group.

  • Questions that ask for solutions in more than one system get closed for being too broad. Please focus down on one software system per question. You can ask multiple almost identical questions but replacing R with QGIS etc, just do that with separate posts.
    – Spacedman
    Nov 21, 2018 at 22:03
  • 1
    I have to respectfully disagree. My question isn't broad, it's simply : "How can I cluster polygons spatially a way to have a fix cluster size". I'm not asking for an answer of that question in every language, I'm (not so much) limiting the solution to be doable in any of 4 different software listed (pretty much excluding only ArcGIS, PostGIS and Python, even if the last two would have done...). Maybe I should not have used the qgis, R, saga and grass tag, that can be confusing. Should I remove them?
    – Bastien
    Nov 22, 2018 at 18:43
  • "I'm open for a solution is R, QGIS, SAGA or GRASS.". There are lots of examples here of Qs that ask for solution in one of N systems, and are closed for being too broad, mostly by @PolyGeo - eg: gis.stackexchange.com/questions/303285/… - it doesn't matter how specific the question is, if you ask for solutions in any of two systems, its too broad. See PolyGeo's comments on such questions for more guidance.
    – Spacedman
    Nov 22, 2018 at 23:29

1 Answer 1


I finally found a way, I scripted my own clustering function in R which works relatively great. The idea is to start from the neighbor list and iteratively make groups. The algorithm looks a little bit like that:

  1. Make one group
    1. select first polygon
    2. find it's neighbors
    3. group them
    4. find the neigbors neighbors
    5. group them
    6. repeat until target size reach
  2. Remove all polygons grouped before from the neighbors list
  3. repeat 1 and 2 until no more polygons are availables
  4. Fuse the small leftover groups to majors groups

I coded this algorithm in R in the form of 4 functions. They are available on my github. It may not be perfect but it works well enough for my need and I manage to get it to work on a fairly large shape (~13000 polygons). It's weakness are potentially:

  • Can produce group of size quite different from target size (e.g. aiming at 1000 finishing qith 1400)
  • the shape of the groups are not optimized (like bubble shaped)

Here is the code for the 4 main functions (however, github version is more recent):

## function to make individual groups
make_a_group <- function(neighb_list, gr_size) {

  # make first group
  gr <- list(gr = c(as.numeric(names(neighb_list)[1]), neighb_list[[1]]),
             eval = as.numeric(names(neighb_list)[1]))
  nn <- length(gr$gr)

  # Until you reach target cluster size, do:
  while(nn < gr_size){

    # find neighbours to evaluate
    to_eval <- gr$gr[!gr$gr %in% gr$eval]
    if(length(to_eval)==0) break() # if none, stop
    i = 1
    nb_to_eval <- length(to_eval)

    # and evaluate them sequentially until you did them al or reach target size
    while(nn < gr_size & i<=nb_to_eval){
      gr$gr <- unique(c(gr$gr, neighb_list[[as.character(to_eval[i])]]))
      i = i+1
      nn <- length(gr$gr)

    # update the info on which polygon was evaluated
    gr$eval <- c(gr$eval, to_eval[1:(i-1)])
  lapply(gr, sort)

## function that group the leftover polygons
##  (sometime, a polygon looses all it's neighbors because they were assign other groups)
group_empty_pol <- function(pol_name, neighb_list, the_groups){

  # find the original neighbors of this polygon
  all_ne <-neighb_list[[pol_name]]

  # Find the group which has the most polygons in common with its neighbors
  better_gr <- which.max(sapply(the_groups, function(xx) sum(xx %in% all_ne)))

  # add him to the group
  the_groups[[better_gr]] <- sort(c(as.numeric(pol_name), the_groups[[better_gr]]))

## Function to split the whole shape in groups from neighbours list

cluster_neighbours <- function(neighb_list, gr_size){

  # initialize object
  ll_group <- list()
  neighb_temp <- neighb_list

  # until all are neighbors list is empty, do:

    # Make a group
    res1 <- make_a_group(neighb_list = neighb_temp, gr_size = gr_size)

    # save the group
    ll_group <- append(ll_group, list(res1$gr))

    # Remove from neighbor list all polygons from the just made group
    neighb_temp <- neighb_temp[!names(neighb_temp)%in%as.character(res1$gr)] %>% 
      lapply(function(xx) xx[!xx%in%res1$gr])                                       # remove the used neighbors from neighbor list of left neighbors

    ## managing polygons with no more neighbors
    lost_all_neighb <- neighb_temp[sapply(neighb_temp, length)==0] %>% names

    # if some polygons lost all their neighbors, do:
      for(i in 1:length(lost_all_neighb)){
        #assign the lonely polygon to a group
        ll_group <- group_empty_pol(pol_name = lost_all_neighb[i], neighb_list = neighb_list, the_groups = ll_group)

        # Remove it from the neighbor list
        neighb_temp <- neighb_temp[!names(neighb_temp)%in%lost_all_neighb[i]] %>% 
          lapply(function(xx) xx[!xx%in%as.numeric(lost_all_neighb[i])])       

## final grouping of tiny groups
##  (sometime, some tiny groups are left, we want to fuse them to bigger ones)
fuse_tiny_group <- function(initial_gr, init_sh, hard_min_size){

  # initialize data
  initial_gr <- initial_gr %>% 

  # assign good group to polygon, and dissolve the shape by group
  sh_with_gr <- bind_cols(init_sh,
                          do.call("rbind", mapply(function(gr, index) cbind(index, gr),1:length(initial_gr), initial_gr)) %>% 
                            as.data.frame() %>% 
                          ) %>% 
    select(gr) %>% 
    mutate(gr = as.character(gr)) %>% 
    group_by(gr) %>% 
    summarise(n = n()) %>% 

  # Calculate new neighbors
  adj_gr <- sh_with_gr %>% 
    st_touches() %>% 

  # find the small groups that are not island (which cannot be group to anything)
  small_gr <- sh_with_gr$gr[sh_with_gr$n<hard_min_size]
  not_island <- names(adj_gr)[sapply(adj_gr, length)!=0]
  small_gr <- intersect(small_gr, not_island)

  # for every small group, fuse it with the smaller number of polygons
  for(igr in small_gr) {
    groups_close_by <- adj_gr[[igr]]

    group_selected <- initial_gr %>% 
      subset(.,names(.)%in%groups_close_by) %>% 
      sapply(length) %>% 
      sort() %>%
      names %>% 

    initial_gr[[group_selected]] <- sort(c(initial_gr[[group_selected]], initial_gr[[igr]]))
    initial_gr[[igr]] <- NULL

And here is a example to use it:


DA <- st_read("C:/Users/DXD9163/Desktop/DA_regrouping/DA_CAN_2016_reduced.shp", stringsAsFactors = F) %>% 

DA_qc <- DA[grepl("Quebec",DA$PRNAME),]

DA_qc$area <- st_area(DA_qc)
DA_qc <- DA_qc %>% 
  group_by(DAUID) %>% 

neighb <- st_touches(DA_qc, sparse = T) %>% 

gr1 <- cluster_neighbours(neighb_list = neighb, gr_size = 1000)

sapply(gr1, length)

final_group <- fuse_tiny_group(initial_gr = gr1, init_sh = DA_qc, hard_min_size = 400)

          do.call("rbind", mapply(function(gr, index) cbind(index, gr),1:length(final_group), final_group)) %>% 
            as.data.frame() %>% 
          ) %>% 
  select(gr) %>% 
  mutate(gr = as.character(gr)) %>% 

Which gives:

enter image description here

and zoomed:

enter image description here

  • Given that this algorithm "Can produce group of size quite different from target size" but your Q was "all the clusters formed contained approximately the same number of polygons" I don't see how this is an answer. But then it may be impossible to evenly partition the space subject to the connectivity constraint.
    – Spacedman
    Nov 23, 2018 at 8:49
  • Why not just use the knearneigh and knn2nb functions in the spdep package? This would be quite a bit less convoluted and all you would have to test for would be duplicate neighbors with some sort of inclusion/exclusion criteria. I would add, echoing @Spacedman, that unless you include some sort of optimization algorythm, it will be quite difficult to allocate an even space or uniform cluster size using a connectivity constraint. Nov 26, 2018 at 17:57
  • 2
    @Jeffrey Evans, just tested your proposition. Quickly, it doesn't seems to produce required results. I already got my neighbor list from st_touches. I need groups. Also, it requires points and not polygons. Using centroids in this case with big variance in polygon sizes is not ideal. Anyway, if I didn't understand well your comment, feel free to provide and answer, I'll test it and if it's better than mine I'll accept it and use it.
    – Bastien
    Nov 26, 2018 at 19:33
  • @Bastien this is a nice script. I would like to modify it and use it to select regions for an analysis I am working on, how should I give you credit if we decide to publish the analysis?
    – qdread
    Apr 8, 2019 at 19:10
  • 1
    Hi @qdread, glad it's useful for you. I've put the code on my github under the GPL3 licence: (github.com/BastienFR/clusteRneighbors/blob/master/LICENSE), so I guess you do pretty much as you want. If you find bugs or make improvements, don't hesitate to let me know or create pull request.
    – Bastien
    Apr 9, 2019 at 11:43

Your Answer

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

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