Week 6: APIS: Gis en R

1. Checklist

Antes de iniciar con esta lectura asegúrese de…

☑ Lectures previas

Puede darle una mirada a Introducción a GIS en R

☑ Versión de R

Tener la versión R version 4.1.1 (2021-08-10) instalada:

R.version.string

## [1] "R version 4.1.1 (2021-08-10)"

☑ Librerías

Instale/llame la librería pacman, y use la función p_load() para instalar/llamar las librerías de esta sesión:

## Llamar pacman (contiene la función p_load)
require(pacman) 

## Llama/instala-llama las librerías listadas
p_load(tidyverse,rio,
       sf, # Leer/escribir/manipular datos espaciales
       leaflet, # Visualizaciones dinámicas
       tmaptools, # geocode_OSM()
       osmdata) # Get OSM's data

2. OpenStreetMap

OpenStreetMap (OSM) es un proyecto de mapeo de acceso abierto global, gratuito y licenciado bajo ODbL Licence. OSM recopila información geográfica capturada con dispositivos GPS móviles, ortofotografías y otras fuentes de información libre.

2.1. Geocodificar direcciones

La función geocode_OSM() de la librería tmaptools se conecta a la API de OpenStreetMap y retorna el vértice con la coordenada geográfica del sitio/dirección buscado.

## Buscar un lugar público por el nombre
geocode_OSM("Casa de Nariño, Bogotá")

## $query
## [1] "Casa de Nariño, Bogotá"
## 
## $coords
##          x          y 
## -74.077549   4.595506 
## 
## $bbox
##       xmin       ymin       xmax       ymax 
## -74.077962   4.595103 -74.077078   4.595872

Puede adicionar el argumento as.sf=T para generar un objeto de clase sf:

## geocode_OSM no reconoce el caracter #, en su lugar se usa %23% 
point = geocode_OSM("Cra. 8 %23% 7-26, Bogotá", as.sf=T) 
point

## Simple feature collection with 1 feature and 7 fields
## Active geometry column: point
## Geometry type: POINT
## Dimension:     XY
## Bounding box:  xmin: -74.07907 ymin: 4.593874 xmax: -74.07907 ymax: 4.593874
## Geodetic CRS:  WGS 84
##                      query      lat       lon  lat_min  lat_max   lon_min
## 1 Cra. 8 %23% 7-26, Bogotá 4.593874 -74.07907 4.593004 4.594861 -74.07957
##     lon_max                           bbox                      point
## 1 -74.07849 POLYGON ((-74.07957 4.59300... POINT (-74.07907 4.593874)

Puede visualizar el objeto point usando la librería leaflet:

## la función addTiles adiciona la capa de OpenStreetMap
leaflet() %>% addTiles() %>% addCircles(data=point)

2.2. Librería `osmdata`

osmdata es una librería de R que permite descargar y usar datos de OpenStreetMap (OSM).

2.2.1. Features disponibles

Puede acceder a la lista de features disponibles en OSM aquí. En R puede obtener un vector con los nombres de los features usando la función available_features():

available_features() %>% head(20)

##  [1] "4wd_only"                "abandoned"              
##  [3] "abutters"                "access"                 
##  [5] "addr"                    "addr:city"              
##  [7] "addr:conscriptionnumber" "addr:country"           
##  [9] "addr:county"             "addr:district"          
## [11] "addr:flats"              "addr:full"              
## [13] "addr:hamlet"             "addr:housename"         
## [15] "addr:housenumber"        "addr:inclusion"         
## [17] "addr:interpolation"      "addr:place"             
## [19] "addr:postbox"            "addr:postcode"

Cada feature contiene una lista de tags. Puede acceder al vector de tags usando la función available_tags(). Por ejemplo, puede acceder a la lista de amenity así:

available_tags("amenity") %>% head(20)

##  [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"

2.2. Descargar features

Para descargar features desde OSM, primero debe definir un espacio geográfico y obtener la caja de de coordenadas que lo contiene:

## obtener la caja de coordenada que contiene el polígono de Bogotá
opq(bbox = getbb("Bogotá Colombia"))

## $bbox
## [1] "4.4711754,-74.2235137,4.8331695,-74.0102577"
## 
## $prefix
## [1] "[out:xml][timeout:25];\n(\n"
## 
## $suffix
## [1] ");\n(._;>;);\nout body;"
## 
## $features
## NULL
## 
## attr(,"class")
## [1] "list"           "overpass_query"
## attr(,"nodes_only")
## [1] FALSE

Puede almacenar un objeto osm de clase list y overpass_query con las estaciones de autobús para la ciudad de Bogotá:

## objeto osm
osm = opq(bbox = getbb("Bogotá Colombia")) %>%
      add_osm_feature(key="amenity" , value="bus_station") 
class(osm)

## [1] "list"           "overpass_query"

Puede acceder a los Simple Features del objeto osm usando la función osmdata_sf(). El objeto contiene una lista de objetos con los puntos, líneas y polígonos disponibles.

## extraer Simple Features Collection
osm_sf = osm %>% osmdata_sf()
osm_sf

## Object of class 'osmdata' with:
##                  $bbox : 4.4711754,-74.2235137,4.8331695,-74.0102577
##         $overpass_call : The call submitted to the overpass API
##                  $meta : metadata including timestamp and version numbers
##            $osm_points : 'sf' Simple Features Collection with 2656 points
##             $osm_lines : 'sf' Simple Features Collection with 98 linestrings
##          $osm_polygons : 'sf' Simple Features Collection with 370 polygons
##        $osm_multilines : NULL
##     $osm_multipolygons : 'sf' Simple Features Collection with 91 multipolygons

Puede acceder a los elementos de la lista y crear un objeto de clase sf:

## Obtener un objeto sf
bus_station = osm_sf$osm_points %>% select(osm_id,amenity) 
bus_station

## Simple feature collection with 2656 features and 2 fields
## Geometry type: POINT
## Dimension:     XY
## Bounding box:  xmin: -74.22302 ymin: 4.530743 xmax: -74.03537 ymax: 4.812363
## Geodetic CRS:  WGS 84
## First 10 features:
##              osm_id amenity                   geometry
## 261630505 261630505    <NA> POINT (-74.14693 4.631519)
## 261630506 261630506    <NA> POINT (-74.14694 4.631385)
## 261630508 261630508    <NA> POINT (-74.14423 4.631132)
## 261630509 261630509    <NA> POINT (-74.14422 4.631255)
## 266537186 266537186    <NA> POINT (-74.11613 4.654154)
## 266537187 266537187    <NA> POINT (-74.11531 4.652883)
## 266537189 266537189    <NA> POINT (-74.11439 4.652935)
## 266537456 266537456    <NA>  POINT (-74.11531 4.65398)
## 266537457 266537457    <NA>  POINT (-74.11507 4.65568)
## 292774696 292774696    <NA> POINT (-74.14562 4.631682)

Visualiza el objeto bus_station:

## Pintar las estaciones de autobus
leaflet() %>% addTiles() %>% addCircleMarkers(data=bus_station , col="red")

3. Operaciones geometrías

Puede acceder a las viñetas de la librería sf así:

## Help
vignette("sf3")
vignette("sf4")

3.1. Obtener datos

## apartamentos
housing = import("https://eduard-martinez.github.io/data/house_points.rds")
housing %>% head()

## Simple feature collection with 4 features and 1 field
## Geometry type: POINT
## Dimension:     XY
## Bounding box:  xmin: -74.04945 ymin: 4.67549 xmax: -74.04716 ymax: 4.67857
## Geodetic CRS:  WGS 84
## # A tibble: 4 × 2
##   cod_apto            geometry
##   <chr>            <POINT [°]>
## 1 Apto 101 (-74.04926 4.67737)
## 2 Apto 102 (-74.04944 4.67549)
## 3 Apto 103 (-74.04716 4.67678)
## 4 Apto 104 (-74.04945 4.67857)

## precios medianos de las viviendas
mnz = import("https://eduard-martinez.github.io/data/mnz_prices_house.rds")
mnz %>% head()

## Simple feature collection with 6 features and 2 fields
## Geometry type: POLYGON
## Dimension:     XY
## Bounding box:  xmin: -74.05689 ymin: 4.67863 xmax: -74.04354 ymax: 4.682574
## Geodetic CRS:  WGS 84
##   MANZ_CAG median_price                       geometry
## 2   144838    1072.8725 POLYGON ((-74.04457 4.68209...
## 3   145114     879.1825 POLYGON ((-74.04434 4.68064...
## 4   145134     911.0521 POLYGON ((-74.05272 4.68059...
## 5   144834     805.2859 POLYGON ((-74.05675 4.68239...
## 7   144808    1230.4512 POLYGON ((-74.04722 4.68257...
## 8   144814     919.7815 POLYGON ((-74.04799 4.68249...

## bares
bar = opq(bbox = st_bbox(mnz)) %>%
      add_osm_feature(key = "amenity", value = "bar") %>%
      osmdata_sf() %>% .$osm_points %>% select(osm_id,name)
bar %>% head()

## Simple feature collection with 6 features and 2 fields
## Geometry type: POINT
## Dimension:     XY
## Bounding box:  xmin: -74.05089 ymin: 4.675395 xmax: -74.04696 ymax: 4.679601
## Geodetic CRS:  WGS 84
##                osm_id                  name                   geometry
## 595803379   595803379              La Vasca POINT (-74.05089 4.676612)
## 1484781209 1484781209 Galería, Café y Libro POINT (-74.04795 4.675395)
## 1661484725 1661484725             Full 80´S POINT (-74.04759 4.679601)
## 4214169507 4214169507                  <NA> POINT (-74.04718 4.676927)
## 4214169509 4214169509                  <NA> POINT (-74.04696 4.676819)
## 4214169511 4214169511                  <NA> POINT (-74.04704 4.676671)

## parque de la 93
park = getbb(place_name = "Parque de la 93", 
             featuretype = "amenity",
             format_out = "sf_polygon")

## old-style crs object detected; please recreate object with a recent sf::st_crs()
## old-style crs object detected; please recreate object with a recent sf::st_crs()
## old-style crs object detected; please recreate object with a recent sf::st_crs()
## old-style crs object detected; please recreate object with a recent sf::st_crs()

park %>% head()

## Simple feature collection with 1 feature and 0 fields
## Geometry type: POLYGON
## Dimension:     XY
## Bounding box:  xmin: -74.04913 ymin: 4.676096 xmax: -74.04744 ymax: 4.67744
## Geodetic CRS:  WGS 84
##                         geometry
## 1 POLYGON ((-74.04913 4.67681...

Visualizar todos los objetos:

leaflet() %>% addTiles() %>% 
addPolygons(data=mnz) %>% # manzanas
addPolygons(data=park , col="green") %>%  # parque de la 93
addCircles(data=housing , col="red" , weight=2) %>% # apartamentos
addCircles(data=bar , col="black" , weight=2) # bares

3.2. Precio de las viviendas

Puede adicionar el precio mediano de las viviendas (mnz) a cada vivienda (housing):

## Afinar las transformaciones
st_crs(mnz) == st_crs(housing)

## [1] TRUE

## unir dos conjuntos de datos basados en la geometría
housing_p = st_join(x=housing , y=mnz)
housing_p

## Simple feature collection with 4 features and 3 fields
## Geometry type: POINT
## Dimension:     XY
## Bounding box:  xmin: -74.04945 ymin: 4.67549 xmax: -74.04716 ymax: 4.67857
## Geodetic CRS:  WGS 84
## # A tibble: 4 × 4
##   cod_apto            geometry MANZ_CAG median_price
## * <chr>            <POINT [°]> <chr>           <dbl>
## 1 Apto 101 (-74.04926 4.67737) 145431          1041.
## 2 Apto 102 (-74.04944 4.67549) 145651           925.
## 3 Apto 103 (-74.04716 4.67678) 145555           874.
## 4 Apto 104 (-74.04945 4.67857) 145333          1054.

Puede adicionar el precio máximo de las manzanas vecinas a cada vivienda:

## precio promedio
max_price = st_join(x=st_buffer(x=housing , dist=50) , y=mnz)
st_geometry(max_price) = NULL

## precio máximo
max_price = max_price %>% group_by(cod_apto) %>% summarise(price_max=max(median_price))
max_price

## # A tibble: 4 × 2
##   cod_apto price_max
##   <chr>        <dbl>
## 1 Apto 101     1054.
## 2 Apto 102     1051.
## 3 Apto 103     1013.
## 4 Apto 104     1054.

## join data
housing_p = left_join(x=housing_p , y=max_price , by="cod_apto")
housing_p

## Simple feature collection with 4 features and 4 fields
## Geometry type: POINT
## Dimension:     XY
## Bounding box:  xmin: -74.04945 ymin: 4.67549 xmax: -74.04716 ymax: 4.67857
## Geodetic CRS:  WGS 84
## # A tibble: 4 × 5
##   cod_apto            geometry MANZ_CAG median_price price_max
##   <chr>            <POINT [°]> <chr>           <dbl>     <dbl>
## 1 Apto 101 (-74.04926 4.67737) 145431          1041.     1054.
## 2 Apto 102 (-74.04944 4.67549) 145651           925.     1051.
## 3 Apto 103 (-74.04716 4.67678) 145555           874.     1013.
## 4 Apto 104 (-74.04945 4.67857) 145333          1054.     1054.

3.3. Distancia a amenities

Distancia al parque de la 93:

## Distancia a parque
housing_p$dist_park = st_distance(x=housing_p , y=park)
housing_p

## Simple feature collection with 4 features and 5 fields
## Geometry type: POINT
## Dimension:     XY
## Bounding box:  xmin: -74.04945 ymin: 4.67549 xmax: -74.04716 ymax: 4.67857
## Geodetic CRS:  WGS 84
## # A tibble: 4 × 6
##   cod_apto            geometry MANZ_CAG median_price price_max
## * <chr>            <POINT [°]> <chr>           <dbl>     <dbl>
## 1 Apto 101 (-74.04926 4.67737) 145431          1041.     1054.
## 2 Apto 102 (-74.04944 4.67549) 145651           925.     1051.
## 3 Apto 103 (-74.04716 4.67678) 145555           874.     1013.
## 4 Apto 104 (-74.04945 4.67857) 145333          1054.     1054.
## # … with 1 more variable: dist_park [m]

Matrix de distancias a todos los bares:

## Distancia a bares
dist_bar = st_distance(x=housing_p , y=bar)
dist_bar

## Units: [m]
##          [,1]     [,2]     [,3]      [,4]      [,5]      [,6]      [,7]
## [1,] 199.0130 263.3124 309.4507 235.96082 261.77220 258.17845 231.95678
## [2,] 203.1238 165.5768 500.9486 297.29290 311.72986 296.77635 281.56409
## [3,] 413.4238 177.0878 317.2686  16.45325  22.21122  18.12409  10.29698
## [4,] 269.7492 390.2617 235.7994 311.11455 337.39084 340.61366 314.59799
##          [,8]      [,9]    [,10]    [,11]     [,12]     [,13]
## [1,] 1423.976 1030.5068 748.6555 139.0405 1007.1293  76.17941
## [2,] 1330.587  867.8151 588.9216 117.0754  829.8639 274.29384
## [3,] 1199.550 1153.9565 877.4660 181.3302 1108.8251 217.59272
## [4,] 1521.987 1122.4758 842.8956 271.0421 1109.0516 106.75646

Puede calcular la distancia mínima a cada bar:

## Distancia minima
min_dist = apply(dist_bar , 1 , min)
min_dist

## [1]  76.17941 117.07540  10.29698 106.75646

housing_p$dist_bar = min_dist
housing_p

## Simple feature collection with 4 features and 6 fields
## Geometry type: POINT
## Dimension:     XY
## Bounding box:  xmin: -74.04945 ymin: 4.67549 xmax: -74.04716 ymax: 4.67857
## Geodetic CRS:  WGS 84
## # A tibble: 4 × 7
##   cod_apto            geometry MANZ_CAG median_price price_max
## * <chr>            <POINT [°]> <chr>           <dbl>     <dbl>
## 1 Apto 101 (-74.04926 4.67737) 145431          1041.     1054.
## 2 Apto 102 (-74.04944 4.67549) 145651           925.     1051.
## 3 Apto 103 (-74.04716 4.67678) 145555           874.     1013.
## 4 Apto 104 (-74.04945 4.67857) 145333          1054.     1054.
## # … with 2 more variables: dist_park [m], dist_bar <dbl>

Referencias

Lovelace, R., Nowosad, J., & Muenchow, J. (2019). Geocomputation with R. [Ver aquí]
- Cap. 4: Spatial Data Operations
- Cap. 5: Geometry Operations
- Cap. 6: Reprojecting geographic data
- Cap. 11: Statistical learning
Bivand, R. S., Pebesma, E. J., Gómez-Rubio, V., & Pebesma, E. J. (2013). Applied spatial data analysis with R. [Ver aquí]
- Cap. 4: Spatial Data Import and Export
- Cap. 7: Spatial Point Pattern Analysis
- Cap. 8: Interpolation and Geostatistics

Week 6: Obtener datos GIS en R

Big Data and Machine Learning for Applied Economics | MECA-4107

Eduard-Martínez & Lucas Antonio II

1. Checklist

☑ Lectures previas

☑ Versión de R

☑ Librerías

2. OpenStreetMap

2.1. Geocodificar direcciones

2.2. Librería `osmdata`

2.2.1. Features disponibles

2.2. Descargar features

3. Operaciones geometrías

3.1. Obtener datos

3.2. Precio de las viviendas

3.3. Distancia a amenities

Referencias

Week 6: Obtener datos GIS en R

Big Data and Machine Learning for Applied Economics | MECA-4107

Eduard-Martínez & Lucas Antonio II

1. Checklist

☑ Lectures previas

☑ Versión de R

☑ Librerías

2. OpenStreetMap

2.1. Geocodificar direcciones

2.2. Librería osmdata

2.2.1. Features disponibles

2.2. Descargar features

3. Operaciones geometrías

3.1. Obtener datos

3.2. Precio de las viviendas

3.3. Distancia a amenities

Referencias

2.2. Librería `osmdata`