Calculating distance between spatial object and the nearest border

dist_to_border.md

Calculating distance between spatial object and the nearest border

Quick example showing how to calculate the distance from each census tract to the nearest muncipality border

library(sf)
library(furrr)
library(geobr)
library(mapview)
library(magrittr)
library(data.table)
library(dplyr)
mapview::mapviewOptions(platform = 'leafgl')

# RMSP
rmsp <- read_metro_area(year=2018)
rmsp <- subset(rmsp, name_metro == "RM São Paulo")

# municipalities
mun_all <- read_municipality(code_muni = 'SP', year=2018, simplified = F)
mun_sp <- subset(mun_all, code_muni %in% rmsp$code_muni)

# census tracts
cts <- read_census_tract(code_tract = 'SP', year=2010, simplified = F)
cts <- subset(cts, code_muni %in% rmsp$code_muni)

  # select tracts closer to border
  mun_buff <- st_geometry(mun_sp) %>% st_cast("MULTILINESTRING") %>% st_buffer(dist=.02)
  mun_buff2 <- st_as_sf(mun_buff)
  mun_buff2$id <- '1'
  mun_buff2 = mun_buff2 %>% group_by(id) %>% summarise()
  head(mun_buff2)
  plot(mun_buff2)
  
  cts_candidates <- st_intersection(cts, mun_buff2)
  cts <- cts_candidates
    

mapview(cts_candidates) + mun_buff2 

# function to get distance to nearest muni border
dist_border2 <- function(i){

  message(i)
  temp_ct <- cts[i,]
  # temp_ct <- subset(cts, code_tract=='354680105000074')

  # remove municipality of corresponding census tract
  temp_mun <- subset(mun_all, code_muni != temp_ct$code_muni )

  
  # limit crop extent to improve speed
  buff <- st_buffer(temp_ct, dist=.15)
  crop <- st_crop(temp_mun, buff)
  ## view
  # mapview( cts[1,] ) + buff + temp_mun + crop

  code_tract = temp_ct$code_tract
  crop$distances <- as.numeric(st_distance( temp_ct, crop))
  min_dist <- min( crop$distances )
  closest_muni <- crop$code_muni[which.min(crop$distances)]

  df <- data.frame(rowid = i,
                   code_tract = code_tract,
                   min_dist = min_dist,
                   closest_muni = closest_muni)
  return(df)
}

gc(reset = T, verbose = T, full = T)
# apply function in parallel
future::plan(future::multiprocess)
system.time( output <-  furrr::future_map( .x = 1:nrow(cts), .f = dist_border2, .progress = T) %>% rbindlist() )


# check output
head(output)

Determine which census tracts are in the borders

ind <- sf::st_within(cts$geom, mun$geom)
ind <- sapply(1:length(ind),function(i){length(ind[[i]])})

borders <- ifelse(ind==0,1,0)

cts$borders <- borders