2

I have a grid of type SpatialPolygonsDataFrame (each cell has a set of data attributes). I want to merge each cell which has a surface < X with one of its neighbours. An example is depicted in the picture.

enter image description here

I loop over all cells and see which one has to be merged, then I get its neighbors with gTouches, but I couldn't get it work with functions like gUnaryUnion and aggregate.

Here is my code:

for(ix in(1:length(mygridshp))){
   if(mygridshp@polygons[[ix]]@area < 20000){
      neighbor_cells <- gTouches(mygridspolygons[ix], my_other_cells, byid=TRUE);
      neighbor_cells_ids <- my_other_cells$unit_id[neighbor_cells == TRUE];
      if(length(neighbor_cells_ids) > 0){
          first_neighbor <- mygridshp[mygridshp@data$unit_id==neighbor_cells_ids[1],];
          new_cell <- aggregate(rbind(first_neighbor, mygridshp[ix,]));
          new_cell <- gUnaryUnion(spgeom=union(first_neighbor,mygridshp[ix,]));
      }
   }
}

Any better idea how to do this?

0

I wrote a code that can do this. It's a little bit complicated but it works !

##### function #####
      get_my_first_father <- function(x){
         if(list_of_my_father[x] == x){ return (x); }
         else{ return (get_my_first_father(list_of_my_father[x])); }
      }
      #####
      # for each cell, we look for other cells who are supposed to merge with it (children)
      list_of_my_children <- vector("list", length(mygridspolygons));
      list_of_my_father <- seq.int(length(mygridspolygons));
      for(ix in(1:length(mygridspolygons))){                    
      #for each cell smaller than X, we merge it to its neighbor of the same district
      if(mygridshp@polygons[[ix]]@area < grid_cell_min_area){
         my_district <- mygridshp@data[ix,'dist_code'];
         # all cells of my district
         my_district_cells <- mygridshp[mygridshp@data$dist_code == my_district,];
         # who are my neighbors
         neighbor_cells <- gTouches(mygridspolygons[ix], my_district_cells, byid=TRUE);
         my_district_cells@data$is_neighbor = as.vector(neighbor_cells);
         # my neighbors cells
         my_neighbor_cells <- my_district_cells[my_district_cells@data$is_neighbor == TRUE,];
         # who are not already my children
         my_neighbor_cells <- my_neighbor_cells[!my_neighbor_cells$unit_id %in% 
         unlist(lapply(list_of_my_children[[ix]], function(e) e$unit_id)),];
         # for each cell we calculate the length of the boundaries
         if(length(my_neighbor_cells) > 0){
            line_length <- 0;
            for(lx in(1:length(my_neighbor_cells))){
               my_line <- rgeos::gIntersection(mygridshp[ix,], my_neighbor_cells[lx,], byid = TRUE);
               if(class(my_line) == 'SpatialLines'){
                  new_line_length <- SpatialLinesLengths(my_line);
                  if(new_line_length > line_length){
                     selected_neighbor <- my_neighbor_cells[lx,];
                     line_length <- new_line_length;
                  }
               }
            }
            # getting the big father
            my_father <- get_my_first_father(selected_neighbor$unit_id);
            list_of_my_father[ix] <- my_father;
            list_of_my_children[[my_father]] <- c(list_of_my_children[[my_father]], mygridshp[ix,]);
            # if I have children, they also goes to my father
            if(length(list_of_my_children[[mygridshp[ix,]$unit_id]]) > 0){
               list_of_my_children[[my_father]] <- c(list_of_my_children[[my_father]], list_of_my_children[[mygridshp[ix,]$unit_id]]);
               list_of_my_children[mygridshp[ix,]$unit_id] <- list(NULL);
            }
         }                  
      }
   }


# merge spatial polygons
   my_new_grid <- mygridshp[1,]; # to avoid raster::bind NULL issue
   for(ix in(1:length(list_of_my_children))){
      if(length(list_of_my_children[[ix]]) > 0){
         new_polysp <- as(mygridshp[ix,], 'SpatialPolygons');
         for(lx in(1:length(list_of_my_children[[ix]]))){
            new_polysp <- gUnaryUnion(spgeom=union(new_polysp,list_of_my_children[[ix]][[lx]]));
         }
         new_polyspdf <- SpatialPolygonsDataFrame(new_polysp, mygridshp[ix,]@data, match.ID=FALSE);
         my_new_grid <- raster::bind(my_new_grid, new_polyspdf);    
      }
      else if(mygridshp@polygons[[ix]]@area >= grid_cell_min_area){
         my_new_grid <- raster::bind(my_new_grid, mygridshp[ix,]);
      }
   }

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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