Mapas
Teresa Ortiz
Paquetes a utilizar:
library(ggplot2)
library(plyr)
library(dplyr)
library(scales)
library(maptools)
library(rgdal)
library(ggmap)
library(gridExtra)
gpclibPermit()## [1] TRUE¿Por qué mapas?
Gran parte de los datos que se estudian en estadística espacial consisten en observaciones en la Tierra por lo que es importante realizar mapas donde grafiquemos las observaciones, o realizar mapas que incorporen estimaciones de un modelo. Los mapas pueden ser parte de la presentación de resultados, del análisis exploratorio o pueden servir para evaluar el ajuste de un modelo.
La creación de mapas requiere conocimientos técnicos (uso de software) y entendimiento de los principios de sistemas de referencia de coordenadas (CRS). Es por ello que comenzaremos definiendo algunos conceptos de los CRS.
Elipses, datums, proyecciones
- Elipse: Describe la forma de la Tierra, es una aproximación que no ajusta a la perfección. Existen distintos elipsoides, algunos están diseñados para ajustar toda la Tierra (WGS84, GRS80) y algunos ajustan únicamente por regiones (NAD27). Los locales son más precisos para el área para la que fueron diseñados pero no funcionan en otras partes del mundo.
projInfo(type = "ellps")[1:10, ]## name major ell description
## 1 MERIT a=6378137.0 rf=298.257 MERIT 1983
## 2 SGS85 a=6378136.0 rf=298.257 Soviet Geodetic System 85
## 3 GRS80 a=6378137.0 rf=298.257222101 GRS 1980(IUGG, 1980)
## 4 IAU76 a=6378140.0 rf=298.257 IAU 1976
## 5 airy a=6377563.396 b=6356256.910 Airy 1830
## 6 APL4.9 a=6378137.0. rf=298.25 Appl. Physics. 1965
## 7 NWL9D a=6378145.0. rf=298.25 Naval Weapons Lab., 1965
## 8 mod_airy a=6377340.189 b=6356034.446 Modified Airy
## 9 andrae a=6377104.43 rf=300.0 Andrae 1876 (Den., Iclnd.)
## 10 aust_SA a=6378160.0 rf=298.25 Australian Natl & S. Amer. 1969- Datum: Define un punto de origen para el sistema de coordenadas de la Tierra y la dirección de los ejes. El datum define el elipsoide que se usará pero el sentido contrario no es cierto.
projInfo(type = "datum")## name ellipse
## 1 WGS84 WGS84
## 2 GGRS87 GRS80
## 3 NAD83 GRS80
## 4 NAD27 clrk66
## 5 potsdam bessel
## 6 carthage clrk80ign
## 7 hermannskogel bessel
## 8 ire65 mod_airy
## 9 nzgd49 intl
## 10 OSGB36 airy
## definition
## 1 towgs84=0,0,0
## 2 towgs84=-199.87,74.79,246.62
## 3 towgs84=0,0,0
## 4 nadgrids=@conus,@alaska,@ntv2_0.gsb,@ntv1_can.dat
## 5 towgs84=598.1,73.7,418.2,0.202,0.045,-2.455,6.7
## 6 towgs84=-263.0,6.0,431.0
## 7 towgs84=577.326,90.129,463.919,5.137,1.474,5.297,2.4232
## 8 towgs84=482.530,-130.596,564.557,-1.042,-0.214,-0.631,8.15
## 9 towgs84=59.47,-5.04,187.44,0.47,-0.1,1.024,-4.5993
## 10 towgs84=446.448,-125.157,542.060,0.1502,0.2470,0.8421,-20.4894
## description
## 1
## 2 Greek_Geodetic_Reference_System_1987
## 3 North_American_Datum_1983
## 4 North_American_Datum_1927
## 5 Potsdam Rauenberg 1950 DHDN
## 6 Carthage 1934 Tunisia
## 7 Hermannskogel
## 8 Ireland 1965
## 9 New Zealand Geodetic Datum 1949
## 10 Airy 1830- Proyección: Proyectar el elipsoide en un espacio de dos dimensiones. Es decir pasar de longitud y latitud en el elipsoide a coordenadas de un mapa, esto conlleva necesariamente una distorsión de la superficie. La estrategia usual es utilizar superficies de desarrollo y después aplanar la superficie. Todas las proyecciones inducen alguna distorsión y las proyecciones de los mapas serán diferentes.
Las distorsiones resultan en pérdidas de información que pueden ser en área, forma, distancia o dirección. Por ejemplo Mercator preserva dirección y es útil para navegación:
states <- map_data("state")
usmap <- ggplot(states, aes(x=long, y=lat, group=group)) +
geom_polygon(fill="white", colour="black")
usmap + coord_map("mercator")
Por su parte Azimuthal Equal Area preserva area pero no dirección:
usmap + coord_map("azequalarea")
En particular en investigación es común usar Universal Transverse Mercator (UTM) porque preserva ángulos y dirección pero distorsiona distancia. Para minimizar la distorsión se divide la Tierra en 60 regiones y utiliza una proyección (de secantes) en cada una.
En R podemos ver la lista de proyecciones con la siguiente instrucción:
# 135 proyecciones distintas
projInfo(type="proj")[1:20, ]## name description
## 1 aea Albers Equal Area
## 2 aeqd Azimuthal Equidistant
## 3 airy Airy
## 4 aitoff Aitoff
## 5 alsk Mod. Stererographics of Alaska
## 6 apian Apian Globular I
## 7 august August Epicycloidal
## 8 bacon Bacon Globular
## 9 bipc Bipolar conic of western hemisphere
## 10 boggs Boggs Eumorphic
## 11 bonne Bonne (Werner lat_1=90)
## 12 calcofi Cal Coop Ocean Fish Invest Lines/Stations
## 13 cass Cassini
## 14 cc Central Cylindrical
## 15 cea Equal Area Cylindrical
## 16 chamb Chamberlin Trimetric
## 17 collg Collignon
## 18 crast Craster Parabolic (Putnins P4)
## 19 denoy Denoyer Semi-Elliptical
## 20 eck1 Eckert IEs posible trabajar con datos no proyectados (longitud/latitud) pero se requieren métodos para medir distancias: gran círculo, naive euclideana, o cordal (esto último pensando en modelado).
Sistemas de referencia de coordenadas (CRS)
Los CRS nos dan una manera estándar de describir ubicaciones en la Tierra. Si combinamos información de distintos CRS debemos transformarlos para poder alinear la información.
En R, la notación que se utiliza para describir un CRS es proj4string de la libraría PROJ.4 y se ven como:
+init=epsg:4121 +proj=longlat +ellps=GRS80 +datum=GGRS87 +no_defs +towgs84=-199.87,74.79,246.62
Los shapefiles suelen tener asociada una proyección.
ogrInfo(dsn = "shapes", layer = "Municipios")## Source: "shapes", layer: "Municipios"
## Driver: ESRI Shapefile; number of rows: 2456
## Feature type: wkbPolygon with 2 dimensions
## Extent: (907821.8 319149.1) - (4082997 2349615)
## CRS: +proj=lcc +lat_1=17.5 +lat_2=29.5 +lat_0=12 +lon_0=-102 +x_0=2500000 +y_0=0 +ellps=GRS80 +units=m +no_defs
## LDID: 87
## Number of fields: 6
## name type length typeName
## 1 OBJECTID 0 9 Integer
## 2 CVE_ENT 4 2 String
## 3 CVE_MUN 4 3 String
## 4 NOM_MUN 4 70 String
## 5 Shape_Leng 2 19 Real
## 6 Shape_Area 2 19 RealogrInfo(dsn = "shapes", layer = "Mex_Edos")## Source: "shapes", layer: "Mex_Edos"
## Driver: ESRI Shapefile; number of rows: 32
## Feature type: wkbPolygon with 2 dimensions
## Extent: (-118.3919 14.58918) - (-86.72716 32.70473)
## CRS: +proj=longlat +datum=WGS84 +no_defs
## LDID: 89
## Number of fields: 1
## name type length typeName
## 1 NOM_ENT 4 70 StringEn el caso de la proyección en los shapes de Municipios se usó North America Lambert Conformal Conic Projection.
Hay varias maneras de transformar el CRS de un objeto espacial, una de ellas es usar transform.
mun_shp <- readOGR("shapes" , "Municipios")## OGR data source with driver: ESRI Shapefile
## Source: "shapes", layer: "Municipios"
## with 2456 features
## It has 6 fieldsmun_shp_2 <- spTransform(mun_shp, CRS("+proj=longlat +datum=WGS84"))
mun_shp@proj4string## CRS arguments:
## +proj=lcc +lat_1=17.5 +lat_2=29.5 +lat_0=12 +lon_0=-102
## +x_0=2500000 +y_0=0 +ellps=GRS80 +units=m +no_defsmun_shp_2@proj4string## CRS arguments:
## +proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0Shapefile
Un shapefile es un grupo de archivos que contienen geometrías e información de acuerdo a un estándar especificado por el Insituto de Investigación en Sistemas de Ecosistemas (ESRI). Nosotros tenemos los siguientes grupos de archivos (para estados):
Estados.shp: es el archivo principal y contiene la geometría correspondiente.
Estados.dbf: es la base de datos y almacena la información de los atributos de los objetos
- Estados.shx: almacena el índice de las entidades geométricas.
- Estados.prj (opcional): información del CRS.
Estados.sbn y .sbx (opcional): almacena índice espacial.
Veamos las librerías que usaremos
library(ggplot2)
library(plyr)
library(dplyr)
library(scales)
library(maptools)
library(rgdal)
library(ggmap)
library(gridExtra)Usamos la función readOGR para leer los archivos de estados.
edo_shp <- readOGR("shapes" , "Mex_Edos")## OGR data source with driver: ESRI Shapefile
## Source: "shapes", layer: "Mex_Edos"
## with 32 features
## It has 1 fieldsel “shapes” indica que los archivos están en el directorio shapes.
Notemos que el objeto edo_shp no es un data frame,
class(edo_shp)## [1] "SpatialPolygonsDataFrame"
## attr(,"package")
## [1] "sp"Lo podemos graficar directamente con plot.
plot(edo_shp)
pero para poder graficarlo con ggplot debemos convertirlo en data frame.
Formamos el data frame (con fortify) que convierte los polígonos en puntos, les asigna el id (NOM_ENT) correspondiente, y asigna también un grupo que indica a que polígono corresponde cada punto.
edo_df <- fortify(edo_shp, region = "NOM_ENT")
class(edo_df)## [1] "data.frame"head(edo_df)## long lat order hole piece group id
## 1 -102.1068 22.29674 1 FALSE 1 Aguascalientes.1 Aguascalientes
## 2 -102.1054 22.29639 2 FALSE 1 Aguascalientes.1 Aguascalientes
## 3 -102.0971 22.29499 3 FALSE 1 Aguascalientes.1 Aguascalientes
## 4 -102.0768 22.29267 4 FALSE 1 Aguascalientes.1 Aguascalientes
## 5 -102.0649 22.28985 5 FALSE 1 Aguascalientes.1 Aguascalientes
## 6 -102.0573 22.28697 6 FALSE 1 Aguascalientes.1 AguascalientesYa estamos listos para graficar, usaremos la geometría polígono.
ggplot(data = edo_df, aes(long, lat, group=group)) +
geom_polygon(colour='darkgrey', fill='white') +
coord_map(projection="mercator") +
theme_nothing()
En la instrucción de arriba el mapa se genera uniendo los polígonos (geom_polygon), las siguientes dos líneas son opcionales:
coord_map asegura que el cociente entre el ancho y alto del gráfico sea el correspondiente a la proyección de los datos y así evitar deformaciones.
theme_nothing limpia el gráfico eliminando grids y etiquetas.
Añadir variables al mapa
Nuestro objetivo final es hacer mapas para representar una variable, veamos como haríamos para representar el índice de carencias, el índice objetivo esta almacenado en el archivo indice_carencias.csv.
library(Hmisc)## Loading required package: lattice
## Loading required package: survival
## Loading required package: Formula
##
## Attaching package: 'Hmisc'
##
## The following object is masked from 'package:maptools':
##
## label
##
## The following objects are masked from 'package:dplyr':
##
## combine, src, summarize
##
## The following objects are masked from 'package:plyr':
##
## is.discrete, summarize
##
## The following objects are masked from 'package:base':
##
## format.pval, round.POSIXt, trunc.POSIXt, units# leer base de datos a representar
conapo <- read.csv("data/indice_carencias.csv", stringsAsFactors=FALSE)
# leemos la clave como caracter (es una manera redundante)
conapo$CVE <- read.table("data/conapo_2010.csv", header=TRUE, quote="\"",
stringsAsFactors=FALSE, colClasses = "character")$CVE
# Creamos una variable categórica (en este ejemplo queremos una cont y una cat.)
conapo <- mutate(conapo, indice_cat = cut2(PC1, 5))
# seleccionamos las variables que necesitamos para tener una base más limpia
# CVE es la clave del municipio y la usaremos para unir la variable
# a graficar con el mapa
indice <- select(conapo, indice = PC1, indice_cat, CVE)
# Filtramos las regiones de interés (DF, edo. de México y Morelos)
centro <- indice$CVE[conapo$NOM_ABR %in% c("DF", "Mex.", "Mor.")]
# Leemos los shapes de municipio
mun_shp <- readOGR("shapes" , "Municipios")## OGR data source with driver: ESRI Shapefile
## Source: "shapes", layer: "Municipios"
## with 2456 features
## It has 6 fieldsnames(mun_shp)## [1] "OBJECTID" "CVE_ENT" "CVE_MUN" "NOM_MUN" "Shape_Leng"
## [6] "Shape_Area"mun_shp@data$CVE <- paste(as.character(mun_shp@data$CVE_ENT),
as.character(mun_shp@data$CVE_MUN), sep = "")
mun_df <- fortify(mun_shp, region = "CVE")Como vemos en el código de arriba, para incluir variables en el mapa las añadimos a la base de datos mun_df. Para este ejemplo graficaremos únicamente una región. Podemos crear el subconjunto directamente del objeto SpatialPolygonsDataFrame (mun_df) pero en este ejemplo filtramos la base de datos data.frame ().
# unimos los polígonos (en data.frame) y las variables a graficar
mun_ind <- mun_df %>%
mutate(CVE = id) %>%
left_join(indice) %>%
filter(CVE %in% centro)## Joining by: "CVE"ggplot() +
geom_polygon(data = mun_ind, aes(long, lat, group = group, fill = indice)) +
labs(title = "Índice de carencias", fill = "Índice") +
coord_fixed()
En el siguiente mapa los colores son a nivel municipio pero dibujamos las fronteras a nivel estado para evitar que los bordes opaquen los datos.
# los municipios están en una proyección distinta a los estados por lo que
# si queremos que coincidan hay que reproyectarlos
mun_shp_2 <- spTransform(mun_shp, CRS("+proj=longlat +datum=WGS84"))
mun_df <- fortify(mun_shp_2, region = "CVE")
mun_ind <- mun_df %>%
mutate(CVE = id) %>%
left_join(indice)## Joining by: "CVE"ggplot() +
geom_polygon(data = mun_ind, aes(long, lat, group = group, fill = indice)) +
geom_polygon(data = edo_df, aes(x = long, y = lat, group = group),
fill = NA, color = "darkgray", size = 0.25) +
labs(title = "Índice de carencias", fill = "Índice") +
theme_nothing(legend = TRUE) + #fondo blanco
guides(fill = guide_legend(reverse = TRUE)) +
scale_fill_distiller(palette = "GnBu", breaks = pretty_breaks(n = 10)) + #paleta
coord_map()
En la instrucción de arriba, a partir de la cuarta línea los comandos son opcionales:
labs: título y etiqueta del gráfico.
theme_nothing(legend=TRUE): limpiar el gráfico (quitar grids, escalas, ejes,…) dejando la etiqueta.
guides(fill = guide_legend(reverse = TRUE)): cambiar las leyendas para que el valor más bajo quede en la parte inferior.
scale_fill_distiller: seleccionar la paleta de color que va de verde a azúl, categorías limpias.
Ahora, para usar mapas de ggplot como fondo usamos el paquete ggmap es importante considerar que estos mapas usan elipsoide World Geodesic Siystem 1984 y longitud/latitud.
map_Cuer <- get_map("Cuernavaca", maptype = "roadmap")## Map from URL : http://maps.googleapis.com/maps/api/staticmap?center=Cuernavaca&zoom=10&size=640x640&scale=2&maptype=roadmap&language=en-EN&sensor=false
## Information from URL : http://maps.googleapis.com/maps/api/geocode/json?address=Cuernavaca&sensor=falseggmap(map_Cuer)
Veamos un ejemplo de Ryan Peek
hotsprings <- read.csv("data/hotspringsCA.csv")
# Solo nos interesan algunas columnas
df <- hotsprings[, c(2:6, 10, 11)]
head(df)## LATITUDE LONGITUDE SpringName TF TC NOAA
## 1 39.022 -122.444 ABBOTT MINE 86 30 UKIAH
## 2 38.654 -122.484 AETNA SPRINGS 91 33 SANTA ROSA
## 3 34.538 -119.560 AGUA CALIENTE SPRING 133 56 LOS ANGELES
## 4 32.947 -116.305 AGUA CALIENTE SPRINGS 101 38 SAN DIEGO
## 5 33.363 -117.017 AGUA TIBIA SPRING 92 33 SANTA ANA
## 6 37.503 -121.904 ALAMEDA WARM SPRINGS 80 27 SAN JOSE
## USGS.Quadrangle
## 1 WILBUR SPRINGS 15
## 2 AETNA SPRINGS 7.5
## 3 HILDRETH PEAK 7.5
## 4 AGUA CALIENTE 7.5
## 5 (PALA 7.5)
## 6 NILES 7.5Podemos seleccionar la región usando lat/lon.
ca <- get_map(location = c(lon = -120.5, lat = 37.5), zoom = 6, crop = T,
scale = "auto", color = "color", maptype = "satellite") ## Map from URL : http://maps.googleapis.com/maps/api/staticmap?center=37.5,-120.5&zoom=6&size=640x640&scale=2&maptype=satellite&language=en-EN&sensor=falsegg <- ggmap(ca, extent = "panel", padding = 0)
gg_ca <- gg +
geom_point(data = df, aes(x = LONGITUDE, y = LATITUDE), size = 4,
pch = 21, fill = "orange", alpha = 0.6)
grid.arrange(gg, gg_ca, ncol = 2, main = "Ubicación de aguas termales")
Grafica los polígonos (de los shapes de Municipios) de tu estado (o DF) de nacimiento sobre un google maps.
También hay herramientas para hacer mapas dinámicos.
library(googleVis)
# creamos una columna latitiud-longitud que escribimos como la primera del df
df$LatLong <- paste(round(df$LATITUDE, 3), round(df$LONGITUDE, 3), sep = ":")
head(df)## LATITUDE LONGITUDE SpringName TF TC NOAA
## 1 39.022 -122.444 ABBOTT MINE 86 30 UKIAH
## 2 38.654 -122.484 AETNA SPRINGS 91 33 SANTA ROSA
## 3 34.538 -119.560 AGUA CALIENTE SPRING 133 56 LOS ANGELES
## 4 32.947 -116.305 AGUA CALIENTE SPRINGS 101 38 SAN DIEGO
## 5 33.363 -117.017 AGUA TIBIA SPRING 92 33 SANTA ANA
## 6 37.503 -121.904 ALAMEDA WARM SPRINGS 80 27 SAN JOSE
## USGS.Quadrangle LatLong
## 1 WILBUR SPRINGS 15 39.022:-122.444
## 2 AETNA SPRINGS 7.5 38.654:-122.484
## 3 HILDRETH PEAK 7.5 34.538:-119.56
## 4 AGUA CALIENTE 7.5 32.947:-116.305
## 5 (PALA 7.5) 33.363:-117.017
## 6 NILES 7.5 37.503:-121.904dfvis <- df[, c(8, 1:7)]
# Creamos una variable donde se almacena la información a desplegar
dfvis$Tip <- as.character(paste("Hot Spring =", dfvis$SpringName, "<BR>",
"Water Temp (F) =", dfvis$TF, " <BR>", "NOAA =", dfvis$NOAA, " <BR>",
"USGS Quad =", dfvis$USGS.Quadrangle, "<BR>"))
# Grafica
CAsprings <- gvisMap(dfvis, locationvar = "LatLong", tipvar = "Tip",
options = list(enableScrollWheel = TRUE,
height = 600, mapType = c("terrain", "satellite"), useMapTypeControl = TRUE))
#print(CAsprings) # abre en un nuevo htmlMás de googleVis en su github, y en el blog de Markus Gesmann.
Otros paquetes populares son RGoogleMaps y plotKML.
Leaflet Mapas interactivos, hay paquete para controlar mapas de Leaflet desde R.