Hasta acá trabajamos con datos descargados de portales open data oficiales de diferentes Ciudades del mundo, pero cabe destacar que esa no es la única fuente de información que existe para analizar Ciudades, ya que por suerte hay varios repositorio online de donde podemos descargarnos información georreferenciada, cómo por ejemplo: OpenStreetMap (OSM, https://www.openstreetmap.org/)
Pero se estarán preguntando ¿Qué es y cómo funciona OSM? La respuesta: Es una plataforma creada por una gran comunidad de colaboradores que aportan datos geoespaciales de muchas Ciudades del mundo (calles, rutas, restaurants, museos, hospitales, escuelas, etc) y los mantienen actualizados. Esta información es consumida por miles de sitios web y aplicaciones ya que puede ser utilizada libremente para cualquier propósito siempre y cuando se respeten los derechos de autor y licencia.
Por suerte, en R existe un paquete llamado osmdata
que nos permitirá conectarnos directamente a OSM y descargarnos toda la información que necesitemos. Así que, comencemos activando las ya conocidas tidyverse
, sf
y ggmap
e instalemos y activemos osmdata
y leaflet
que nos permitirá realizar mapas interactivos:
#install.packages("tidyverse")
library(tidyverse)
#install.packages("sf")
library(sf)
#install.packages("ggmap")
library(ggmap)
#install.packages("osmdata")
library(osmdata)
#install.packages("leaflet")
library(leaflet)
Antes de comenzar a descargar información, tomemos unos minutos para analizar y conocer que información tenemos disponible. En este este link podrán encontrar una lista detallada de las categorías y tipo de elementos que existen en OSM, pero si quisiésemos también podríamos hacer algunas consultas desde acá con available_features()
y available_tags()
. La primer función nos permitirá obtener el encabezado de la información que contiene la base de datos que nos podemos descargar:
available_features()
## [1] "4wd_only" "abandoned"
## [3] "abutters" "access"
## [5] "addr" "addr:city"
## [7] "addr:conscriptionnumber" "addr:country"
## [9] "addr:district" "addr:flats"
## [11] "addr:full" "addr:hamlet"
## [13] "addr:housename" "addr:housenumber"
## [15] "addr:inclusion" "addr:interpolation"
## [17] "addr:place" "addr:postcode"
## [19] "addr:province" "addr:state"
## [21] "addr:street" "addr:subdistrict"
## [23] "addr:suburb" "admin_level"
## [25] "aeroway" "agricultural"
## [27] "alt_name" "amenity"
## [29] "area" "atv"
## [31] "backward" "barrier"
## [33] "basin" "bdouble"
## [35] "bicycle" "bicycle_road"
## [37] "biergarten" "boat"
## [39] "border_type" "boundary"
## [41] "bridge" "building"
## [43] "building:fireproof" "building:flats"
## [45] "building:levels" "building:min_level"
## [47] "building:soft_storey" "bus_bay"
## [49] "busway" "charge"
## [51] "construction" "covered"
## [53] "craft" "crossing"
## [55] "crossing:island" "cuisine"
## [57] "cutting" "cycleway"
## [59] "denomination" "destination"
## [61] "diet" "direction"
## [63] "dispensing" "disused"
## [65] "disused:shop" "drink"
## [67] "drive_in" "drive_through"
## [69] "ele" "electric_bicycle"
## [71] "electrified" "embankment"
## [73] "embedded_rails" "emergency"
## [75] "end_date" "entrance"
## [77] "est_width" "fee"
## [79] "fire_object:type" "fire_operator"
## [81] "fire_rank" "foot"
## [83] "footway" "ford"
## [85] "forestry" "forward"
## [87] "frequency" "fuel"
## [89] "gauge" "golf_cart"
## [91] "goods" "hazmat"
## [93] "healthcare" "healthcare:counselling"
## [95] "healthcare:speciality" "height"
## [97] "hgv" "highway"
## [99] "historic" "horse"
## [101] "ice_road" "incline"
## [103] "industrial" "inline_skates"
## [105] "inscription" "internet_access"
## [107] "junction" "kerb"
## [109] "landuse" "lanes"
## [111] "lanes:bus" "lanes:psv"
## [113] "layer" "leaf_cycle"
## [115] "leaf_type" "leisure"
## [117] "lhv" "lit"
## [119] "location" "man_made"
## [121] "maxaxleload" "maxheight"
## [123] "maxlength" "maxspeed"
## [125] "maxstay" "maxweight"
## [127] "maxwidth" "military"
## [129] "minspeed" "mofa"
## [131] "moped" "motor_vehicle"
## [133] "motorboat" "motorcar"
## [135] "motorcycle" "motorroad"
## [137] "mountain_pass" "mtb:description"
## [139] "mtb:scale:imba" "mtb_scale"
## [141] "name" "name:left"
## [143] "name:right" "narrow"
## [145] "natural" "noexit"
## [147] "non_existent_levels" "note"
## [149] "nudism" "office"
## [151] "official_name" "old_name"
## [153] "oneway" "opening_hours"
## [155] "operator" "organic"
## [157] "oven" "overtaking"
## [159] "parking:condition" "parking:lane"
## [161] "passing_places" "place"
## [163] "power" "priority_road"
## [165] "produce" "proposed"
## [167] "protected_area" "psv"
## [169] "public_transport" "railway"
## [171] "railway:preserved" "railway:track_ref"
## [173] "recycling_type" "ref"
## [175] "religion" "residential"
## [177] "roadtrain" "route"
## [179] "sac_scale" "service"
## [181] "service_times" "shelter_type"
## [183] "shop" "sidewalk"
## [185] "site" "ski"
## [187] "smoothness" "social_facility"
## [189] "sorting_name" "speed_pedelec"
## [191] "start_date" "step_count"
## [193] "substation" "surface"
## [195] "tactile_paving" "tank"
## [197] "tidal" "toilets:wheelchair"
## [199] "toll" "tourism"
## [201] "tracks" "tracktype"
## [203] "traffic_calming" "traffic_sign"
## [205] "trail_visibility" "trailblazed"
## [207] "trailblazed:visibility" "tunnel"
## [209] "turn" "type"
## [211] "usage" "vehicle"
## [213] "vending" "voltage"
## [215] "water" "wheelchair"
## [217] "wholesale" "width"
## [219] "winter_road" "wood"
Dentro de las 200 variables que tiene OSM incluye tanto las características primarias (primary features), como las adicionales (aditional properties). En cada una podremos encontrar keys (por ejemplo, amenity, shop, railway, etc.) y tags que podemos consultarlos de la siguiente forma:
available_tags("amenity")
## [1] "animal_boarding" "animal_breeding" "animal_shelter"
## [4] "arts_centre" "atm" "baby_hatch"
## [7] "baking_oven" "bank" "bar"
## [10] "bbq" "bench" "bicycle_parking"
## [13] "bicycle_rental" "bicycle_repair_station" "biergarten"
## [16] "boat_rental" "boat_sharing" "brothel"
## [19] "bureau_de_change" "bus_station" "cafe"
## [22] "car_rental" "car_sharing" "car_wash"
## [25] "casino" "charging_station" "childcare"
## [28] "cinema" "clinic" "clock"
## [31] "college" "community_centre" "conference_centre"
## [34] "courthouse" "crematorium" "dentist"
## [37] "dive_centre" "doctors" "dog_toilet"
## [40] "drinking_water" "driving_school" "embassy"
## [43] "events_venue" "fast_food" "ferry_terminal"
## [46] "fire_station" "food_court" "fountain"
## [49] "fuel" "funeral_hall" "gambling"
## [52] "give_box" "grave_yard" "grit_bin"
## [55] "gym" "hospital" "hunting_stand"
## [58] "ice_cream" "internet_cafe" "kindergarten"
## [61] "kitchen" "kneipp_water_cure" "language_school"
## [64] "library" "lounger" "love_hotel"
## [67] "marketplace" "monastery" "motorcycle_parking"
## [70] "music_school" "nightclub" "nursing_home"
## [73] "parking" "parking_entrance" "parking_space"
## [76] "pharmacy" "photo_booth" "place_of_mourning"
## [79] "place_of_worship" "planetarium" "police"
## [82] "post_box" "post_depot" "post_office"
## [85] "prison" "pub" "public_bath"
## [88] "public_bookcase" "public_building" "ranger_station"
## [91] "recycling" "refugee_site" "restaurant"
## [94] "sanitary_dump_station" "school" "shelter"
## [97] "shower" "social_centre" "social_facility"
## [100] "stripclub" "studio" "swingerclub"
## [103] "taxi" "telephone" "theatre"
## [106] "toilets" "townhall" "toy_library"
## [109] "university" "vehicle_inspection" "vending_machine"
## [112] "veterinary" "waste_basket" "waste_disposal"
## [115] "waste_transfer_station" "water_point" "watering_place"
Podemos ver que, por ejemplo dentro de amenity se encuentran bares, restaurantes, heladerías, hospitales, bibliotecas, etc.
Listo, hasta acá ya tenemos un poco más de información sobre el contenido posible de descargar. Ahora elijamos de que lugar del mundo queremos descargar la información.
Llegó el momento de decidir el sitio/lugar de donde queremos traernos información ya que no podremos descargar todos los datos que tiene OSM porque son millones! Para esto debemos crear un cuadro delimitador o bounding box con la función getbb()
:
bbox_mdp <- getbb("Mar del Plata, Buenos Aires, Argentina")
La función anterior nos devolverá 4 coordenadas (long max, long min, lat max, lat min) que son los límites de la información a descargar. Revisemos el resultado:
bbox_mdp
## min max
## x -57.65911 -57.51784
## y -38.12733 -37.91325
Aprovechemos lo que aprendimos en la clase anterior y descarguemos un mapa de fondo a partir de get_stamenmap()
:
mapa_mdp <- get_stamenmap(bbox = bbox_mdp,
zoom=13)
ggmap(mapa_mdp)
Cabe destacar que, si ajustamos el parámetro format_out="sf_polygon en getbb()
podremos obtener el polígono del límite territorial de la Ciudad de la que vamos a descargar datos:
mdp_polygon <- getbb("Mar del Plata, Buenos Aires, Argentina", format_out = "sf_polygon")
Veamos todo junto en un mapa:
ggmap(mapa_mdp)+
geom_sf(data=mdp_polygon, fill=NA, size=1, color="firebrick3", inherit.aes = FALSE)+
labs(title="Mar del Plata",
caption="Fuente: Open Street Map")+
theme_void()
## Coordinate system already present. Adding new coordinate system, which will replace the existing one.
Listo ya tenemos nuestro área de interés.
Ahora con la función opq()
asignemos el bounding box correspondiente y con add_osm_feature()
especifiquemos que datos queremos, en este caso probemos con las calles (líneas) de la Ciudad de Mar del Plata, que en OSM se las denomina highway:
mdp_calles <- opq(bbox_mdp) %>%
add_osm_feature(key = "highway")
Hasta acá solo tenemos armada la búsqueda que vamos a realizar en OSM pero aún no hemos descargado los datos que necesitamos ya que para eso debemos usar osmdata_sf()
(puede tardar algunos minutos):
mdp_calles <- osmdata_sf(mdp_calles)
Veamos que nos trajimos:
mdp_calles
## Object of class 'osmdata' with:
## $bbox : -38.1273304,-57.6591066,-37.9132456,-57.5178445
## $overpass_call : The call submitted to the overpass API
## $meta : metadata including timestamp and version numbers
## $osm_points : 'sf' Simple Features Collection with 31171 points
## $osm_lines : 'sf' Simple Features Collection with 6001 linestrings
## $osm_polygons : 'sf' Simple Features Collection with 298 polygons
## $osm_multilines : NULL
## $osm_multipolygons : 'sf' Simple Features Collection with 14 multipolygons
Nos pudimos descargar toda la información relacionada a highway que se encontró dentro de los límites establecidos: 31171 puntos, 6001 lineas y 298 polígonos. Sin embargo ahora solo queremos quedarnos con osm_lines que representa a las calles:
mdp_calles <- mdp_calles$osm_lines
Ahora si ya tenemos nuestra información y podemos explorar el tamaño con dim()
:
dim(mdp_calles)
## [1] 6001 64
Efectivamente tenemos las 6001 calles junto a 64 columnas de información. Sumemos esto a nuestro mapa:
ggmap(mapa_mdp)+
geom_sf(data=mdp_polygon, fill=NA, size=1, color="firebrick3", inherit.aes = FALSE)+
geom_sf(data = mdp_calles, color="deepskyblue4", alpha=0.5, inherit.aes = FALSE)+
labs(title="Mar del Plata",
caption="Fuente: Open Street Map")+
theme_void()
## Coordinate system already present. Adding new coordinate system, which will replace the existing one.
Podrán ver que las calles se descargaron de acuerdo al límite del bounding box y no por los límites geográficos de Mar del Plata. Para solucionarlo alcanza con hacer una intersección entre ambos datos geográficos utilizando st_intersection()
:
mdp_calles <- st_intersection(mdp_calles, mdp_polygon)
## although coordinates are longitude/latitude, st_intersection assumes that they are planar
## Warning: attribute variables are assumed to be spatially constant throughout all
## geometries
ggmap(mapa_mdp)+
geom_sf(data=mdp_polygon, fill=NA, size=1, color="firebrick3", inherit.aes = FALSE)+
geom_sf(data = mdp_calles, color="deepskyblue4", alpha=0.5, inherit.aes = FALSE)+
labs(title="Mar del Plata",
caption="Fuente: Open Street Map")+
theme_void()
## Coordinate system already present. Adding new coordinate system, which will replace the existing one.
También podemos colorear las calles de acuerdo a algún atributo que tengan, como por ejemplo la velocidad máxima, la velocidad mínima, la cantidad de carriles o el tipo de calle. Probemos con la velocidad máxima (maxspeed):
ggmap(mapa_mdp)+
geom_sf(data = mdp_calles, aes(color=maxspeed), inherit.aes = FALSE)+
labs(title="Mar del Plata",
caption="Fuente: Open Street Map")+
theme_void()
## Coordinate system already present. Adding new coordinate system, which will replace the existing one.
Como verán, el campo maxspeed está en formato character a pesar de ser numérico ya que en los dataset que se descargan de OSM vienen todos los campos en formato texto. Vamos a tener que modificarlo:
mdp_calles <- mdp_calles %>%
mutate(maxspeed = as.numeric(maxspeed))
ggmap(mapa_mdp)+
geom_sf(data = mdp_calles, aes(color=maxspeed), inherit.aes = FALSE)+
scale_color_viridis_c()+
labs(title="Mar del Plata",
subtitle="Velocidad Máxima Permitida",
caption="Fuente: Open Street Map")+
theme_void()
## Coordinate system already present. Adding new coordinate system, which will replace the existing one.
Y ahora dejemos solo las Avenidas:
ggmap(mapa_mdp)+
geom_sf(data = filter(mdp_calles, str_detect(name, "Avenida")), color="red", inherit.aes = FALSE)+
scale_color_viridis_c()+
labs(title="Mar del Plata",
subtitle="Avenidas",
caption="Fuente: Open Street Map")+
theme_void()
## Coordinate system already present. Adding new coordinate system, which will replace the existing one.
Ahora probemos descargar datos con geometría de polígonos y aprovechando que estamos analizando Mar del Plata, que mejor que descargar las superficies de playas, ¿No?
mdp_playas <- opq(bbox_mdp) %>%
add_osm_feature(key = "natural", value = "beach")
mdp_playas <- osmdata_sf(mdp_playas)
mdp_playas
## Object of class 'osmdata' with:
## $bbox : -38.1273304,-57.6591066,-37.9132456,-57.5178445
## $overpass_call : The call submitted to the overpass API
## $meta : metadata including timestamp and version numbers
## $osm_points : 'sf' Simple Features Collection with 1287 points
## $osm_lines : 'sf' Simple Features Collection with 4 linestrings
## $osm_polygons : 'sf' Simple Features Collection with 43 polygons
## $osm_multilines : NULL
## $osm_multipolygons : 'sf' Simple Features Collection with 2 multipolygons
Tenemos 43 polígonos, usemos esos!
mdp_playas <- mdp_playas$osm_polygons
Y por las dudas recortemoslos por el polígono de Mar del Plata:
mdp_playas <- st_intersection(mdp_playas, mdp_polygon)
## although coordinates are longitude/latitude, st_intersection assumes that they are planar
## Warning: attribute variables are assumed to be spatially constant throughout all
## geometries
Veamos los polígonos en el mapa!
ggmap(mapa_mdp)+
geom_sf(data=mdp_playas, fill="darkorange", color=NA, inherit.aes = FALSE)+
labs(title="Mar del Plata",
subtitle="Ubicación de las Playas",
caption="Fuente: Open Street Map")+
theme_void()
## Coordinate system already present. Adding new coordinate system, which will replace the existing one.
Por último, descarguemos geometrías de tipo puntos, como por ejemplo los restaurant y bares de Mar del Plata:
mdp_restaurant <- opq(bbox_mdp) %>%
add_osm_feature(key = "amenity", value = c("restaurant", "bar"))
mdp_restaurant <- osmdata_sf(mdp_restaurant)
mdp_restaurant
## Object of class 'osmdata' with:
## $bbox : -38.1273304,-57.6591066,-37.9132456,-57.5178445
## $overpass_call : The call submitted to the overpass API
## $meta : metadata including timestamp and version numbers
## $osm_points : 'sf' Simple Features Collection with 284 points
## $osm_lines : NULL
## $osm_polygons : 'sf' Simple Features Collection with 14 polygons
## $osm_multilines : NULL
## $osm_multipolygons : NULL
Nos quedamos con los 284 puntos:
mdp_restaurant <- mdp_restaurant$osm_points
Y realizamos la intersección con los límites de Mar del Plata:
mdp_restaurant <- st_intersection(mdp_restaurant, mdp_polygon)
## although coordinates are longitude/latitude, st_intersection assumes that they are planar
## Warning: attribute variables are assumed to be spatially constant throughout all
## geometries
A mapear:
ggmap(mapa_mdp)+
geom_sf(data=mdp_restaurant, aes(color=amenity), inherit.aes = FALSE)+
labs(title="Mar del Plata",
subtitle="Bares y Restaurants",
color="Tipo",
caption="Fuente: Open Street Map")+
theme_void()
## Coordinate system already present. Adding new coordinate system, which will replace the existing one.
¿Y si queremos hacer zoom en la zona céntrica de Mar del Plata?
bbox_mdp_centro <- getbb('Centro, Mar del Plata, Buenos Aires, Argentina')
mapa_mdp_centro <- get_stamenmap(bbox = bbox_mdp_centro,
zoom=15)
ggmap(mapa_mdp_centro)
mdp_centro_polygon = getbb('Centro, Mar del Plata, Buenos Aires, Argentina', format_out = "sf_polygon")
mdp_restaurant_pm <- st_intersection(mdp_restaurant, mdp_centro_polygon)
## although coordinates are longitude/latitude, st_intersection assumes that they are planar
## Warning: attribute variables are assumed to be spatially constant throughout all
## geometries
ggmap(mapa_mdp_centro)+
geom_sf(data=filter(mdp_restaurant_pm, !is.na(amenity)), aes(color=amenity), size=3, inherit.aes = FALSE)+
labs(title="Mar del Plata - Zona Centro",
subtitle="Bares y Restaurants",
color="Tipo",
caption="Fuente: Open Street Map")+
scale_color_manual(values=c("deeppink", "darkgreen"))+
theme_void()
## Coordinate system already present. Adding new coordinate system, which will replace the existing one.
¿Y si queremos ubicar todos los bares y restaurantes en un mapa interactivo así no tenemos que hacer zooms estáticos (como el anterior) en todas las zonas de interés?
Bueno, para esto sirve leaflet()
y podemos utilizarlo de la siguiente forma:
leaflet(mdp_restaurant) %>%
addTiles() %>%
addMarkers()
Agreguemos información (tipo y nombre) en un pop up, pero antes filtremos aquellos que tienen valores nulos en dichos campos:
mdp_restaurant <- mdp_restaurant %>%
filter(!is.na(name), !is.na(amenity))
Ahora si:
factpal <- colorFactor(palette = c("seagreen4","deeppink4"),
levels = mdp_restaurant$amenity)
## Warning in colorFactor(palette = c("seagreen4", "deeppink4"), levels =
## mdp_restaurant$amenity): Duplicate levels detected
leaflet(mdp_restaurant) %>%
addTiles() %>%
addCircleMarkers(popup = paste("Tipo:", mdp_restaurant$amenity, "<br>",
"Nombre:", mdp_restaurant$name),
color = ~factpal(amenity))%>%
addLegend("bottomright", pal = factpal, values = ~amenity,
title = "Tipo",
opacity = 1)
Aprovechemos que tenemos la ubicación exacta de bares y restaurants de Mar del Plata, y armemos un mapa coroplético por barrio, como el que vimos la clase anterior. Para eso carguemos el geojson que podemos descargar de https://data.world/angie-scetta/mdp-barrios:
barrios_mdp <- st_read("barrios_mdp.geojson")
## Reading layer `barrios_mdp' from data source `C:\Users\27356214477\Documents\MEU-CDDPC2\barrios_mdp.geojson' using driver `GeoJSON'
## Simple feature collection with 124 features and 6 fields
## Geometry type: MULTIPOLYGON
## Dimension: XY
## Bounding box: xmin: -57.84755 ymin: -38.24135 xmax: -57.5225 ymax: -37.75424
## Geodetic CRS: WGS 84
Hagamos un mapa para visualizar ambos datos antes de unificarlos:
ggmap(mapa_mdp)+
geom_sf(data=barrios_mdp, color="firebrick3", fill=NA, inherit.aes = FALSE)+
geom_sf(data=mdp_restaurant, inherit.aes = FALSE)+
labs(title="Mar del Plata - Barrios",
subtitle="Bares y Restaurants",
caption="Fuente: Open Street Map")+
theme_void()
## Coordinate system already present. Adding new coordinate system, which will replace the existing one.
Y hagamos el join espacial:
mdp_restaurant <- st_join(mdp_restaurant, barrios_mdp)
Agrupemos por barrio:
mdp_restaurant <- mdp_restaurant %>%
group_by(barrio) %>%
summarise(cantidad=n()) %>%
st_set_geometry(NULL)
Veamos el top 5 de barrios con más bares y restaurantes:
mdp_restaurant %>%
filter(!is.na(barrio)) %>%
arrange(desc(cantidad)) %>%
head(5)
## # A tibble: 5 x 2
## barrio cantidad
## <chr> <int>
## 1 AREA CENTRO 23
## 2 LOS PINARES 10
## 3 PUNTA MOGOTES 9
## 4 ESTRADA JOSE MANUEL 8
## 5 SAN JOSE 5
El barrio que más bares y restaurantes tiene es el Centro, seguido por Los Troncos. Ahora agreguemosle la columna “cantidad” al dataset de barrios:
barrios_mdp <- barrios_mdp %>%
left_join(mdp_restaurant, by="barrio")
Finalmente hagamos nuestro mapa coroplético:
ggmap(mapa_mdp)+
geom_sf(data=barrios_mdp %>%
filter(!is.na(cantidad)), aes(fill=cantidad), alpha=0.75, color=NA, inherit.aes = FALSE)+
scale_fill_viridis_c(option="magma", direction=-1)+
labs(title="Mar del Plata - Barrios",
subtitle="Bares y Restaurants",
fill="Cantidad",
caption="Fuente: Open Street Map")+
theme_void()
## Coordinate system already present. Adding new coordinate system, which will replace the existing one.
Listo! Ahora les toca practicar a uds.
EJERCICIOS
Descargar de OpenStreetMap la grilla de calles para la Ciudad elegida en el ejercicio 1 y mapearla por uno de sus atributos (velocidad mínima, velocidad máxima, cantidad de carriles, etc).
Descargar de OpenStreetMap una (o más) capas de datos de tipo puntos o polígonos. (Ver catálogo de categorías en https://wiki.openstreetmap.org/wiki/Map_Features)