Inspirado num antigo post de BlakeRMills, criador do pacote {MetBrewer}, criei um mapa com a altitude das ruas em Brasília.
Em breve farei um post com tutorial detalhado e também pretendo replicar este tipo de mapa para outras cidades interessantes. O código para replicar o gráfico está abaixo.
Code
library(dplyr)library(sf)library(ggplot2)library(osmdata)library(purrr)sf::sf_use_s2(FALSE)sysfonts::font_add_google("Roboto Condensed", "Roboto Condensed")showtext::showtext_auto()url ="https://pt.wikipedia.org/wiki/Lista_de_municípios_do_Brasil_por_população"tab = xml2::read_html(url) |> rvest::html_table() |> purrr::pluck(2)as_numeric_char =Vectorize(function(x) { ls = stringr::str_extract_all(x, "[[:digit:]]") y =paste(ls[[1]], collapse ="")as.numeric(y)})clean_tab = tab |> janitor::clean_names() |>rename(code_muni = codigo_ibge,name_muni = municipio,rank = posicao,name_state = unidade_federativa,pop = populacao ) |>filter(name_muni !="Brasil") |>mutate(code_muni =as.numeric(code_muni),pop =as_numeric_char(pop),rank =rank(-pop) )top20 =slice_max(clean_tab, pop, n =20)get_border =function(city) {#> Encontra o código do município code_muni = top20 |>filter(name_muni == city) |>pull(code_muni) |>unique()stopifnot(length(code_muni) ==1)#> Baixa o shapefile do município border = geobr::read_municipality(code_muni, showProgress =FALSE)return(border)}get_streets =function(city, border) {#> Encontra o nome da Unidade Federativa nome_uf = top20 |>filter(name_muni == city) |>pull(name_state)#> Monta o nome do local name_place = stringr::str_glue("{city}, {nome_uf}, Brazil")#> Monta a query place =opq(bbox =getbb(name_place))#> Importa todas as principais vias da cidade# streets = add_osm_feature(# place,# key = "highway",# value = c(# "motorway", "primary", "motorway_link", "primary_link",# "secondary", "tertiary", "secondary_link", "tertiary_link",# "residential"# )# ) streets =add_osm_feature(place, key ="highway")#> Converte o dado streets = streets %>%osmdata_sf() %>% .$osm_lines %>%select(osm_id, name) %>%st_transform(crs =4674)#> Enconrtra a intersecção entre as estradas e o limites do município streets_border =st_intersection(streets, border)# out = list(streets = streets, streets_border = streets_border)return(streets_border)}get_elevation =function(border, z =8) { altitude = elevatr::get_elev_raster(border, z = z, clip ="bbox") altitude = raster::rasterToPolygons(altitude) altitude =st_as_sf(altitude)names(altitude)[1] ="elevation" altitude =st_transform(altitude, crs =4674) altitude =suppressWarnings(st_intersection(altitude, border)) altitude =filter(altitude, st_is_valid(altitude))return(altitude)}add_jenks_breaks =function(shp, k =7, round =TRUE, r =0) {#> Classifica os dados de altitude em k grupos segundo o algo. de Jenks jbreaks = BAMMtools::getJenksBreaks(shp$elevation, k = k)#> Arredonda os números para chegar numa legenda menos quebradaif (round) { jbreaks =round(jbreaks, r) }#> Cria a coluna 'jenks_group' que classifica cada valor num grupo shp =mutate(shp, jenks_group =cut(elevation, jbreaks))#> Verifica se todas as observações tem um grupo check =any(is.na(shp$jenks_group))if (check) {warning("Some observations have failed to be grouped") }#> Transforma os groups em legendas labels =get_jenks_labels(jbreaks)#> Retorna o output numa lista out =list(shp = shp, labels = labels)return(out)}get_jenks_labels =function(x) { labels =paste(x, x[-1], sep ="–") labels[1] =paste("<", x[2]) labels[length(labels)] =paste(">", max(x))return(labels)}get_streets_altitude =function(altitude, streets) {stopifnot(any(colnames(altitude) %in%"jenks_group"))#> Get all groups groups =levels(altitude$jenks_group)#> For each group get the full polygon and join with streets join_streets =function(group) { poly = altitude %>%filter(jenks_group == group) %>%st_union(.) %>%st_as_sf() %>%st_make_valid() joined =suppressWarnings(st_intersection(streets, poly))return(joined) }#> Apply the function to all groups street_levels = purrr::map(groups, join_streets)#> Bind all results together out =bind_rows(street_levels, .id ="level")return(out)}map_plot =function(shp, labels, title, showtext =TRUE) { colors = viridis::plasma(n =length(labels) +1) colors = colors[-length(colors)] font =ifelse(showtext ==TRUE, "Roboto Condensed", "sans") plot =ggplot(data = shp) +geom_sf(aes(color = level, fill = level), linewidth =0.2) +scale_color_manual(name ="Altitude",labels = labels,values = colors) +scale_fill_manual(name ="Altitude",labels = labels,values = colors) +guides(fill =guide_legend(nrow =1), color =guide_legend(nrow =1)) +ggtitle(title) + ggthemes::theme_map() +coord_sf() +theme(plot.title =element_text(size =30,hjust =0.5,family = font ),legend.title =element_text(size =20,family = font,color ="gray10" ),legend.text =element_text(size =14,family = font,color ="gray10" ),legend.position ="top",legend.direction ="horizontal",plot.background =element_rect(color =NA, fill ="#f6eee3"),panel.background =element_rect(color =NA, fill ="#f6eee3"),legend.background =element_rect(color =NA, fill ="#f6eee3") )return(plot)}map_altitude =function(city, k, z) {#> Importa o shape do limite do municípiomessage("Importando os limites do município: ", city) city_border =get_border(city)#> Importa as principais vias da cidade e junta com o limite do munimessage("Importando as vias.") city_street =get_streets(city, city_border)#> Importa a altitude da cidademessage("Importando a altitude.") city_elevation =suppressMessages(get_elevation(city_border, z = z))#> Classifica a altitude em gruposmessage("Classificando e juntando os shapefiles.") jenks =add_jenks_breaks(city_elevation, k = k) city_elevation = jenks[["shp"]] labels = jenks[["labels"]]#> Junta a altitude (agrupada) com as vias city_street_elevation =get_streets_altitude(city_elevation, city_street)#> Monta o mapa finalmessage("Gerando o mapa final.") plot =map_plot(city_street_elevation, labels = labels, title = city)message("Feito.")#> Retorna o output numa lista out =list(shp = city_street_elevation,streets = city_street,elevation = city_elevation,plot = plot )return(out)}export_citymap =function(city, w =14, h =16, ...) { plot =map_altitude(city, ...)$plotif (is.numeric(city)) { name_city = cities_brasil |>filter(code_muni == city) |>pull(name_muni) } elseif (is.character(city)) { name_city = city } name_file = glue::glue("elevation_{janitor::make_clean_names(name_city)}.svg" )ggsave( here::here("graphics/altitude", name_file), plot,width = w,height = h ) name_file = glue::glue("elevation_{janitor::make_clean_names(name_city)}.png" )ggsave( here::here("graphics/altitude", name_file), plot,width = w,height = h )}safe_export = purrr::safely(export_citymap)params =tribble(~city, ~z, ~k,#----------------#---#---#"São Paulo", 8, 7,"Rio de Janeiro", 9, 7,"Brasília", 8, 7,"Fortaleza", 11, 7,"Salvador", 10, 7,"Belo Horizonte", 10, 7,"Manaus", 6, 7,"Curitiba", 10, 7,"Recife", 10, 7,"Goiânia", 8, 7,"Porto Alegre", 8, 7,"Belém", 9, 7,"Guarulhos", 9, 7,"Campinas", 9, 7,"São Luís", 9, 7,"Maceió", 11, 8,"Campo Grande", 7, 7,"São Gonçalo", 8, 7,"Teresina", 10, 7,"João Pessoa", 9, 7,"Joinville", 9, 7 )#> Exportar todas as cidades listadas acima#pmap(params, safe_export)#> Exportar uma cidade em particular#safe_export("Brasília", z = 8, k = 7)
