An example of the process required to join objects by most area covered

sf-join-by-most-area-covered.R

options(
  reprex.advertise = FALSE,
  reprex.si = TRUE,
  reprex.style = TRUE,
  reprex.comment = "##",
  reprex.tidyverse_quiet = TRUE
)

reprex({
  # SETUP ----
  library(magrittr)
  library(tidyverse) 
  library(ggplot2)     # devtools::install_github('tidyverse/ggplot2) for `geom_sf()`
  library(sf)
  library(scales)
  library(gridExtra)
  
  # LOAD DATA ----
  
  nc_sf <- 
    st_read(system.file("shape/nc.shp", package="sf")) %>% 
    st_transform(2264) # NC state plane, US foot
  
  # CREATE GRID ----
  
  grid_sf <- 
    tibble(geometry = st_make_grid(nc_sf)) %>% 
    st_sf %>% 
    bind_cols(crossing(X=1:10,Y=1:10)) %>% 
    mutate(CENTROID = map(geometry, st_centroid),
           COORDS = map(CENTROID, st_coordinates),
           COORDS_X = map_dbl(COORDS,1),
           COORDS_Y = map_dbl(COORDS,2),
           COORD_COMB = map2_dbl(COORDS_X,COORDS_Y,~ sum(abs(.x)*1e5,abs(.y)))) %>% 
    arrange(desc(COORD_COMB)) %>% 
    mutate(GRID_ID = purrr::cross2(LETTERS[1:10],str_pad(1:10,2,side = "left",pad = "0")) %>% map_chr(str_c, collapse = "")) %>% 
    as_tibble() %>% 
    st_as_sf()
  
  # SPATIAL JOIN FUNCTION ----
  
  st_drop_geometry <- function (x) 
  {
    if (inherits(x, "sf")) {
      x <- st_set_geometry(x, NULL)
      class(x) <- "data.frame"
      x <- as_tibble(x)
    }
    return(x)
  }
  
  st_intersects_most <- function(x, y){ 
    ints <- st_intersects(x,y)
    
    tibble(
      ROW_ID = imap(ints, ~rep(.y, times = length(.x))) %>% flatten_chr(),
      INTERSECT_DF = flatten_int(ints) %>% map(~y[.x,] %>% st_drop_geometry),
      INTERSECT_AREA = st_intersection(x,y) %>% st_area %>% map_dbl(1)
    ) %>%  
      unnest %>%  
      mutate(geometry = st_geometry(x[ROW_ID,])) %>% 
      st_as_sf
  }
  
  # SPATIAL JOIN ----
  
  nc_grid_sf <- 
    st_intersects_most(nc_sf, grid_sf) %>% 
    group_by(ROW_ID) %>% 
    top_n(1, INTERSECT_AREA) %>% 
    ungroup %>% 
    select(-CENTROID,-COORDS,-COORDS_X,-COORDS_Y,-COORD_COMB) %>%   
    mutate(CENTROID = map(geometry, st_centroid),
           COORDS = map(CENTROID, st_coordinates),
           COORDS_X = map_dbl(COORDS,1),
           COORDS_Y = map_dbl(COORDS,2),
           COORD_COMB = map2_dbl(COORDS_X,COORDS_Y,~ sum(abs(.x)*1e5,abs(.y))))
  
  
  # PLOT RESULTS ----
  
  # Create a discrete color scale for each grid rectangle 
  
  scale_fill_grid <- function(...){
    
    cols <- rep(brewer_pal(type = "qual", palette = 3)(12),times = 9)[1:100]
    
    vals <- grid_sf$GRID_ID
    
    ggplot2::scale_fill_manual(..., 
                               values = setNames(cols,
                                                 vals))
  }
  
  
  # Grid beside counties
  g1 <- ggplot()
  g1 <- g1 + geom_sf(data = grid_sf,  color = "black", fill = "transparent") 
  g1 <- g1 + geom_text(data = grid_sf, mapping = aes(COORDS_X, COORDS_Y, label = GRID_ID))
  g1 <- g1 + coord_sf(crs = st_crs(grid_sf), datum = NA)
  g1 <- g1 + theme_void()
  g1 <- g1 + theme(legend.position = "none")
  
  g2 <- ggplot() 
  g2 <- g2 + geom_sf(data = nc_sf, color = "white")  
  g2 <- g2 + coord_sf(crs = st_crs(nc_sf), datum = NA)
  g2 <- g2 + theme_void()
  g2 <- g2 + theme(legend.position = "none")
  
  # Grid over counties
  
  g3 <- ggplot()
  g3 <- g3 + geom_sf(data = grid_sf, mapping = aes(fill = GRID_ID), alpha = .5, color = "transparent")
  g3 <- g3 + geom_sf(data = nc_sf, color = "black",  fill = "transparent")
  g3 <- g3 + scale_fill_grid()
  g3 <- g3 + geom_text(data = grid_sf, mapping = aes(COORDS_X, COORDS_Y, label = GRID_ID))
  g3 <- g3 + coord_sf(crs = st_crs(grid_sf), datum = NA)
  g3 <- g3 + theme_void()
  g3 <- g3 + theme(legend.position = "none")
  
  g4  <- ggplot()
  g4  <- g4  + geom_sf(data = nc_grid_sf, aes(fill = GRID_ID), alpha = .5, color = "white")
  g4  <- g4  + scale_fill_grid() 
  g4  <- g4  + geom_text(data = nc_grid_sf, mapping = aes(COORDS_X, COORDS_Y, label = GRID_ID), size = 2)
  g4  <- g4  + geom_sf(data = grid_sf, fill = "transparent", color = "black")
  g4  <- g4  + coord_sf(crs = st_crs(grid_sf), datum = NA)
  g4  <- g4  + theme_void()
  g4  <- g4  + theme(legend.position = "none")
  
  
  grid.arrange(g2, g1, ncol=1)
  grid.arrange(g3, g4, ncol=1)
})