library(data.table)
library(tidygeocoder)
library(sf)
library(stringr)
library(osrm)
library(mapview)
library(ggplot2)
library(aopdata)
Acessibilidade à Saúde
Neste post, vou mapear a acessibilidade a hospitais e leitos em São Paulo. Para avaliar quantitativamente o nível de acessibilidade vou montar uma métrica bastante simples: o tempo mínimo necessário que se leva para chegar no hospital/leito mais próximo, considerando um deslocamento de bicicleta1.
Para tornar o problema tratável, divido a cidade em hexágonos, no padrão H3, em resolução 9. Esta resolução tem um tamanho aproximado de 1km2 e estratifica a cidade em cerca de 15 mil subáreas. Seria possível reduzir o número de hexágonos cruzando este grid com dados de população do Censo. Contudo, como os dados do Censo a nível de setor censitário estão consideravelmente defasados, utilizo o grid inteiro.
Para cada hospital/leito, contruo isócronas de 10, 15, 20, 25 e 30 minutos, considerando um deslocamento de bicicleta. A partir desta métrica simples, pode-se construir refinamentos como ponderar o número de leitos acessíveis em relação a população que reside em cada região. Evidentemente, pode-se considerar também outros modos de transporte e intervalos de tempo.
Setup
Como os dados são relativamente grandes vou usar o data.table
para a manipulação dos dados. O pacote tidygeocoder
é utilizado para georeferenciar alguns dos endereços que estão com lat/long ausente.
Dados
As bases de dados são de livre acesso no site dados.gov.br. Especificamente, vamos usar a base de Cadastros de CNES e de Hospitais e Leitos. Normalmente, eu importaria estes dados usando webscrapping, mas neste caso acho mais simples baixar os arquivos manualmente. Como são apenas duas tabelas que precisam ser baixadas, acho que seria um exagero fazer um rotina para importar os dados via scraping. Além disso, olhando o código-fonte da página, o download é feito via uma URL dinâmica, o que exigiria o uso do RSublime
, ou, alternativamente, o uso do BeautifulSoup
em Python2.
= fread("cnes_estabelecimentos.csv")
cnes = fread("Leitos_2024.csv", encoding = "Latin-1") leitos
Limpeza
Os códigos abaixo mostram a rotina de limpeza necessária. Como os dados brutos são relativamente bem estruturados não é necessário muito esforço neste passo. O primeiro código abaixo seleciona apenas as colunas necessárias da base de Cadastro de CNES para São Paulo.
# Simplifica e limpa o nome das colunas
= janitor::clean_names(cnes)
cnes # Vetor para renomear e selecionar as colunas
<- c("code_cnes", "lng", "lat")
sel_cols # Altera o nome das colunas
setnames(cnes, c("co_cnes", "nu_longitude", "nu_latitude"), sel_cols)
# Seleciona as colunas apenas para Sao Paulo
<- cnes[co_ibge == 355030, ..sel_cols] cnes
O segundo código abaixo mostra os passos para limpar a base de hospitais e leitos. Por fim, as duas bases são combinadas via o identificador code_cnes
.
# Vetor para renomear as colunas
# Remove o prefixo do nome das colunas
<- str_remove_all(names(leitos), "^[A-Z]{2}_")
new_names # Cria novos nomes paras as colunas
<- janitor::make_clean_names(new_names)
new_names setnames(leitos, names(leitos), new_names)
setnames(leitos, c("cnes", "endereco"), c("code_cnes", "numero"))
# Vetor para selecionar colunas
<- c(
sel_cols "code_cnes", "nome_estabelecimento", "logradouro", "numero", "bairro", "cep",
"leitos_existentes"
)
<- leitos[uf == "SP" & municipio == "SAO PAULO", ..sel_cols]
leitos
<- merge(leitos, cnes, by = "code_cnes") dat
Geocoding
Há poucos casos que não foram identificados no cruzamento dos dados. Especificamente, apenas três hospitais não foram identificados. Note que a função unique
remove as linhas duplicadas pois a base de leitos não unicamente identificada por estabelecimento.
<- dat[is.na(lng) | is.na(lat)]
missing_leitos <- unique(missing_leitos[, .(code_cnes, logradouro, numero, bairro)])
missing_leitos
print(missing_leitos)
Key: <code_cnes>
code_cnes logradouro numero bairro
<int> <char> <char> <char>
1: 955736 RUA GENERAL CHAGAS SANTOS 314 VILA DA SAUDE
2: 2688603 PAULISTA 200 BELA VISTA
3: 7252455 VERGUEIRO 235 LIBERDADE
Para georeferenciar estas observações uso o tidygeocoder
.
missing_leitos[,:= glue::glue_data(.SD, "{numero}, {logradouro}, {bairro}, São Paulo, SP")
endereco
]
<- geocode(missing_leitos, address = endereco) geo_leitos
Pelo mapa abaixo, fica claro que o processo teve sucesso.
= st_as_sf(geo_leitos, coords = c("long", "lat"), crs = 4326)
test mapview(test)
Por fim, o código abaixo remove as linhas da base original, que tinham lat/long ausente, e substitui elas pelas novas linhas encontradas acima. Algum cuidado é necessário já que a base original não é identificada por estabelecimento. Por fim, faço um teste simples para verificar que todas as observações tem lat/lng.
setnames(geo_leitos, "long", "lng")
setDT(geo_leitos)
<- geo_leitos[, .(code_cnes, lng, lat)]
geo_leitos
<- dat[is.na(lng) | is.na(lat)]
fix_dat set(fix_dat, j = c("lng", "lat"), value = NULL)
<- merge(fix_dat, geo_leitos, all.x = TRUE, by = "code_cnes")
fix_dat
<- dat[!(is.na(lng) | is.na(lat))]
dat <- rbind(dat, fix_dat)
dat
# Verifica se há lat/lngs faltando
nrow(dat[is.na(lng) | is.na(lat)])
[1] 0
Isócronas
Isócronas são polígonos que representam áreas de alcance dentro de um período de tempo pré-determinado. Isto é, para um determinado ponto no espaço \(s_{i}\) a isócrona \(I(s_{i}, m, t)\) representa todos os pontos que se pode chegar, partindo de \(s_{i}\), usando o modo de transporte \(m\) em \(t\) minutos. Para este exemplo considero apenas deslocamentos de bicicleta em intervalos de 10, 15, 20, 25, e 30 minutos.
O código abaixo exemplifica a construção de uma isócrona em torno de um hospital.
= osrmIsochrone(
iso_test 1, .(lng, lat)],
dat[breaks = seq(from = 10, to = 30, by = 5),
osrm.profile = "bike"
)
mapview(iso_test, zcol = "isomax")
Ao todo, temos 224 pontos que representam um estabelecimento único. Ao todo, calcula-se 1140 isócronas.
<- unique(dat[, .(code_cnes, lng, lat)])
locations nrow(locations)
[1] 224
<- st_as_sf(
geo_locations
locations[, .(lng, lat)],coords = c("lng", "lat"),
crs = 4326
)
mapview(geo_locations)
O código abaixo importa as isócronas para todos os pontos. Note que este código leva um tempo considerável para rodar.
= function(dat) {
get_isochrone
= osrm::osrmIsochrone(
iso
dat[, .(lng, lat)],breaks = seq(from = 10, to = 30, by = 5),
osrm.profile = "bike"
)
return(iso)
}
<- split(locations, by = "code_cnes")
locations
<- lapply(locations, get_isochrone) geo_locations
Grid
Para montar o grid, uso o pacote h3jsr
que implementa as funções da biblioteca H3 da Uber dentro do R (via JavaScript). Importo o shapefile da cidade de São Paulo do IBGE via pacote geobr
.
Vale notar que é possível conseguir um grid H3 em resolução 9 diretamente do pacote aopdata::read_grid()
, que oferece os dados do projeto de Acesso a Oportunidades do IPEA. Este atalho, contudo, está disponível apenas para as cidades que entraram no estudo. Assim, o código abaixo é mais geral, pois funciona para qualquer cidade do Brasil.
library(h3jsr)
library(geobr)
= read_municipality(3550308, simplified = FALSE, showProgress = FALSE)
spo = polygon_to_cells(spo, res = 9, simple = FALSE)
grid
= polygon_to_cells(spo, res = 9, simple = FALSE)
polyfill = unlist(str_split(unlist(polyfill$h3_addresses), ", "))
index_h3
= data.frame(
grid id_hex = index_h3
)
= cell_to_polygon(grid)
h3grid = st_as_sf(h3grid, crs = 4326)
h3grid $id_hex = index_h3 h3grid
Fazendo a interseção espacial entre o grid H3 e as isócronas, calcula-se para cada hexágono o tempo mínimo necessário para chegar no hospital/leito mais próximo.
<- dplyr::bind_rows(geo_locations, .id = "code_cnes") isocronas
= function(dest) {
match_h3_destination
= isocronas |>
iso ::filter(code_cnes == dest)
dplyr
= lapply(1:5, \(i) {
hex_codes
= dplyr::filter(iso, id == i) |>
idhex ::polygon_to_cells(res = 9) |>
h3jsrunlist()
return(data.frame(id_hex = idhex))
})
names(hex_codes) = as.character(c(10, 15, 20, 25, 30))
= rbindlist(hex_codes, idcol = "isomax")
h3iso
}
<- unique(isocronas$code_cnes)
cnes_codes
<- parallel::mclapply(cnes_codes, match_h3_destination) od_table
Agora, para cada hexágono, encontra-se o leito/hospital mais próximo
# Consolida tabela que faz o match das origens e destinos (hospitais)
names(od_table) <- cnes_codes
<- rbindlist(od_table, idcol = "code_cnes")
od # Encontra o hospital mais próximo para cada ponto
:= as.numeric(isomax)]
od[, isomax <- od[, .SD[which.min(isomax)], by = "id_hex"]
min_table
# Grid H3
<- setDT(st_drop_geometry(h3grid))
dtgrid <- merge(dtgrid, min_table, by = "id_hex", all.x = TRUE)
timetable
# Junta os dados com o grid
<- dplyr::left_join(h3grid, timetable, by = "id_hex")
acesso_saude # Troca NAs por > 30
<- acesso_saude |>
acesso_saude ::mutate(
dplyrisomax = ifelse(is.na(isomax), 30, isomax)
)
Resultado
O mapa abaixo mostra o resultado final. De maneira geral, a área central da cidade, o Centro Expandido, está relativamente bem atendido de hospitais numa distância de 10 a 15 minutos até o hospital mais próximo. As áreas periféricas da cidade tem indicadores bastante piores, contudo, vale notar que estas áreas potencialmente podem ser atendidas por hospitais em municípios vizinhos, que fazem conurbação com São Paulo.
Como comentado inicialmente, esta métrica é bastante rudimentar. Uma melhoria interessante seria ponderar o número total de leitos disponíveis em relação a população atendida em cada região.
Code
ggplot(acesso_saude) +
geom_sf(aes(fill = isomax, color = isomax)) +
scale_fill_distiller(
name = "",
palette = "RdBu",
labels = c("Até 10 min.", "15 min.", "20 min", "25 min", "30 min. ou mais")) +
scale_color_distiller(
name = "",
palette = "RdBu",
labels = c("Até 10 min.", "15 min.", "20 min", "25 min", "30 min. ou mais")) +
labs(
title = "Acesso a Hospitais",
subtitle = "Tempo mínimo necessário para chegar ao hospital/leito mais próximo.",
caption = "Fonte: Open Source Routing Machine (osrm), Dados Abertos (Ministério da Saúde)."
+
) ::theme_map() +
ggthemestheme(
legend.position = "top",
legend.justification = 0.5,
legend.key.size = unit(1, "cm"),
legend.key.width = unit(1.75, "cm"),
legend.text = element_text(size = 12),
legend.margin = margin(),
plot.margin = margin(10, 5, 5, 10),
plot.title = element_text(size = 18, hjust = 0.5),
plot.subtitle = element_text(size = 10, hjust = 0.5)
)
Footnotes
Seria possível também fazer o deslocamento de carro, mas, na minha experiência, o deslocamento de bicicleta aproxima bem o deslocamento de carro com trânsito. As isócronas de deslocamento de carro em São Paulo - sem ajustes - costumam ser bastante “otimistas” com relação às possibilidades de deslocamento.↩︎
Caso fosse necessário seria possível combinar a rotina de importação em
Python
com a rotina de análise emR
usandoreticulate
.↩︎