sábado, 16 de junio de 2012

1/1000 Maneras de Graficar un Mapa en R

Que tal chibchombianos!

Así como existen 1000 maneras de morir tambien podemos hablar de que existen 1000 formas de hacer todo en R y en este caso graficar un mapa no es la excepción. Aunque personalmente me tomo tiempo encontrar una forma rápida y sencilla para lograrlo, aqui les ahorro esfuerzo y sufrimiento!...

Simplemente deben seguir las siguientes pasos para ser felices:

1. Descargar el shapefile.
El shapefile es un formato para sistemas de información geográfica. De acuerdo con la Wiki es un formato vectorial de almacenamiento digital donde se guarda la localización de los elementos geográficos y los atributos asociados a ellos. Si bien la definición es básica, es increible la cantidad de información descriptiva de un mapa que se puede almacenar en un shapefile.

Existe un shapefile disponible casi que para cualquier región y es posible encontrar varios sitios en la web que permiten acceder gratuitamente a sus shapefiles, en este caso les recomiendo la base de datos del GADM a la cual podrán acceder con el siguiente link http://www.gadm.org/ 

Otros links de interés 

Para el ejemplo utilizaremos la base de datos del GADM, asi que vayan a la sección de descarga y seleccionen el mapa de Uzbekistan (de paso se enteran que existe este país) y descarguen el shapefile respectivo. 

1a. Instalen y carguen en R las bibliotecas maptools, sp y RColorBrewer

Para instalar los paquetes:

> install.packages("maptools",lib="/directorio")

Para cargar:

> library(maptools) 
# Provee herramientas para leer y manejar objetos espaciales
> library(sp) 
#Provee métodos para manejar datos espaciales
> library(RcolorBrewer) 
# Permite acceder a paletas de colores para cartografía
http://colorbrewer2.org/

1b. Ahora lean el shapefile y conviertanlo en un vector, en este caso lo llamé ohsi
> ohsi<-readShapeSpatial("UZB_adm1.shp") 
#Corresponde al mapa politico de Uzbekistan

1c. Comprueben que se cargo el mapa
> plot(ohsi)

1d. Pueden incluso revisar la estructura interna del vector 
> str(ohsi)
> ohsi@data
1.e. Asignen el vector como una variable
> mapa=ohsi

2. El siguiente paso es construir una tabla que contenga en una columna las localidades "LOCAL" y en la siguiente columna algún criterio para la selección de las localidades, en este caso el criterio fue "lugares que quiero visitar en Uzbekistan" y la llame "VISITA"  :P. Es importante guardar la tabla en formato csv o delimitado por espacios.

Aqui tienen la tabla
"lOCAL"    "VISITA"
"Andijon"    "NA"
"Bukhoro"    1
"Ferghana"    "NA"
"Jizzakh"    "NA"
"Karakalpakstan"    "NA"
"Kashkadarya"    "NA"
"Khorezm"    "NA"
"Namangan"    "NA"
"Navoi"    "NA"
"Samarkand"    1
"Sirdaryo"    "NA"
"Surkhandarya"    "NA"
"Tashkent City"    "NA"
"Tashkent"    1

La codificación corresponde a: NA=No aplica y 1=lugares que quiero visitar
2a. Ahora deben abrir el archivo de la tabla y crear otro vector  el cual llamé "yque"

> yque<-read.csv("tabla.csv",header =TRUE, sep = "\t")

2b. Revisen los datos de la columna visita en la tabla. Deberian ser visibles en la consola.

> yque$VISITA

2c. Reasignen la información de vector "yque" variable "VISITA" a "mapa" 
> mapa@data=data.frame(yque$VISITA)
> mapa

3. Aunque hay múltiples paquetes para asignar colores en un mapa, en este caso me base en la paleta de colores de colorBrewer. El color pueden seleccionarlo guiándose con el siguiente link http://colorbrewer2.org/

Ahora grafiquen el mapa con las localidades seleccionadas en "VISITA" y el color lo deben especificar en col.regions=brewer.pal() como sigue a continuación:


> spplot(mapa,c("yque.VISITA"),col.regions=brewer.pal(3, "RdYlBu"), scales=list(draw = TRUE))

Listo!! ahora tenemos el mapa de Uzbekistan y las áreas seleccionadas corresponden a las localidades que me gustarian visitar :P.

otro ejemplo


Fácil! Ahora lo pueden modicar y adaptar a sus datos. Esta es solo una forma sencilla, les recuerdo hay cientos de opciones por explorar!

Buena energia!










10 comentarios:

  1. y si tengo la coordenada y no se el pa'is ¿qu'e hago? obviamente desde R

    ResponderEliminar
  2. si tenes la coordenada(osea puntos)podes llamar un mapa mundi (o alguna region en particular) y poner tus puntos:

    library(maptools)
    data(wrld_simpl)

    plot(wrld_simpl,axes=T,col='lightgrey',xlim=c(-180,20), ylim=c(-50,11))# con xlim y ylim le das un "zoom" a la región deseada


    plot(wrld_simpl,add=TRUE, border='dark grey',xlim=c(-80,10), ylim=c(-60,10))
    points(x=-73.153477,y=6.015178,col=2,pch=16, cex=0.3)

    points(x=c(-103.214156,-102.879113),y=c(28.073906,21.220430), col=1,pch=22, cex=0.6)
    points(x=c(-73.666888,-73.924321),y=c(6.723821,3.969741), col=2,pch=20, cex=0.6)
    points(x=c(77.537207,77.482795),y=c(23.517155,17.222947), col=3,pch=20, cex=0.6)
    points(x=c(114.872530),y=c(4.269109), col=4,pch=20, cex=0.6)
    points(x=c(-86.265691),y=c( 15.199517), col=5,pch=20, cex=0.6)
    points(x=c(105.061967),y=c(12.565841), col=6,pch=20, cex=0.6)
    points(x=c(102.599696),y=c(19.557999), col=7,pch=20, cex=0.6)
    points(x=c( 46.371590),y=c(22.480598), col="lightgreen",pch=20, cex=0.6)
    points(x=c(8.709582),y=c(9.090636), col="orange",pch=20, cex=0.6)
    points(x=c(108.190736),y=c( 14.055180), col="purple",pch=20, cex=0.6)

    ResponderEliminar
    Respuestas
    1. ¿hay manera de cambiar los símbolos por números?

      Eliminar
  3. excelente post, lo realice con el mapa de Argentina pero por algun motivo me dio error en el data.frame:
    Error en data.frame(..., check.names = FALSE) :
    arguments imply differing number of rows: 503, 9

    el codigo que utilice, adaptando del ejemplo, es el siguiente:
    # carga de mapa politico e impresion del mismo
    cartografia<-readShapeSpatial("d:/mapas/ARG_adm2.shp")
    plot(cartografia)
    # asigno vector con referencias al mapa
    str(cartografia)
    cartografia@data
    mapa=cartografia

    # leer datos de una tabla
    datos_mani<-read.csv("d:/mapas/tabla.csv",header =TRUE, sep = "\t")
    datos_mani
    # reasignar informacion al mapa y vector referencias, ojo al PIOJO en el vector!
    mapa@data=data.frame(datos_mani$LOCAL, datos_mani$MANI, check.names=TRUE)

    # colorear mapa
    spplot(mapa,c("datos_mani.MANI"),col.regions=brewer.pal(3, "RdYlBu"), scales=list(draw = TRUE))

    ¿Cual podra ser el error? Agradeceria una orientacion

    ResponderEliminar
  4. Hola! :)
    Asegúrate de tener la misma cantidad de columnas en el archivo csv y que la columna que contiene las provincias o regiones de Argentina corresponda a el numero real en el shapefile.

    Para ello con cartografia@data buscas las columna y la copias tal cual.

    2) Aquí "mapa@data=data.frame(datos_mani$LOCAL, datos_mani$MANI, check.names=TRUE)"
    Sólo debes indicar "una" columna y esta debe ser la correspondiente a
    la que contiene los valores para asignar color a la columna de las provincias.
    por ejemplo
    "mapa@data=data.frame(datos_mani$MANI, check.names=TRUE)"

    Espero te ayude ;)

    ResponderEliminar
    Respuestas
    1. Este comentario ha sido eliminado por el autor.

      Eliminar
    2. (corrijo la respuesta anterior) Gracias Susana, el problema era tal cual mencionas en tu respuesta, ademas tenia un error en la tabla de datos (menos filas de datos que las de la base del SHP) por lo que corregi la tabla y ahora me funciona perfectamente.
      Aprovecho para preguntarte si hay manera de tomar lso nombres de las provincias del ejemplo para imprimir dentro de cada poligono, intente pero no logro comprenderlo. Gracias nuevamente y un cordial saludo desde Corrientes/Argentina

      Eliminar
  5. Holaa yo estoy intentando hacer el mapa de china pero no me sale ... al leer el archivo shp me sale el siguiente error:

    Error in getinfo.shape(fn) : Error opening SHP file

    Que hago !!!!!

    ResponderEliminar
  6. Que genial tu aporte ... Muchisimas gracias, sin embargo tengo un problema ... tengo el shape de municipios de Colombia descargado del sigot, hago todo hasta el momento de darle color, pero cuando termino en mi mapa no están coloreadas las áreas sino las lineas de división municipal ... alguna sugerencia??

    ResponderEliminar
  7. Hola

    Estuve leyendo mi archivo shape pero presenta este Error in ogrListLayers(dsn = dsn) : Cannot open data source, ¿alguien sabe por que se presenta este error?

    ResponderEliminar