Desde Escuela de Datos, Sebastián Oliva, fellow 2017, enseña cómo usar Python para generar mapas a partir de datos georreferenciados.

Mapas y Python

Es obvia la importancia de los mapas, para la visualizacion de datos. Las coordenadas, latitud y longitud, pueden describir un punto sobre la tierra. Utilizamos estandares como WGS-84 para atar esas coordenadas a un punto real.

Utilizando MatPlotlib, podemos aprovechar Basemap, una libreria que provee funcionalidad básica de mapa, con la cual podemos construir y componer. Agregar poligonos, puntos, areas, barras, colores, etc; se hace mediante estas librerias.

Librerias

La libreria mas utilizada en el ecosistema Jupyter-Matplotlib es Basemap. Tambien existen otras, entre ellas, Plotly, que son muy poderosas y convenientes pero tienen dependencias externas.

• Basemap: el «industry standard». Un poco complicado para el setup, bastante poderoso e integrado con matplotlib.
• Plotly: Una libreria que permite desplegar datos mediante Javascript en la Web
• QGis: Una aplicación completa para el manejo de datos GIS, tiene bindings en Python y es posible utilizar como libreria para aplicaciónes tanto GUI como en notebooks o para analizar en scripts y exportar.
• OSMnx: Basada en Open Street Map, Permite tanto analisis como visualización de mapas a nivel de calle, region, ciudad y más.

Para este webinar, vamos a indagar mas en Basemap, que es el mas accesible.

Basemap

Basemap es una extensión de la funcionalidad disponible

Existen varias formas de instalarlo, así que puede ser un poco confuso. Dependiendo de el método en el cual tengas instalado matplotlib hace variar la forma apropiada de instalarlo.

Ambiente de Trabajo

$ #ESTE_ENV = midevenviroment
$ source ~/miniconda3/envs/$ESTE_ENV/bin/activate

$ conda install jupyter-notebook
$ conda install gdal -c conda-forge
$ conda install basemap -c conda-forge
$ conda install pandas seaborn

## En caso hayan instalado basemap en algun directorio no standard: utiliza un link para la carpeta data.
$ ln -s /home/tian/miniconda3/pkgs/basemap-1.1.0-py36_2/lib/python3.6/site-packages/mpl_toolkits/basemap/data/ /usr/share/basemap

In [16]:

# Importamos lo ya usual.
import matplotlib.pyplot as plt
import matplotlib.cm
import pandas as pd
import numpy as np
import seaborn

# Algunas librerias extra que usaremos
from matplotlib.colors import Normalize
import matplotlib.colors as colors
from numpy import array
from numpy import max

# Aqui cargamos Basemap
from mpl_toolkits.basemap import Basemap
from matplotlib.patches import Polygon
from matplotlib.collections import PatchCollection


sns.set(style="white", color_codes=True)

%matplotlib inline

Mapeando los terremotos globales de la ultima semana

Vamos a usar la feed de datos del US Geological, ellas tienen disponibles datos referenciados de actividad geologica a nivel mundial, regional y de EEUU.

In [17]:

quakes = pd.read_csv("http://earthquake.usgs.gov/earthquakes/feed/v1.0/summary/1.0_week.csv")

# Creamos la lista de latitudes y longitudes.
lats, lons = list(quakes['latitude']), list(quakes['longitude'])

In [60]:

mags = list(quakes['mag'])
quakes.head()

Out[60]:

5 rows × 22 columns

Iniciemos con el mapa

In [25]:

eq_map = Basemap(projection='robin', resolution = 'l', area_thresh = 1000.0,
              lat_0=0, lon_0=-130)
eq_map.drawcoastlines()
eq_map.drawcountries()
eq_map.fillcontinents(color = 'gray')
eq_map.drawmapboundary()
plt.show()
/usr/lib64/python3.6/site-packages/mpl_toolkits/basemap/__init__.py:1631: MatplotlibDeprecationWarning: The get_axis_bgcolor function was deprecated in version 2.0. Use get_facecolor instead.
  fill_color = ax.get_axis_bgcolor()
/usr/lib64/python3.6/site-packages/mpl_toolkits/basemap/__init__.py:1775: MatplotlibDeprecationWarning: The get_axis_bgcolor function was deprecated in version 2.0. Use get_facecolor instead.
  axisbgc = ax.get_axis_bgcolor()
 

In [64]:

figu, ax = plt.subplots(figsize=(20,10))
eq_map = Basemap(projection='robin', resolution = 'l', area_thresh = 1000.0,
              lat_0=0, lon_0=-130)
eq_map.drawcoastlines()
eq_map.drawcountries()
eq_map.fillcontinents(color = 'gray')
eq_map.drawmapboundary()
eq_map.drawmeridians(np.arange(0, 360, 30))
eq_map.drawparallels(np.arange(-90, 90, 30))

## Coordenadas a posiciones
x,y = eq_map(lons, lats)
eq_map.plot(x, y, '^', markersize=6)
 
plt.show()
/usr/lib64/python3.6/site-packages/mpl_toolkits/basemap/__init__.py:1631: MatplotlibDeprecationWarning: The get_axis_bgcolor function was deprecated in version 2.0. Use get_facecolor instead.
  fill_color = ax.get_axis_bgcolor()
/usr/lib64/python3.6/site-packages/mpl_toolkits/basemap/__init__.py:1775: MatplotlibDeprecationWarning: The get_axis_bgcolor function was deprecated in version 2.0. Use get_facecolor instead.
  axisbgc = ax.get_axis_bgcolor()
/usr/lib64/python3.6/site-packages/mpl_toolkits/basemap/__init__.py:3298: MatplotlibDeprecationWarning: The ishold function was deprecated in version 2.0.
  b = ax.ishold()
/usr/lib64/python3.6/site-packages/mpl_toolkits/basemap/__init__.py:3307: MatplotlibDeprecationWarning: axes.hold is deprecated.
    See the API Changes document (http://matplotlib.org/api/api_changes.html)
    for more details.
  ax.hold(b)

In [33]:

x[0]

Out[33]:

17740352.000926033

Veamos esto mas a detalle.

In [28]:

df = pd.read_csv('hiv_cr_data.csv')
df.columns
df.shape # (71, 8)
df.describe()
df.columns
df.loc[df.coordenadas == df.coordenadas]
subset = df.loc[df.coordenadas == df.coordenadas] 
coordenadas = subset[['sitio','latitud', 'longitud', 'coordenadas']]
coordenadas.head()

Out[28]:

In [66]:

coordenadas.count()

Out[66]:

sitio          67
latitud        67
longitud       67
coordenadas    67
dtype: int64

In [67]:

coordenadas.coordenadas.head()

Out[67]:

0      (10.01801 , -84.21648)
1    (10.016787 , -84.213914)
2    (10.020168 , -84.214064)
3      (10.00602 , -84.21274)
4       (9.98512 , -84.19554)
Name: coordenadas, dtype: object

In [35]:

 (10.01801 , -84.21648)

Out[35]:

(10.01801, -84.21648)

In [30]:

fig, ax = plt.subplots(figsize=(10,20))

mapa = Basemap(projection='merc', 
    lat_0 = 9.74, lon_0 = -83.5,
    resolution = 'i',
    llcrnrlon=-88.1, llcrnrlat=5.5,
    urcrnrlon=-80.1, urcrnrlat=11.8)

mapa.drawmapboundary(fill_color='#479EE0')
mapa.drawcoastlines()


from ast import literal_eval as make_tuple
def unpac(t):
    # haciendo trampa en la vida
    return pd.Series(make_tuple(t))

def plot_area(pos):
    ps = unpac(pos)
    x, y = mapa(ps[1], ps[0])
    mapa.plot(x, y, 'o', markersize=7, color='#444444', alpha=0.8)
    
    
coordenadas.coordenadas.apply(plot_area)


plt.show()
/usr/lib64/python3.6/site-packages/mpl_toolkits/basemap/__init__.py:3298: MatplotlibDeprecationWarning: The ishold function was deprecated in version 2.0.
  b = ax.ishold()
/usr/lib64/python3.6/site-packages/mpl_toolkits/basemap/__init__.py:3307: MatplotlibDeprecationWarning: axes.hold is deprecated.
    See the API Changes document (http://matplotlib.org/api/api_changes.html)
    for more details.
  ax.hold(b)

In [12]:

clox = array(coordenadas['longitud'])
clay = array(coordenadas['latitud'])
clo = list(clox)
cla = list(clay)

In [51]:

clay.mean()
clay

Out[51]:

array([ 10.01801 ,  10.016787,  10.020168,  10.00602 ,   9.98512 ,
        10.001528,   9.998438,   9.99943 ,   9.998952,   9.996179,
         9.98495 ,   9.99961 ,   9.935734,   9.93335 ,   9.93284 ,
         9.93355 ,   9.9356  ,   9.9359  ,   9.93454 ,   9.927243,
         9.93387 ,   9.93191 ,   9.93378 ,   9.937275,   9.937206,
         9.93281 ,   9.868255,   9.864336,   9.864255,   9.86715 ,
         9.97685 ,   9.99725 ,   9.974695,   9.61626 ,   9.39646 ,
         9.42387 ,   9.43062 ,   9.930423,   9.930036,   9.934636,
         9.929361,   9.937733,   9.930169,   9.927714,   9.934579,
         9.927496,   9.93141 ,   9.938098,   9.927755,   9.933922,
         9.936659,   9.932065,   9.927739,   9.930635,   9.932147,
         9.93535 ,   9.93286 ,   9.927324,  10.018506,  10.018993,
        10.002973,   9.408455,   9.39838 ,   9.403425,   9.40677 ,
         9.866258,   9.865848])

In [58]:

plt.figure(2)
#fig.add_subplot(223)
fig2, ax2 = plt.subplots(figsize=(20,10))

mapa2 = Basemap(projection='merc', 
    lat_0 = 9.74, lon_0 = -83.5,
    resolution = 'i',
    llcrnrlon=-88.1, llcrnrlat=7.5,
    urcrnrlon=-80.1, urcrnrlat=11.8)
pos_x, pos_y = mapa2(clox, clay)
mapa2.drawmapboundary(fill_color='#A6CAE0', linewidth=0)
mapa2.fillcontinents(color='darkgrey', alpha=0.3)
mapa2.drawcoastlines(linewidth=0.1, color="white")

paleta = seaborn.diverging_palette(10, 220, sep=80, as_cmap=True)

#sns.cubehelix_palette(8, start=2, rot=0, dark=0, light=.95, reverse=True)

hb = plt.hexbin(pos_x, pos_y, gridsize=4, mincnt=1,
    edgecolor='none', cmap = paleta)

cb = fig2.colorbar(hb, ax=ax2)


plt.show()
/usr/lib64/python3.6/site-packages/mpl_toolkits/basemap/__init__.py:1775: MatplotlibDeprecationWarning: The get_axis_bgcolor function was deprecated in version 2.0. Use get_facecolor instead.
  axisbgc = ax.get_axis_bgcolor()
<matplotlib.figure.Figure at 0x7f8646f5d7f0>

Cool links

• http://beneathdata.com/how-to/visualizing-my-location-history/
• https://pypi.python.org/pypi/descartes *

Desde Bogotá D.C., Andrés Forero  , especialista en GIS y datos geográficos, nos envía este post como parte de una serie de tutoriales  sobre datos geográficos y cómo podemos aprovecharlos desde distintos software para su edición y publicación.

En estos últimos años he tenido la oportunidad de ejercer diferentes roles como consultor y asesor para el desarrollo territorial en Colombia, un rol que, dadas las deficiencias tecnológicas y gubernamentales existentes, ha sido un enorme reto profesional.

En Colombia existen 1.123 municipios, 900 de ellos requieren la actualización de su Plan de Ordenamiento Territorial, el cual es el reglamento oficial que regula y genera normas para el desarrollo adecuado del territorio y la utilización del suelo. Solo 100 de ellos, de acuerdo con el Departamento Nacional de Planeación – DNP, han apostado por modernizarse, concentrados en resolver grandes desafíos asociados a sistemas viales, transporte multimodal, gestión del riesgo de desastres, cambio climático y protección de recursos naturales, entre otros aspectos necesarios para el logro de una mejor calidad de vida y de ciudades sostenibles al servicio de sus habitantes.

Empoderar al ciudadano, al especialista y al tomador de decisiones implica encontrar herramientas que integren soluciones para dichos desafíos donde se garantice la participación conjunta. Partiendo de esta premisa, la construcción que hemos empezado a generar nace desde la cartografía social donde reunidos con líderes comunitarios y habitantes de un territorio se han empezado a identificar todos aquellos actores involucrados, desde la espacialización de los puntos de interés hasta la identificación de estrategias llevadas a cabo para mitigar riesgos.

Estas iniciativas han trascendido de lo análogo a lo digital, y es apenas lógico, es un salto necesario y obligatorio el cual facilita una toma de decisiones ideal. En Colombia, por ejemplo, la Organización Internacional para las Migraciones – OIM  ha realizado diferentes proyectos desde la cartografía social, a lo largo del territorio nacional enfocados en comunidades afectadas por el conflicto armado y la formulación de estrategias aplicadas al posconflicto. Desde septiembre del 2018, se ven reflejados en una biblioteca virtual georeferenciada conocida como BiblioDat, la biblioteca sobre construcción de paz y el conflicto en Colombia. México es otro gran ejemplo del empoderamiento de la tecnología y generación de información comunitaria en tiempo real a través de mapas colaborativos para la identificación con precisión de infraestructura, centros de acopio y capas de información posdesastre.

Mapas táctiles para personas con discapacidad visual. (Foto: Andrés Forero)

¡Empecemos!

Creando mapas utilizando ArcGIS Online

 ArcGIS Online es un sistema SIG en la web de colaboración online que permite usar, crear y compartir mapas, escenas, aplicaciones y datos listo para usar. Para poder acceder a sus funcionalidades es necesario entrar allí creando una cuenta, de tal manera, que puedas tener un portal SIG propio. Sí, para nadie es un secreto que Esri fabrica software licenciado y eso no lo pondremos en discusión, pero déjenme decirles que como desarrolladores podemos acceder a sus funcionalidades para poderlas utilizar con estos fines… ¡No tiene costo!

Al momento de seguir los pasos que te solicitan allí, como el llenado de datos como nombre, apellido, correo, etc., te pedirá que le asignes un nombre al portal SIG que acabas de crear, el cual tendrá como dominio el siguiente:

http://tuportal.maps.arcgis.com

Listo, tienes tu portal geográfico desde el cual accederás con un usuario y contraseña donde podrás crear los mapas, publicar información y disponer la información que elijas hacia la comunidad.

Primero, ¿qué queremos mostrar en el mapa?

ArcGIS for Developers brinda la facilidad de agregar contenido proveniente de diferentes fuentes de información, o crearla desde allí, lo cual puede ser mucho mejor. Aquí empieza el trabajo de ustedes, dateros, la información espacial también mantiene unas buenas prácticas…. ¡y vaya si son extensas! dan para todo un estudio universitario de pregrado y hasta posgrado, de esto depende poder tener una representación geográfica adecuada. Comencemos creando una capa geográfica definiendo algunos parámetros utilizando la opción «Create an Empty Layer».

La primera capa geográfica que crearemos corresponde a la identificación de (algunos) puntos de interés al interior de una zona de estudio. Villa de Leyva es un municipio colombiano ubicado en la provincia de Ricaurte del departamento de Boyacá, está ubicado a 40 km al occidente de Tunja, capital del departamento. Sobre este municipio agregaremos inicialmente tres puntos: la Plaza Mayor, el Museo Paleontológico y Pozos Azules. Para poder hacerlo, esta capa la denominaremos “puntos_de_interes” el cual será su título y nombre del servicio disponible para ser utilizado por otros desarrolladores. Agregamos algunas etiquetas para identificar con facilidad esta capa, en este caso: puntos de interés, Villa de Leyva, Escuela de Datos.

Detalles de la capa geográfica

Al hacer click sobre “Geometry” definimos cómo queremos representar la información, sea como punto, línea o polígono y bajo qué sistema de referencia estará dicha información. Aquí definimos “punto” y el valor por defecto “4326” el cual corresponde al sistema de referencia WGS 84.

Geometría de la capa geográfica

Luego definimos los campos que deseamos incluir sobre la capa; se incluirán el nombre del campo, el alias y el tipo de dato (string, integer, double, date). Agregaremos cuatro de ellos los cuales se nombrarán así:

Campos de la capa geográfica

Por último, se realizan las configuraciones necesarias para hacer de la capa un elemento editable, compartido y actualizado (offline y online) las cual se desplegarán por defecto. Para este caso, la capa será compartida para el público. Haz click en “Create layer” y ¡voilá!, se ha creado una capa geográfica.

La ventana que aparecerá por defecto corresponde a los metadatos de la capa geográfica la cual está sujeta a modificaciones como:

• Imagen miniatura
• Resumen
• Descripción
• Fecha de creación
• Fecha de modificación
• Etiquetas
• Capas incluidas
• Previsualización
• Créditos
• Detalles
• Acceso y restricciones de uso
• URL del servicio

Si te das cuenta, desde este momento ArcGIS establece como prioridad describir la capa que se está creando y es una buena práctica que nosotros debemos empezar a implementar en nuestro flujo de trabajo. Aclaro que es opcional, ya que este proceso puedes hacerlo luego de agregar los datos.

 Descripción de la capa geográfica

Representando datos sobre un mapa

Ya van 1110 palabras y agradezco que hayas llegado hasta esta sección. Aquí empezarás a agregar información en ArcGIS Online con la capa creada online, también puedes agregarla desde distintas fuentes de información.

 Descripción de la capa geográfica

Al hacer click sobre “Open in Map Viewer” verás el que será tu espacio de trabajo desde ahora, aquí desplegarás, procesarás y publicarás la información. Como ves, el espacio por defecto es realmente simple y facilita realizar algunas configuraciones.

1. Primero cambia el mapa base. Selecciona la pestaña “Basemap” y elige “Open StreetMap”
2. En la barra de búsqueda escribe “Plaza Mayor, Villa de Leyva”
3. Ahora haz click sobre la pestaña “Edit” y “New Feature” y agregarás la información pedida en los campos creados anteriormente

Agregando información al mapa

Los datos que se agregarán son los siguientes:

Visualización sobre el mapa

Puedes agregar los puntos que necesites de acuerdo a tus requerimientos, esta es la forma en que visualizar la información a través de ventanas emergentes. En este enlace encontrarás más información sobre cómo configurar el pop-up. Luego de editarla se desplegará así:

Ventana emergente sobre el mapa

Agregando información desde un archivo separado por comas (.csv)

 Como todos en todos los archivos, aquí se debe tener un trabajo de carpintería por debajo para desplegar correctamente la información; es esencial que contenga entre sus campos, datos geográficos los cuales estén relacionados con coordenadas (X, Y) o una dirección que pueda geocodificar. Para desplegarla únicamente tengo que arrastrar el archivo .CSV el cual contiene más puntos de interés hacia el mapa. Puedes descargarla haciendo click aquí.

Cargue de información desde un .CSV

Agregando información desde un shapefile (.shp)

El shapefile es el formato de datos vectorial más popular y extendido en el trabajo con un SIG, se trata de un formato vectorial que guarda la localización de elementos geográficos y atributos a ellos asociados, pero no es capaz de almacenar relaciones topológicas. Estos elementos geográficos se pueden representar a partir de una capa de tipo punto, línea o polígono (áreas).

Un shapefile se compone de varios archivos informáticos, sin los cuales no podríamos visualizarlo en un software SIG, siendo tres el número imprescindible:

• .shp – almacena las entidades geométricas de los objetos.
• .shx – almacena el índice de los elementos geométricos.
• .dbf – tabla dBASE donde se almacenan los atributos de los elementos geométricos.

Por otro lado y, opcionalmente, un shapefile puede estar compuesto por otros archivos que lo hacen más completo para la explotación de la información, fundamentalmente son:

• .prj – almacena la información relativa al Sistema de Coordenadas.
• .sbn y .sbx – almacena el índice espacial de las entidades.
• .shp.xml – almacena los metadatos de la capa.

Teniendo claro el concepto, lo que haremos ahora es agregar algunas vías que se encuentran allí, esto debe estar en un archivo comprimido .ZIP el cual será cargado en ArcGIS Online así:

1. Ve a la pestaña “Add” y selecciona “Add Layer from File”
2. Selecciona el archivo “vias.zip” el cual podrás descargar desde aquí

Carga de información como un archivo

Agregando información desde ArcGIS Online como un servicio

Miles de usuarios a diario publican información como un servicio, de tal manera que pueda ser consumido por otras personas. Esto facilita que no haya redundancia de información, que se pueda procesar datos actualizados y que pueda llegarse a utilizar información oficial.

Para este ejercicio llamaremos la división administrativa del municipio de Villa de Leyva únicamente buscando la capa que ya está creada.

1. Ve a la pestaña “Add” y selecciona “Search for Layers”
2. Busca “Veredas Villa de Leyva” en “ArcGIS Online” y verifica que la casilla “Within map area” esté deshabilitada

Carga la capa que aparece allí

Carga de información

Cambia el nombre de las capas

Viste algunas de las formas en que se agrega y se edita información a un conjunto de capas, ahora es necesario personalizar la visualización de tal manera que configurando una simbología adecuada llegue a ser mucho más amigable para el usuario final. Para ello es necesario modificar primero el nombre de las capas.

1. Dirígete a la pestaña “Details” y selecciona la capa “veredas villa de Leyva”.
2. En “More options” selecciona “Rename” y nómbralo ahora como “Veredas”.

Renombrando capas

Cambio de simbología sobre la capa “Vereda”

Personalizando simbología sobre la capa “Vereda”

De aquí para adelante eres tú quien decide cómo mostrarás estas capas, ¿cómo representarías adecuadamente los drenajes y los puntos de interés?

Guarda y comparte el mapa…

1. En la pestaña “Save” podrás asignarle un título, un resumen, etiquetas y guardar el mapa que has creado.
2. Una vez guardado, ve a la pestaña “Share” y elige la audiencia a quien mostrarás el mapa.

Edición final simbología

Hay más…

Puedes representar esta información de una manera mucho más amigable y te invito a que revises todos los storymaps que Esri ofrece donde existen diferentes plantillas de aplicaciones web responsive que pueden servirte para contrar más historias. Revisa https://storymaps.arcgis.com/en/ y escoge cómo quieres representar lo que quieres contar.