hugo-ox-mpvd/publicaciones.org
2023-03-28 11:47:37 +02:00

64 KiB
Raw Permalink Blame History

Artículos

PUBLICADO Artículos

subtitle = "y programación literaria"
view = 2
slug = "publicacion"
pagecat = "publicaciones"
description = "Listado de publicaciones"
url = "/publicacion/"
layout = "publication"
[header]
focal_point = "Right"

PUBLICADO Streetmaps

  [image]
   caption = "[Orion aveugle cherchant le soleil](https://commons.wikimedia.org/wiki/File:Orion_aveugle_cherchant_le_soleil.jpg), Nicolas Poussin  / Public domain"
   focal_point = "BottomLeft"
   placement = 3
   preview_only = false

El título original es mucho más feo ya que incluye al final "y PowerPoint": "Create a streetmap of your favorite city with ggplot2 and powerpoint". En este caso vamos a intentar no usar esa herramienta por no ser software libre.

El propósito del autor era crear mapas/callejeros como el de este cartel:

Para este tutorial se utilizan dos paquetes de R: osmdata y ggplot2. Con osmdata se extraen las calles de OpenStreetMap, una base de datos libre con licencia abierta. Ni siquiera requiere de una clave para trabajar con su API. Para crear las visualizaciones y trabajar con los datos se utiliza tidyverse.

Típica imagen de cartel

El autor no había trabajado mucho con ggplot2 y se inspiró especialmente en Dominic Royé y su detallado trabajo de datos OpenStreetMap con R.

Se puede empezar por la ciudad en la que estás, en la que vives, en la que naciste o una que te gusta. En su caso empezó por Freiburg, en el sur de Alemania. Para elegir la ciudad tan solo hay que ajustar la longitud y latitud de los puntos de inicio.

Instalar librerías

Lo primero es instalar las librerías necesarias, tidyverse y osmdata:

library(tidyverse)
library(osmdata)

OpenStreetMap

Antes de trabajar con R hay que entender cómo almacena OpenStreetMap los datos de los mapas y, por tanto, las calles.

OpenStreetMap describe las cosas físicas como "características" o "rasgos", del inglés features. En este caso son map_features Estas características se almacenan como pares de clave-valor.

Por ejemplo, hay highways que son vías principales que pueden conectar tanto la ciudad consigo misma como una ciudad con otras.

Para obtener todas las etiquetas de una característica espacial vía osmdata introducimos la función available_tags()

available_tags("highway")
 [1] "bridleway"              "bus_guideway"           "bus_stop"              
 [4] "construction"           "corridor"               "crossing"              
 [7] "cycleway"               "elevator"               "emergency_access_point"
[10] "emergency_bay"          "escape"                 "footway"               
[13] "give_way"               "living_street"          "milestone"             
[16] "mini_roundabout"        "motorway"               "motorway_junction"     
[19] "motorway_link"          "passing_place"          "path"                  
[22] "pedestrian"             "platform"               "primary"               
[25] "primary_link"           "proposed"               "raceway"               
[28] "residential"            "rest_area"              "road"                  
[31] "secondary"              "secondary_link"         "service"               
[34] "services"               "speed_camera"           "steps"                 
[37] "stop"                   "street_lamp"            "tertiary"              
[40] "tertiary_link"          "toll_gantry"            "track"                 
[43] "traffic_mirror"         "traffic_signals"        "trailhead"             
[46] "trunk"                  "trunk_link"             "turning_circle"        
[49] "turning_loop"           "unclassified"

También hay vías de agua waterway, ya sean ríos, canales, acequias o tuberías. Se hace la misma consulta:

available_tags("waterway")
 [1] "boatyard"      "canal"         "dam"           "ditch"        
 [5] "dock"          "drain"         "fairway"       "fuel"         
 [9] "lock_gate"     "pressurised"   "river"         "riverbank"    
[13] "stream"        "tidal_channel" "turning_point" "water_point"  
[17] "waterfall"     "weir"
Obtener taxonomía de map_features

Se pueden consultar las taxonomías en su página web de map_features o bien con la función available_features():

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_scale"              
[139] "mtb:description"         "mtb:scale:imba"         
[141] "name"                    "narrow"                 
[143] "natural"                 "noexit"                 
[145] "non_existent_levels"     "note"                   
[147] "nudism"                  "office"                 
[149] "official_name"           "old_name"               
[151] "oneway"                  "opening_hours"          
[153] "operator"                "organic"                
[155] "oven"                    "overtaking"             
[157] "parking:condition"       "parking:lane"           
[159] "passing_places"          "place"                  
[161] "power"                   "priority_road"          
[163] "produce"                 "proposed"               
[165] "protected_area"          "psv"                    
[167] "public_transport"        "railway"                
[169] "railway:preserved"       "railway:track_ref"      
[171] "recycling_type"          "ref"                    
[173] "religion"                "residential"            
[175] "resource"                "roadtrain"              
[177] "route"                   "sac_scale"              
[179] "service"                 "service_times"          
[181] "shelter_type"            "shop"                   
[183] "sidewalk"                "site"                   
[185] "ski"                     "smoothness"             
[187] "social_facility"         "speed_pedelec"          
[189] "start_date"              "step_count"             
[191] "substation"              "surface"                
[193] "tactile_paving"          "tank"                   
[195] "tidal"                   "toilets:wheelchair"     
[197] "toll"                    "tourism"                
[199] "tracks"                  "tracktype"              
[201] "traffic_calming"         "traffic_sign"           
[203] "trail_visibility"        "tunnel"                 
[205] "turn"                    "type"                   
[207] "usage"                   "vehicle"                
[209] "vending"                 "voltage"                
[211] "water"                   "wheelchair"             
[213] "wholesale"               "width"                  
[215] "winter_road"             "wood"
Obtener coordenadas

Para obtener las coordenadas de la ciudad que elijamos se puede usar la función getbb(). Por ejemplo, con Murcia:

getbb("Murcia Spain")
        min        max
x -1.384834 -0.8508899
y 37.715872 38.1179185

Esto nos da las coordenadas de la ciudad. el valor de la x da la longitud mientras que la y da la altitud.

Exportar datos de OSM

Ahora queremos exportar las carreteras del sistema de coordenadas. Para ello primero pasamos la salida de la función getbb() a la función opq().

A continuación pasamos esta salida a la función add_osm_feature(). La función tiene dos argumentos. Con la clave especificamos la clave de la característica; con el valor especificamos la etiqueta de la característica.

En este caso primero extraemos las principales calles de la ciudad y las pasamos a la función osmdata_sf para insertarlo luego en ggplot2.

  calles_principales <- getbb("Murcia Spain")%>%
    opq()%>%
    add_osm_feature(key = "highway", 
			value = c("motorway", "primary", 
				  "secondary", "tertiary")) %>%
    osmdata_sf()
  calles_principales
Object of class 'osmdata' with:
                 $bbox : 37.7158719,-1.3848335,38.1179185,-0.8508899
        $overpass_call : The call submitted to the overpass API
                 $meta : metadata including timestamp and version numbers
           $osm_points : 'sf' Simple Features Collection with 39961 points
            $osm_lines : 'sf' Simple Features Collection with 4748 linestrings
         $osm_polygons : 'sf' Simple Features Collection with 258 polygons
       $osm_multilines : NULL
    $osm_multipolygons : NULL

Los datos los hemos almacenado en las variables calles_principales, peatonales y agua Estos objeto tienen diferentes objetos hijxs. Los que nos interesan para este mapa son osm_lines ya que son líneas.

Completamos objetos

Ahora añadimos también las calles más pequeñas, peatonales y cursos de agua:

  calles_principales <- getbb("Murcia Spain")%>%
    opq()%>%
    add_osm_feature(key = "highway", 
			  value = c("motorway", "primary", "secondary", "tertiary")) %>%
    osmdata_sf()

  agua <- getbb("Murcia Spain")%>%
    opq()%>%
    add_osm_feature(key = "waterway",
		      value = c("river", "canal", "drain", "stream", "dock", "lock_gate", "boatyard", "riverbank", "stream")) %>%
    osmdata_sf()

  bici <- getbb("Murcia Spain")%>%
    opq()%>%
    add_osm_feature(key = "highway",
			  value = "cycleway") %>%
    osmdata_sf()
  calles_secundarias <- getbb("Murcia Spain")%>%
    opq()%>%
    add_osm_feature(key = "highway",
		      value = c("residential", "living_street", "unclassified", "service")) %>%
    osmdata_sf()

Daba error y he quitado las peatonales:

  peatonales <- getbb("Murcia Spain")%>%
    opq()%>%
    add_osm_feature(key = "highway",
			  value = c("footway", "cycleway")) %>%
    osmdata_sf()

Y las he convertido en bici.

Primer mapa

Creamos nuestro primer mapa con la función ggplot() y nuestro conjunto de datos agua:

ggplot() +
  geom_sf(data = agua$osm_lines,
          inherit.aes = FALSE,
          color = "blue",
          size = 1.4,
          alpha = .8) +
  coord_sf(xlim = c(-1.15, -1.10), 
           ylim = c(37.96, 38.01),
           expand = FALSE)

/flow/hugo-ox-mpvd/src/branch/main/murcia-agua-uno.pdf

ggplot() +
  geom_sf(data = calles_secundarias$osm_lines,
          inherit.aes = FALSE,
          color = "black",
          size = .4,
          alpha = .8) +
  coord_sf(xlim = c(-1.15, -1.10), 
           ylim = c(37.96, 38.01),
           expand = FALSE)

/flow/hugo-ox-mpvd/src/branch/main/murcia-uno.pdf

Primero se añade la geometría geom_sf a la función ggplot. Para los datos añadimos las calles almacenadas en la variable streets$osm_lines. Se puede determinar la anchura de las calles con size. Para que las calles no estén completamente negras se ha creado una pequeña transparencia con el nivel alpha.

Con la función coord_sf se puede determinar el eje de las X y de las Y exactamente. Es mejor jugar con los valores hasta que has definido los límites.

Con expand = FALSE nos aseguramos de que se muestren las coordenadas correctamente.

Segundo mapa

Ahora se añaden calles pequeñas y ríos:

  ggplot() +
    geom_sf(data = calles_principales$osm_lines,
	    inherit.aes = FALSE,
	    color = "yellow",
	    size = .4,
	    alpha = .8) +
    geom_sf(data = calles_secundarias$osm_lines,
	    inherit.aes = FALSE,
	    color = "grey",
	    size = .4,
	    alpha = .6) +
    geom_sf(data = bici$osm_lines,
	    inherit.aes = FALSE,
	    color = "green",
	    size = .5,
	    alpha = .6) +
    geom_sf(data = agua$osm_lines,
	    inherit.aes = FALSE,
	    color = "blue",
	    size = .2,
	    alpha = .5) +
    coord_sf(xlim = c(-1.15, -1.10), 
	     ylim = c(37.96, 38.01),
	     expand = FALSE)

/flow/hugo-ox-mpvd/src/branch/main/murcia-dos.pdf

Colores (para borrar)

En vez de negro se pueden resaltar las calles en color:

ggplot() +
  geom_sf(data = streets$osm_lines,
          inherit.aes = FALSE,
          color = "steelblue",
          size = .4,
          alpha = .8) +
  geom_sf(data = small_streets$osm_lines,
          inherit.aes = FALSE,
          color = "black",
          size = .4,
          alpha = .6) +
  geom_sf(data = river$osm_lines,
          inherit.aes = FALSE,
          color = "black",
          size = .2,
          alpha = .5) +
  coord_sf(xlim = c(-1.38, -0.85), 
           ylim = c(37.71, 38.11),
           expand = FALSE)

/flow/hugo-ox-mpvd/src/branch/main/boceto-tres.pdf

Sin rejilla

Con theme_void() borramos la rejilla de los ejes de coordenadas:

  ggplot() +
    geom_sf(data = calles_principales$osm_lines,
	    inherit.aes = FALSE,
	    color = "yellow",
	    size = .4,
	    alpha = .8) +
    geom_sf(data = calles_secundarias$osm_lines,
	    inherit.aes = FALSE,
	    color = "grey",
	    size = .4,
	    alpha = .6) +
    geom_sf(data = bici$osm_lines,
	    inherit.aes = FALSE,
	    color = "green",
	    size = .5,
	    alpha = .6) +
    geom_sf(data = agua$osm_lines,
	    inherit.aes = FALSE,
	    color = "blue",
	    size = .2,
	    alpha = .5) +
    coord_sf(xlim = c(-1.15, -1.10), 
	     ylim = c(37.96, 38.01),
	     expand = FALSE) +
    theme_void()

/flow/hugo-ox-mpvd/src/branch/main/murcia-sin-rejilla.pdf

ggplot() +
  geom_sf(data = streets$osm_lines,
          inherit.aes = FALSE,
          color = "steelblue",
          size = .4,
          alpha = .8) +
  geom_sf(data = small_streets$osm_lines,
          inherit.aes = FALSE,
          color = "black",
          size = .4,
          alpha = .6) +
  geom_sf(data = river$osm_lines,
          inherit.aes = FALSE,
          color = "black",
          size = .2,
          alpha = .5) +
  coord_sf(xlim = c(-1.38, -0.85), 
           ylim = c(37.71, 38.11),
           expand = FALSE) +
  theme_void()

/flow/hugo-ox-mpvd/src/branch/main/boceto-cuatro.pdf

Colores invertidos

También se pueden ajustar los colores haciendo un invert, es decir, un fondo negro sobre:

    ggplot() +
      geom_sf(data = calles_principales$osm_lines,
	      inherit.aes = FALSE,
	      color = "#ffbe7f",
	      size = .4,
	      alpha = .8) +
      geom_sf(data = calles_secundarias$osm_lines,
	      inherit.aes = FALSE,
	      color = "grey",
	      size = .4,
	      alpha = .6) +
      geom_sf(data = bici$osm_lines,
	      inherit.aes = FALSE,
	      color = "green",
	      size = .5,
	      alpha = .2) +
      geom_sf(data = agua$osm_lines,
	      inherit.aes = FALSE,
	      color = "#7fc0ff",
	      size = .4,
	      alpha = .5) +
      coord_sf(xlim = c(-1.15, -1.10), 
	       ylim = c(37.96, 38.01),
	       expand = FALSE) +
      theme_void() +
    theme(
      plot.background = element_rect(fill = "black")
      )

/flow/hugo-ox-mpvd/src/branch/main/murcia-invertido.pdf

Grande
    ggplot() +
      geom_sf(data = calles_principales$osm_lines,
	      inherit.aes = FALSE,
	      color = "#ffbe7f",
	      size = .4,
	      alpha = .8) +
      geom_sf(data = calles_secundarias$osm_lines,
	      inherit.aes = FALSE,
	      color = "grey",
	      size = .4,
	      alpha = .6) +
      geom_sf(data = bici$osm_lines,
	      inherit.aes = FALSE,
	      color = "green",
	      size = .5,
	      alpha = .2) +
      geom_sf(data = agua$osm_lines,
	      inherit.aes = FALSE,
	      color = "#7fc0ff",
	      size = .4,
	      alpha = .5) +
      coord_sf(xlim = c(-1.21, -0.96), 
	       ylim = c(37.90, 38.03),
	       expand = FALSE) +
      theme_void() +
    theme(
      plot.background = element_rect(fill = "black")
      )

/flow/hugo-ox-mpvd/src/branch/main/murcia-medio-invertido.pdf

Mayor
    ggplot() +
      geom_sf(data = calles_principales$osm_lines,
	      inherit.aes = FALSE,
	      color = "#ffbe7f",
	      size = .4,
	      alpha = .8) +
      geom_sf(data = calles_secundarias$osm_lines,
	      inherit.aes = FALSE,
	      color = "grey",
	      size = .4,
	      alpha = .6) +
      geom_sf(data = bici$osm_lines,
	      inherit.aes = FALSE,
	      color = "green",
	      size = .5,
	      alpha = .2) +
      geom_sf(data = agua$osm_lines,
	      inherit.aes = FALSE,
	      color = "#7fc0ff",
	      size = .4,
	      alpha = .5) +
      coord_sf(xlim = c(-1.30, -1.10), 
	       ylim = c(37.81, 38.08),
	       expand = FALSE) +
      theme_void() +
    theme(
      plot.background = element_rect(fill = "black")
      )

/flow/hugo-ox-mpvd/src/branch/main/murcia-mayor-invertido.pdf

Más grande
    ggplot() +
      geom_sf(data = calles_principales$osm_lines,
	      inherit.aes = FALSE,
	      color = "#ffbe7f",
	      size = .4,
	      alpha = .8) +
      geom_sf(data = calles_secundarias$osm_lines,
	      inherit.aes = FALSE,
	      color = "grey",
	      size = .4,
	      alpha = .6) +
      geom_sf(data = bici$osm_lines,
	      inherit.aes = FALSE,
	      color = "green",
	      size = .5,
	      alpha = .2) +
      geom_sf(data = agua$osm_lines,
	      inherit.aes = FALSE,
	      color = "#7fc0ff",
	      size = .4,
	      alpha = .5) +
      coord_sf(xlim = c(-1.38, -0.85), 
	       ylim = c(37.71, 38.11),
	       expand = FALSE) +
      theme_void() +
    theme(
      plot.background = element_rect(fill = "black")
      )

/flow/hugo-ox-mpvd/src/branch/main/murcia-biggest-invertido.pdf

x -1.384834 -0.8508899
y 37.715872 38.1179185
Antiguo
  ggplot() +
    geom_sf(data = streets$osm_lines,
	    inherit.aes = FALSE,
	    color = "#ffbe7f",
	    size = .4,
	    alpha = .8) +
    geom_sf(data = small_streets$osm_lines,
	    inherit.aes = FALSE,
	    color = "grey",
	    size = .4,
	    alpha = .6) +
    geom_sf(data = river$osm_lines,
	    inherit.aes = FALSE,
	    color = "#7fc0ff",
	    size = 1.2,
	    alpha = .5) +
    coord_sf(xlim = c(-1.15, -1.10), 
	     ylim = c(37.96, 38.00),
	     expand = FALSE)  +
    theme_void() +
  theme(
    plot.background = element_rect(fill = "#f9f9f9")
  )

/flow/hugo-ox-mpvd/src/branch/main/boceto-cinco.pdf

  • Esquina superior izqda.: 38.0001/-1.1576
  • Esquina superior dcha.: 38.0026/-1.1059
  • Esquina inferior izqda.: 37.9657/-1.1541
  • Esquina inferior dcha.: 37.9675/-1.1058

Postproducción

Una vez que están hechas las visualizaciones toca editarlo. En el caso del manual en el que se basa este manual la imagen no se guardaba sino que lo generaba R al vuelo.

En nuestro caso he ido generando los distintos bocetos con Emacs + Babel + R y guardándolos en una carpeta en PDF.

Para generar la imagen finalmente usaba la función ggmap():

ggsave("map.png", width = 6, height = 6)

Nosotrxs, en vez de crear un documento en Powerpoint y redimensionar el mapa hemos abierto un documento en Inkscape donde se edita gráficamente con el estilo de los otros.

Datos de OpenStreetMap con R

  [image]
   caption = "[Orion aveugle cherchant le soleil](https://commons.wikimedia.org/wiki/File:Orion_aveugle_cherchant_le_soleil.jpg), Nicolas Poussin  / Public domain"
   focal_point = "BottomLeft"
   placement = 3
   preview_only = false

Cuando Christian Burkhart creó su mapa del callejero con R se fijó mucho en el trabajo de Dominique Royé con las librerías que trabajan con datos de OpenStreetMap como el de las estaciones de servicio de Europa. Para ello utilizó POI (Point of Interest, puntos de interés en la terminología de las representaciones espaciales de datos, puntos geográficos al fin y al cabo).

Para obtener los datos se utiliza una pasarela de la API.

Instalación de paquetes

Lo primer paso es instalar las librerías necesarias como las librerías tidyverse, que es una colección de distintas librerías entre las que se incluye dplyr para la manipulación de datos o ggplot2 para la visualización; sf que es el nuevo estándar para trabajar con datos espaciales y es compatible con ggplot y dplyr. Finalmente ggmap nos ayuda a crear mapas.

  if(!require("osmdata")) install.packages("osmdata")
  if(!require("tidyverse")) install.packages("tidyverse")
  if(!require("ggmap")) install.packages("ggmap")
  if(!require("sf")) install.packages("sf")
Loading required package: sf
Linking to GEOS 3.9.0, GDAL 3.2.1, PROJ 7.2.1

Se activan las librerías

Y se activan las librerías:

library(tidyverse)
library(osmdata)
library(sf)
library(ggmap)

Hacer una consulta

Antes de hacer una consulta hay que saber lo que se puede filtrar. Para ello se puede probar con la función available_features() que devuelve las características espaciales o features de OSM, algo común en el idioma de la representación espacial geográfica.

Para no cargar mucho la consulta incluimos la función en la función head():

head(available_features())
[1] "4wd_only"  "abandoned" "abutters"  "access"    "addr"      "addr:city"

Estas son las cinco primeras features que devueve la consulta.

Se pueden ver las etiquetas de cualquiera de ellas a través de la función available_tags():

head(available_tags("tourism"))
[1] "alpine_hut" "apartment"  "aquarium"   "artwork"    "attraction"
[6] "camp_pitch"

La primera consulta: dónde están las peluquerías en Murcia

Para construir una consulta se usa el operador %>% que ayuda a encadenar varias funciones sin asignar el resultado a un nuevo objeto. Su uso está muy extendido entre la colección de librerías de tidyverse. Este tutorial de datacamp lo explica más.

En la primera parte de la consulta tenemos que indicar el lugar para el que queremos extraer información. Eso se consigue con la función getbb() que crea unos límites para un lugar según un nombre dado.

La función principal es opq() que crea la consulta final.

Se filtra la información que se busca con la función add_osm_feature().

En esta primera consulta se buscan los cines de Madrid, con amenity como feature y cinema como tag. Se pueden obtener los resultados espaciales en varios formatos. la función osmdata_*() envía la consulta al servidor y según si el sufijo es sf, sp o xml devuelve una característica simple (de single feature), espacial (de spatial) o en formato XML.

consulta <- getbb("Murcia Spain") %>%
  opq() %>%
   add_osm_feature("shop", "hairdresser")
str(consulta)
List of 4
 $ bbox    : chr "37.7158719,-1.3848335,38.1179185,-0.8508899"
 $ prefix  : chr "[out:xml][timeout:25];\n(\n"
 $ suffix  : chr ");\n(._;>;);\nout body;"
 $ features: chr " [\"shop\"=\"hairdresser\"]"
 - attr(*, "class")= chr [1:2] "list" "overpass_query"
 - attr(*, "nodes_only")= logi FALSE
pelus <- osmdata_sf(consulta)
pelus
Object of class 'osmdata' with:
                 $bbox : 37.7158719,-1.3848335,38.1179185,-0.8508899
        $overpass_call : The call submitted to the overpass API
                 $meta : metadata including timestamp and version numbers
           $osm_points : 'sf' Simple Features Collection with 74 points
            $osm_lines : NULL
         $osm_polygons : 'sf' Simple Features Collection with 1 polygons
       $osm_multilines : NULL
    $osm_multipolygons : NULL

El resultado es una lista de diferentes objetos espaciales. En nuestro caso solo nos interesan los puntos osm_points.

Visualizar

Lo bueno de los objetos sf es que para ggplot ya existe una función geométrica llamada geom_sf(). De fondo se puede añadir un mapa con la librería ggmap y la función get_map() que descarga un mapa para un lugar dado. También podría ser una dirección, latitud/longitud o bounding box. El argumento del tipo de mapa permite indicar el estilo o tipo de mapa. Más información sobre esto se puede encontrar en la ayuda ?get_map

Cuando se hace un gráfico con ggplot se suele empezar con la función ggplot(). En este caso se empieza con ggmap() que incluye el objeto con el mapa de fondo. Entonces se añaden los puntos de las peluquerías de Murcia con geom_sf(). Es importante indicar con el argumento inherit.aes = FALSE que use los aesthetic mappings del objeto espacial osm_points. Además, se cambia el color, relleno, transparencia tipo y tamaño de los círculos.

    murcia_se_peina <- get_map(c(left = -1.18, bottom = 37.91, right = -1.05, top = 38.04), maptype = "watercolor")

    ggmap(murcia_se_peina) +
    geom_sf(data = pelus$osm_points, inherit.aes = FALSE, colour = "#238443", fill = "#004529", alpha = .5, size = 4, shape = 21) +
    labs(x = "", y = "")

/flow/hugo-ox-mpvd/src/branch/main/murcia_se_peina.pdf

Dónde hay una zapatería

q_murcia_shoes <- getbb("Murcia Spain") %>%
  opq() %>%
   add_osm_feature("shop", "shoes")
str(q_murcia_shoes)
List of 4
 $ bbox    : chr "37.7158719,-1.3848335,38.1179185,-0.8508899"
 $ prefix  : chr "[out:xml][timeout:25];\n(\n"
 $ suffix  : chr ");\n(._;>;);\nout body;"
 $ features: chr " [\"shop\"=\"shoes\"]"
 - attr(*, "class")= chr [1:2] "list" "overpass_query"
 - attr(*, "nodes_only")= logi FALSE
murcia_shoes <- osmdata_sf(q_murcia_shoes)
murcia_shoes
Object of class 'osmdata' with:
                 $bbox : 37.7158719,-1.3848335,38.1179185,-0.8508899
        $overpass_call : The call submitted to the overpass API
                 $meta : metadata including timestamp and version numbers
           $osm_points : 'sf' Simple Features Collection with 16 points
            $osm_lines : NULL
         $osm_polygons : 'sf' Simple Features Collection with 0 polygons
       $osm_multilines : NULL
    $osm_multipolygons : NULL
  murcia_shoes <- get_map(getbb("Murcia Spain"), maptype = "watercolor", source = "stamen")

  ggmap(murcia_shoes) +
  geom_sf(data = q_murcia_shoes$osm_points, inherit.aes = FALSE, colour = "#238443", fill = "#004529", alpha = .5, size = 4, shape = 21) +
  labs(x = "", y = "")

/flow/hugo-ox-mpvd/src/branch/main/murcia_shoes.pdf

Mercadona

En vez de obtener un bounding box con la función getbb() se puede crear sabiendo los puntos oeste, sur, este y norte. En la consulta hay dos características: el nombre y la tienda para filtrar supermercados de esta marca en particular. Según el area o volumen de la consulta se puede alargar el tiempo de espera. Por defecto el límite se establece en 25 segundos antes del timeout.

Lo que nos interesa son los puntos donde hay supermercados por lo que se usa la geometría a través de geom_sf(). La función theme_void() borra todo menos los puntos.

Se crean los límites:

m <- c(-10, 30, 5, 46)

Se hace la consulta:

    q <- m %>%
    opq (timeout = 25*100) %>%
    add_osm_feature("name", "Mercadona") %>%
    add_osm_feature("shop", "supermarket")
  mercadona <- osmdata_sf(q)
mercadona
Object of class 'osmdata' with:
                 $bbox : 30,-10,46,5
        $overpass_call : The call submitted to the overpass API
                 $meta : metadata including timestamp and version numbers
           $osm_points : 'sf' Simple Features Collection with 6478 points
            $osm_lines : 'sf' Simple Features Collection with 25 linestrings
         $osm_polygons : 'sf' Simple Features Collection with 631 polygons
       $osm_multilines : NULL
    $osm_multipolygons : 'sf' Simple Features Collection with 11 multipolygons

Y el mapa final:

  ggplot(mercadona$osm_points)+
  geom_sf(colour = "#08519c",
	fill = "#08306b",
      alpha = 5,
    size = 1,
  shape = 21) +
  theme_void()

/flow/hugo-ox-mpvd/src/branch/main/mercadonas.pdf

Visualizar el crecimiento urbano

La Dirección General del Catastro de España dispone de información espacial de toda la edificación a excepción del País Vasco y Navarra. Este conjunto de datos forma parte de la implantación de INSPIRE, la Infraestructura de Información Espacial en Europa. Más información podemos encontrar aquí. Utilizaremos los enlaces (urls) en formato ATOM, que es un formato de redifusión de tipo RSS, permitiéndonos obtener el enlace de descarga para cada municipio.

Esta entrada de blog es una versión reducida del caso práctico que podéis encontrar en nuestra reciente publicación - Introducción a los SIG con R - publicado por Dominic Royé y Roberto Serrano-Notivoli.

Paquetes Paquete Descripción tidyverse Conjunto de paquetes (visualización y manipulación de datos): ggplot2, dplyr, purrr,etc. sf Simple Feature: importar, exportar y manipular datos vectoriales fs Proporciona una interfaz uniforme y multiplataforma para las operaciones del sistema de archivos lubridate Fácil manipulación de fechas y tiempos feedeR Importar formatos de redifusión RSS tmap Fácil creación de mapas temáticos classInt Para crear intervalos de clase univariantes sysfonts Carga familias tipográficas del sistema y de Google showtext Usar familias tipográficas más fácilmente en gráficos R

if(!require("tidyverse")) install.packages("tidyverse") if(!require("feedeR")) install.packages("feedeR") if(!require("fs")) install.packages("fs") if(!require("lubridate")) install.packages("lubridate") if(!require("fs")) install.packages("fs") if(!require("tmap")) install.packages("tmap") if(!require("classInt")) install.packages("classInt") if(!require("showtext")) install.packages("showtext") if(!require("sysfonts")) install.packages("sysfonts") if(!require("rvest")) install.packages("rvest")

library(feedeR) library(sf) library(fs) library(tidyverse) library(lubridate) library(classInt) library(tmap) library(rvest) Enlaces de descarga La primera url nos dará acceso a un listado de provincias, sedes territoriales (no siempre coinciden con la provincia), con nuevos enlaces RSS los cuales incluyen los enlaces finales de descarga para cada municipio. En este caso, descargaremos el edificado de Valencia. Los datos del Catastro se actualizan cada seis meses.

url <- "http://www.catastro.minhap.es/INSPIRE/buildings/ES.SDGC.bu.atom.xml"

prov_enlaces <- feed.extract(url) str(prov_enlaces) #estructura es lista ## List of 4 ## $ title : chr "Download service of Buildings. Territorial Office" ## $ link : chr "http://www.catastro.minhap.es/INSPIRE/buildings/ES.SDGC.BU.atom.xml" ## $ updated: POSIXct[1:1], format: "2019-10-26" ## $ items :'data.frame': 52 obs. of 4 variables: ## ..$ title: chr [1:52] "Territorial office 02 Albacete" "Territorial office 03 Alicante" "Territorial office 04 Almería" "Territorial office 05 Avila" ... ## ..$ date : POSIXct[1:52], format: "2019-10-26" "2019-10-26" … ## ..$ link : chr [1:52] "http://www.catastro.minhap.es/INSPIRE/buildings/02/ES.SDGC.bu.atom_02.xml" "http://www.catastro.minhap.es/INSPIRE/buildings/03/ES.SDGC.bu.atom_03.xml" "http://www.catastro.minhap.es/INSPIRE/buildings/04/ES.SDGC.bu.atom_04.xml" "http://www.catastro.minhap.es/INSPIRE/buildings/05/ES.SDGC.bu.atom_05.xml" ... ## ..$ hash : chr [1:52] "d21ebb7975e59937" "bdba5e149f09e9d8" "03bcbcc7c5be2e17" "8a154202dd778143" …

prov_enlaces_tab <- as_tibble(prov_enlaces$items) %>% mutate(title = repair_encoding(title)) ## Best guess: UTF-8 (100% confident) prov_enlaces_tab ## # A tibble: 52 x 4 ## title date link hash ## <chr> <dttm> <chr> <chr> ## 1 Territorial of~ 2019-10-26 00:00:00 http://www.catastro.minhap.es/~ d21ebb79~ ## 2 Territorial of~ 2019-10-26 00:00:00 http://www.catastro.minhap.es/~ bdba5e14~ ## 3 Territorial of~ 2019-10-26 00:00:00 http://www.catastro.minhap.es/~ 03bcbcc7~ ## 4 Territorial of~ 2019-10-26 00:00:00 http://www.catastro.minhap.es/~ 8a154202~ ## 5 Territorial of~ 2019-10-26 00:00:00 http://www.catastro.minhap.es/~ 7d3fd376~ ## 6 Territorial of~ 2019-10-26 00:00:00 http://www.catastro.minhap.es/~ 9c08741f~ ## 7 Territorial of~ 2019-10-26 00:00:00 http://www.catastro.minhap.es/~ ff722b15~ ## 8 Territorial of~ 2019-10-26 00:00:00 http://www.catastro.minhap.es/~ b431aa61~ ## 9 Territorial of~ 2019-10-26 00:00:00 http://www.catastro.minhap.es/~ f79c6562~ ## 10 Territorial of~ 2019-10-26 00:00:00 http://www.catastro.minhap.es/~ d702a6a8~ ## # … with 42 more rows Accedemos y descargamos los datos de Valencia. Para encontrar el enlace final de descarga usamos la función filter() del paquete dplyr buscando el nombre de la sede territorial y posteriormente el nombre del municipio en mayúsculas con la función str_detect() de stringr. La función pull() nos permite extraer una columna de un data.frame.

Actualmente la función feed.extract() no importa correctamente en el encoding UTF-8 (Windows). Por eso, en algunas ciudades pueden aparecer una mala codificación de caracteres especiales “Cádiz”. Para subsanar este problema aplicamos la función repair_encoding() del paquete rvest.

val_atom <- filter(prov_enlaces_tab, str_detect(title, "Valencia")) %>% pull(link)

val_enlaces <- feed.extract(val_atom)

val_enlaces_tab <- val_enlaces$items val_enlaces_tab <- mutate(val_enlaces_tab, title = repair_encoding(title), link = repair_encoding(link)) ## Best guess: UTF-8 (100% confident) ## Best guess: UTF-8 (100% confident)

val_link <- filter(val_enlaces_tab, str_detect(title, "VALENCIA")) %>% pull(link) val_link ## [1] "http://www.catastro.minhap.es/INSPIRE/Buildings/46/46900-VALENCIA/A.ES.SDGC.BU.46900.zip" Descarga de datos La descarga se realiza con la función download.file() que únicamente tiene dos argumentos principales, el enlace de descarga y la ruta con el nombre del archivo. En este caso hacemos uso de la función tempfile(), que nos es útil para crear archivos temporales, es decir, archivos que únicamente existen en la memoría RAM por un tiempo determinado. El archivo que descargamos tiene extensión *.zip, por lo que debemos descomprimirlo con otra función (unzip()), que requiere el nombre del archivo y el nombre de la carpeta donde lo queremos descomprimir. Por último, la función URLencode() codifica una dirección URL que contiene caracteres especiales.

temp <- tempfile()

download.file(URLencode(val_link), temp)

unzip(temp, exdir = "buildings") Importar los datos Para importar los datos utilizamos la función dir_ls() del paquete fs, que nos permite obtener los archivos y carpetas de una ruta concreta al mismo tiempo que filtramos por un patrón de texto (regexp: expresión regular). Aplicamos la función st_read() del paquete sf al archivo espacial de formato Geography Markup Language (GML).

file_val <- dir_ls("buildings", regexp = "building.gml")

buildings_val <- st_read(file_val) ## Reading layer `Building' from data source `C:\Users\xeo19\Documents\GitHub\blogR_update\content\post\es\2019-11-01-visualizar-crecimiento-urbano\buildings\A.ES.SDGC.BU.46900.building.gml' using driver `GML' ## Simple feature collection with 36296 features and 24 fields ## geometry type: MULTIPOLYGON ## dimension: XY ## bbox: xmin: 720608 ymin: 4351287 xmax: 734982.5 ymax: 4382906 ## CRS: 25830 Preparación de los datos Únicamente convertimos la columna de la edad del edificio (beginning) en clase Date. La columna de la fecha contiene algunas fechas en formato 01-01 lo que no corresponde a ninguna fecha reconocible. Por eso, reemplazamos el primer - por 0000.

buildings_val <- mutate(buildings_val, beginning = str_replace(beginning, "^-", "0000") %>% ymd_hms() %>% as_date() ) ## Warning: 4 failed to parse. Gráfico de distribución Antes de crear el mapa de la edad del edificado, lo que reflejará el crecimiento urbano, haremos un gráfico de distribución de la fecha de construcción de los edificios. Podremos identificar claramente períodos de expansión urbana. Usaremos el paquete ggplot2 con la geometría de geom_density() para este objetivo. La función font_add_google() del paquete sysfonts nos permite descargar e incluir familias tipográficas de Google.

#descarga de familia tipográfica sysfonts::font_add_google("Montserrat", "Montserrat")

#usar showtext para familias tipográficas showtext::showtext_auto() #limitamos al periodo posterior a 1750 filter(buildings_val, beginning >= "1750-01-01") %>% ggplot(aes(beginning)) + geom_density(fill = "#2166ac", alpha = 0.7) + scale_x_date(date_breaks = "20 year", date_labels = "%Y") + theme_minimal() + theme(title = element_text(family = "Montserrat"), axis.text = element_text(family = "Montserrat")) + labs(y = "",x = "", title = "Evolución del desarrollo urbano")

Buffer de 2,5 km de Valencia Para poder visualizar bien la distribución del crecimiento, limitamos el mapa a un radio de 2,5 km desde el centro de la ciudad. Usamos la función geocode_OSM() del paquete tmaptools para obtener las coordenadas de Valencia en clase sf. Después proyectamos los puntos al sistema que usamos para el edificado (EPSG:25830). Como último paso creamos con la función st_buffer() un buffer con 2500 m y la intersección con nuestros datos de los edificios. También es posible crear un buffer en forma de un rectángulo indicando el tipo de estilo con el argumento endCapStyle = "SQUARE".

ciudad_point <- tmaptools::geocode_OSM("Valencia", as.sf = TRUE)

ciudad_point <- st_transform(ciudad_point, 25830)

point_bf <- st_buffer(ciudad_point, 2500)

buildings_val25 <- st_intersection(buildings_val, point_bf) ## Warning: attribute variables are assumed to be spatially constant throughout all ## geometries Preparar los datos para el mapas Para poder visualizar bien las diferentes épocas de crecimiento, categorizamos el año en 15 grupos empleando cuartiles.

#encontrar 15 clases br <- classIntervals(year(buildings_val25$beginning), 15, "quantile") ## Warning in classIntervals(year(buildings_val25$beginning), 15, "quantile"): var ## has missing values, omitted in finding classes #crear etiquetas lab <- names(print(br, under = "<", over = ">", cutlabels = FALSE)) ## style: quantile ## < 1890 1890 - 1912 1912 - 1925 1925 - 1930 1930 - 1940 1940 - 1950 ## 940 1369 971 596 1719 1080 ## 1950 - 1957 1957 - 1962 1962 - 1966 1966 - 1970 1970 - 1973 1973 - 1977 ## 1227 1266 1233 1165 1161 932 ## 1977 - 1987 1987 - 1999 > 1999 ## 1337 1197 1190 #categorizar el año buildings_val25 <- mutate(buildings_val25, yr_cl = cut(year(beginning), br$brks, labels = lab, include.lowest = TRUE)) Mapa de Valencia El mapa creamos con el paquete tmap. Es una interesante alternativa a ggplot2. Se trata de un paquete de funciones especializadas en crear mapas temáticos. La filosofía del paquete sigue a la de ggplot2, creando multiples capas con diferentes funciones, que siempre empiezan con tm_* y se combinan con +. La construcción de un mapa con tmap siempre comienza con tm_shape(), donde se definen los datos que queremos dibujar. Luego agregamos la geometría correspondiente al tipo de datos (tm_polygon(), tm_border(), tm_dots() o incluso tm_raster()). La función tm_layout() ayuda a configurar el estilo del mapa.

Cuando necesitamos más colores del máximo permitido por RColorBrewer podemos pasar los colores a la función colorRampPalette(). Esta función interpola para un mayor número más colores de la gama.

#colores col_spec <- RColorBrewer::brewer.pal(11, "Spectral")

#función de una gama de colores col_spec_fun <- colorRampPalette(col_spec)

#crear los mapas tm_shape(buildings_val25) + tm_polygons("yr_cl", border.col = "transparent", palette = col_spec_fun(15), textNA = "Sin dato", title = "") + tm_layout(bg.color = "black", outer.bg.color = "black", legend.outside = TRUE, legend.text.color = "white", legend.text.fontfamily = "Montserrat", panel.label.fontfamily = "Montserrat", panel.label.color = "white", panel.label.bg.color = "black", panel.label.size = 5, panel.label.fontface = "bold")

Mapa dinámico en leaflet Una ventaja muy interesante es la función tmap_leaflet() del paquete tmap para pasar de forma sencilla un mapa creado en el mismo marco a leaflet.

#mapa tmap de Santiago m <- tm_shape(buildings_val25) + tm_polygons("yr_cl", border.col = "transparent", palette = col_spec_fun(15), textNA = "Without data", title = "")

#mapa dinámico tmap_leaflet(m)

Visualizar las anomalías climáticas

When we visualize precipitation and temperature anomalies, we simply use time series as bar graph indicating negative and positive values in red and blue. However, in order to have a better overview we need both anomalies in a single graph. In this way we could more easly answer the question of whether a particular season or month was dry-warm or wet-cold, and even compare these anomalies in the context of previous years.

Packages In this post we will use the following packages:

Package Description tidyverse Collection of packages (visualization, manipulation): ggplot2, dplyr, purrr, etc. lubridate Easy manipulation of dates and times ggrepel Repel overlapping text labels in ggplot2 #we install the packages if necessary if(!require("tidyverse")) install.packages("tidyverse") if(!require("ggrepel")) install.packages("ggrepel") if(!require("lubridate")) install.packages("lubridate")

#packages library(tidyverse) library(lubridate) library(ggrepel) Preparing the data First we import the daily precipitation and temperature data from the selected weather station (download). We will use the data from Tenerife South (Spain) [1981-2020] accessible through Open Data AEMET. In R there is a package called meteoland that facilitates the download with specific functions to access data from AEMET (Spanish State Meteorological Agency), Meteogalicia (Galician Meteorological Service) and Meteocat (Catalan Meteorological Service).

Step 1: import the data We import the data in csv format, the first column is the date, the second column the precipitation (pr) and the last column the average daily temperature (ta).

data <- read_csv("meteo_tenerife.csv") ## Parsed with column specification: ## cols( ## date = col_date(format = ""), ## pr = col_double(), ## ta = col_double() ## ) data ## # A tibble: 14,303 x 3 ## date pr ta ## <date> <dbl> <dbl> ## 1 1981-01-02 0 17.6 ## 2 1981-01-03 0 16.8 ## 3 1981-01-04 0 17.4 ## 4 1981-01-05 0 17.6 ## 5 1981-01-06 0 17 ## 6 1981-01-07 0 17.6 ## 7 1981-01-08 0 18.6 ## 8 1981-01-09 0 19.8 ## 9 1981-01-10 0 21.5 ## 10 1981-01-11 3.8 17.6 ## # … with 14,293 more rows Step 2: preparing the data In the second step we prepare the data to calculate the anomalies. To do this, we create three new columns: the month, the year, and the season of the year. Since our objective is to analyse winter anomalies, we cannot use the calendar year, because winter includes the month of December of one year and the months of January and February of the following. The custom function meteo_yr() extracts the year from a date indicating the starting month. The concept is similar to the hydrological year in which it starts on October 1.

meteo_yr <- function(dates, start_month = NULL) {

dates.posix <- as.POSIXlt(dates)

offset <- ifelse(dates.posix$mon >= start_month - 1, 1, 0)

adj.year = dates.posix$year + 1900 + offset return(adj.year) } We will use many functions of the package collection tidyverse (https://www.tidyverse.org/). The mutate() function helps to add new columns or change existing ones. To define the seasons, we use the case_when() function from the dplyr package, which has many advantages compared to a chain of ifelse(). In case_when() we use two-side formulas, on the one hand the condition and on the other the action when that condition is met. A two-sided formula in R consists of the operator ~. The binary operator %in% allows us to filter several values in a greater set.

data <- mutate(data, winter_yr = meteo_yr(date, 12), month = month(date), season = case_when(month %in% c(12,1:2) ~ "Winter", month %in% 3:5 ~ "Spring", month %in% 6:8 ~ "Summer", month %in% 9:11 ~ "Autum"))

data ## # A tibble: 14,303 x 6 ## date pr ta winter_yr month season ## <date> <dbl> <dbl> <dbl> <dbl> <chr> ## 1 1981-01-02 0 17.6 1981 1 Winter ## 2 1981-01-03 0 16.8 1981 1 Winter ## 3 1981-01-04 0 17.4 1981 1 Winter ## 4 1981-01-05 0 17.6 1981 1 Winter ## 5 1981-01-06 0 17 1981 1 Winter ## 6 1981-01-07 0 17.6 1981 1 Winter ## 7 1981-01-08 0 18.6 1981 1 Winter ## 8 1981-01-09 0 19.8 1981 1 Winter ## 9 1981-01-10 0 21.5 1981 1 Winter ## 10 1981-01-11 3.8 17.6 1981 1 Winter ## # … with 14,293 more rows Step 3: estimate winter anomalies In the next step we create a subset of the winter months. Then we group by the defined meteorological year and calculate the sum and average for precipitation and temperature, respectively. To facilitate the work, the magrittr package introduces the operator called pipe in the form %>% with the aim of combining several functions without the need to assign the result to a new object. The pipe operator passes the output of a function applied to the first argument of the next function. This way of combining functions allows you to chain several steps simultaneously. The %>% must be understood and pronounced as then.

data_inv <- filter(data, season == "Winter") %>% group_by(winter_yr) %>% summarise(pr = sum(pr, na.rm = TRUE), ta = mean(ta, na.rm = TRUE)) Now we only have to calculate the anomalies of precipitation and temperature. The columns pr_mean and ta_mean will contain the climate average, the reference for the anomalies with respect to the normal period 1981-2010. Therefore, we need to filter the values to the period before 2010, which we will do in the usual way of filtering vectors in R. Once we have the references we estimate the anomalies pr_anom and ta_anom. To facilitate the interpretation, in the case of precipitation we express the anomalies as percentage, with the average set at 0% instead of 100%.

In addition, we add three required columns with information for the creation of the graph: 1) labyr contains the year of each anomaly as long as it has been greater/less than -+10% or -+0.5ºC, respectively (this is for reducing the number of labels), 2) symb_point is a dummy variable in order to be able to create different symbols between the cases of (1), and 3) lab_font for highlighting in bold the year 2020.

data_inv <- mutate(data_inv, pr_mean = mean(pr[winter_yr <= 2010]), ta_mean = mean(ta[winter_yr <= 2010]), pr_anom = (pr*100/pr_mean)-100, ta_anom = ta-ta_mean,

labyr = case_when(pr_anom < -10 & ta_anom < -.5 ~ winter_yr, pr_anom < -10 & ta_anom > .5 ~ winter_yr, pr_anom > 10 & ta_anom < -.5 ~ winter_yr, pr_anom > 10 & ta_anom > .5 ~ winter_yr), symb_point = ifelse(!is.na(labyr), "yes", "no"), lab_font = ifelse(labyr == 2020, "bold", "plain") ) Creating the graph We will build the chart adding layer by layer the distinctive elements: 1) the background with the different grids (Dry-Warm, Dry-Cold, etc.), 2) the points and labels, and 3) the style adjustments.

Part 1 The idea is that the points with dry-warm anomalies are located in quadrant I (top-right) and those with wet-cold in quadrant III (bottom-left). Therefore, we must invert the sign in the precipitation anomalies. Then we create a data.frame with the label positions of the four quadrants. For the positions in x and y Inf and -Inf are used, which is equivalent to the maximum panel sides with respect to the data. However, it is necessary to adjust the position towards the extreme points within the panel with the known arguments of ggplot2: hjust and vjust.

data_inv_p <- mutate(data_inv, pr_anom = pr_anom * -1)

bglab <- data.frame(x = c(-Inf, Inf, -Inf, Inf), y = c(Inf, Inf, -Inf, -Inf), hjust = c(1, 1, 0, 0), vjust = c(1, 0, 1, 0), lab = c("Wet-Warm", "Dry-Warm", "Wet-Cold", "Dry-Cold"))

bglab ## x y hjust vjust lab ## 1 -Inf Inf 1 1 Wet-Warm ## 2 Inf Inf 1 0 Dry-Warm ## 3 -Inf -Inf 0 1 Wet-Cold ## 4 Inf -Inf 0 0 Dry-Cold Part 2 In the second part we can start building the chart by adding all graphical elements. First we create the background with different colors of each quadrant. The function annotate() allows adding geometry layers without the use of variables within data.frames. With the geom_hline() and geom_vline() function we mark the quadrants horizontally and vertically using a dashed line. Finally, we draw the labels of each quadrant, using the function geom_text(). When we use other data sources than the main one used in ggplot(), we must indicate it with the argument data in the corresponding geometry function.

g1 <- ggplot(data_inv_p, aes(pr_anom, ta_anom)) + annotate("rect", xmin = -Inf, xmax = 0, ymin = 0, ymax = Inf, fill = "#fc9272", alpha = .6) + #wet-warm annotate("rect", xmin = 0, xmax = Inf, ymin = 0, ymax = Inf, fill = "#cb181d", alpha = .6) + #dry-warm annotate("rect", xmin = -Inf, xmax = 0, ymin = -Inf, ymax = 0, fill = "#2171b5", alpha = .6) + #wet-cold annotate("rect", xmin = 0, xmax = Inf, ymin = -Inf, ymax = 0, fill = "#c6dbef", alpha = .6) + #dry-cold geom_hline(yintercept = 0, linetype = "dashed") + geom_vline(xintercept = 0, linetype = "dashed") + geom_text(data = bglab, aes(x, y, label = lab, hjust = hjust, vjust = vjust), fontface = "italic", size = 5, angle = 90, colour = "white")

g1

Part 3 In the third part we simply add the points of the anomalies and the labels of the years. The geom_text_repel() function is similar to the one known by default in ggplot2, geom_text(), but it repels overlapping text labels away from each other.

g2 <- g1 + geom_point(aes(fill = symb_point, colour = symb_point), size = 2.8, shape = 21, show.legend = FALSE) + geom_text_repel(aes(label = labyr, fontface = lab_font), max.iter = 5000, size = 3.5) g2 ## Warning: Removed 25 rows containing missing values (geom_text_repel).

Part 4 In the last part we adjust, in addition to the general style, the axes, the color type and the (sub)title. Remember that we changed the sign on precipitation anomalies. Hence, we must use the arguments breaks and labels in the function scale_x_continouous() to reverse the sign in the labels corresponding to the breaks.

g3 <- g2 + scale_x_continuous("Precipitation anomaly in %", breaks = seq(-100, 250, 10) * -1, labels = seq(-100, 250, 10), limits = c(min(data_inv_p$pr_anom), 100)) + scale_y_continuous("Mean temperature anomaly in ºC", breaks = seq(-2, 2, 0.5)) + scale_fill_manual(values = c("black", "white")) + scale_colour_manual(values = rev(c("black", "white"))) + labs(title = "Winter anomalies in Tenerife South", caption = "Data: AEMET\nNormal period 1981-2010") + theme_bw()

g3 ## Warning: Removed 25 rows containing missing values (geom_text_repel).

Visualizar el crecimiento urbano

Cuenta Dominic Royé que la Dirección General del Catastro de España dispone de información espacial de toda la edificación a excepción del País Vasco y Navarra. Este conjunto de datos forma parte de la implantación de INSPIRE, la Infraestructura de Información Espacial en Europa. Más información podemos encontrar aquí. Utilizaremos los enlaces (urls) en formato ATOM, que es un formato de redifusión de tipo RSS, permitiéndonos obtener el enlace de descarga para cada municipio.

Instalamos y cargamos los paquetes

tidyverse
metapaquete para visualización y manipulación de datos con ggplot2, dplyr o purrr.
sf
de Simple Feature, sirve para importar, exportar y manipular datos vectoriales.
fs
proporciona una interfaz uniforme y multiplataforma para las operaciones del sistema de archivos.
lubridate
para manipulación de fechas y tiempos.
feedeR
para importar formatos de redifusión RSS.
tmap
para creación de mapas temáticos.
classInt
para crear intervalos de clase univariantes.
sysfonts
carga familias tipográficas del sistema y de Google.
showtext
para usar familias tipográficas en gráficos de R.
if(!require("tidyverse")) install.packages("tidyverse")
if(!require("feedeR")) install.packages("feedeR")
if(!require("fs")) install.packages("fs")
if(!require("lubridate")) install.packages("lubridate")
if(!require("fs")) install.packages("fs")
if(!require("tmap")) install.packages("tmap")
if(!require("classInt")) install.packages("classInt")
if(!require("showtext")) install.packages("showtext")
if(!require("sysfonts")) install.packages("sysfonts")
if(!require("rvest")) install.packages("rvest")

library(feedeR)
library(sf) 
library(fs)
library(tidyverse)
library(lubridate)
library(classInt)
library(tmap)
library(rvest)

COMMENT Local Variables   ARCHIVE