calcular "medida de erro"

[dhersz]

No description provided.

[rafapereirabr]

Aqui uma proposta de como calcular a medida. A ideia é que na etapa de agregação do cnefe, a gente calculte área do poligono concavo.

# geracao dos dados do reprex 

devtools::load_all(".")
library(data.table)
library(dplyr)
library(sf)
library(sfheaders)
library(arrow)
library(mapview)
mapview::mapviewOptions(platform = 'leafgl')



tabela1 <- geocodebr::listar_dados_cache()[7]

dt <- arrow::open_dataset( tabela1 ) |>
  dplyr::filter(estado == 'RJ') |>
  dplyr::filter(municipio == "RIO DE JANEIRO") |>
  dplyr::filter(cep %in% c("22620-110", "20521-470")) |>
  dplyr::collect()


# funcao para calcular area de poligono concavo
get_concave_area <- function(lat_vec, lon_vec){
  
  # lat_vec <- dt$lat
  # lon_vec <- dt$lon
  temp_matrix <- matrix( c(lon_vec, lat_vec), ncol = 2, byrow = FALSE)
  temp_sf <- sfheaders::sf_point(temp_matrix, keep = T)
  sf::st_crs(temp_sf) <- 4674
  
  poly <- concaveman::concaveman(points = temp_sf)
  # poly <- sf::st_convex_hull(x  = st_union(temp_sf))
  # poly <- sf::st_concave_hull(x = st_union(temp_sf), ratio = 0.5)
  
  sf::st_crs(poly) <- 4674
  
  area_m2 <- sf::st_area(poly)
  area_m2 <-  as.numeric(area_m2) |> round()
  
  return(area_m2)
}


# agrega e calcula area
dt[, .(lat = mean(lat),
       lon = mean(lon),
       area_m2 = get_concave_area(lat, lon)
       ), 
   by = cep]

#>          cep      lat       lon area_m2
#>       <char>    <num>     <num>   <num>
#> 1: 22620-110 -23.0111 -43.30139    8763
#> 2: 20521-470 -22.9296 -43.22587    1025

Para uma inspeção visual desses resultados, é só usar o codigo abaixo:

# funcao para calcular poligono concavo
get_concave_poly <- function(lat_vec, lon_vec){
  
  # lat_vec <- dt$lat
  # lon_vec <- dt$lon
  temp_matrix <- matrix( c(lon_vec, lat_vec), ncol = 2, byrow = FALSE)
  temp_sf <- sfheaders::sf_point(temp_matrix, keep = T)
  sf::st_crs(temp_sf) <- 4674
  
  poly <- concaveman::concaveman(points = temp_sf)
  # poly <- sf::st_convex_hull(x  = st_union(temp_sf))
  # poly <- sf::st_concave_hull(x = st_union(temp_sf), ratio = 0.5)
  
  sf::st_crs(poly) <- 4674

  return( sf::st_geometry(poly) )
}

poly <- dt[, .(lat = mean(lat),
               lon = mean(lon),
               geometry = get_concave_poly(lat, lon)
               ), 
           by = cep]

poly <- sf::st_as_sf(poly)

dt_sf <- sfheaders::sf_point(
  obj = dt,
  x = 'lon',
  y = 'lat',
  keep = TRUE
)

sf::st_crs(dt_sf) <- 4674

mapview::mapview(dt_sf, zcol='cep') + poly