¿Por qué Python es el futuro del SIG web?

Tomado de:   Beatriz Ramos López

El mundo GIS presentó a Python como un lenguaje de scripting relativamente fácil, pero con el tiempo se ha convertido en omnipresente, ofreciendo soluciones para muchos usuarios diferentes. Gestión de datos, mapeo, análisis, administración del sistema, lo que usted llama: las posibilidades de Python son infinitas. Siga leyendo para averiguar cómo se puede utilizar a su ventaja y mejorar su carrera como un profesional geoespacial.

1.   Los lenguajes de programación en SIG están aquí para quedarse
En los últimos seis meses, se ha producido una adopción más amplia de los lenguajes de programación en la industria de SIG. Alguien en el DevSummit de Esri lo resumió perfectamente diciendo que Esri ha abrazado múltiples lenguajes de programación en un único flujo de trabajo. ¿Qué significa esto? Esto significa que usted puede hacer su trabajo con uno o varios lenguajes de programación (o ampliarlo si prefiere una GUI). Por supuesto, no tienes que hacer esto. Recientemente recibí un correo electrónico de alguien que me preguntaba acerca de las posibilidades de usar Python para hacer SIG y si realmente era necesario aprenderlo. La respuesta es no: Python no es un reemplazo para un SIG, sino una extensión. Y también lo son R y JavaScript. Sin embargo, muchos trabajos geoespaciales requieren que construya aplicaciones y / o domine varios lenguajes de programación.

2.   Python ahora trasciende la dicotomía de scripting / programación
Mi primer artículo sobre Python para este blog se llamó "Python programming: a beginners guide". Alguien me preguntó si el uso de herramientas preconfiguradas aplicadas a los datos GIS te convertía en un programador. Bastante, no lo creo. El uso de una herramienta de búfer en arcpy no te convierte en programador, pero si utilizas tal herramienta dentro de un script de aplicación más grande, esto califica como escribir una aplicación en lugar de un script. Y la creación de aplicaciones es lo que la programación se trata. Sin embargo, con respecto a Python como un lenguaje, la distinción entre el guion y la programación no siempre es muy clara.
Lo que es interesante, sin embargo, es que en los últimos seis meses, he empezado a usar herramientas Python que superan el GIS de escritorio (y por lo tanto la dicotomía de scripting / programación). Especialmente la API de Python de Esri combinada con los portátiles de Jupyter ha sido una gran nueva herramienta que es mucho más que una serie de herramientas de SIG preconfiguradas y ofrece una gran introducción al concepto de web de GIS de Esri, ya sea ArcGIS Online o Portal de ArcGIS. Lo que me lleva al siguiente punto:

3.   Python es ideal para operaciones de sistema, gestión de datos, análisis y visualización
Lo que quiero decir con esto es que puedes usar Python como lo harías con Window's File Explorer para administrar datos, cargarlo, procesarlo y visualizarlo como quieras. Después de trabajar con los portátiles Jupyter, que le permiten hacer todas estas operaciones a través de un navegador, me pregunté si se trataba de los nuevos IDE para escribir scripts y aplicaciones. Los portátiles Jupyter ofrecen capacidades impresionantes que superan a las de los IDEs, principalmente porque integran código con imágenes, como gráficos y, en el caso de Python API de Esri, mapas 2D y 3D. Son ideales para organizar y probar el código en trozos pequeños y hacer que los resultados se devuelvan al instante, literalmente en la misma página en la que escribes el código. Esto no es posible con un IDE, que devuelve los resultados de código en un shell y no se puede eliminar a menos que reinicie el IDE.

Añadir leyenda

4.   SIG ha abrazado a Python 3
Una diferencia importante entre arcpy para ArcMap y arcpy para ArcGIS Pro es la versión de Python utilizada en el fondo. La mayoría de las discusiones sobre las diferencias entre los dos ignoran el hecho de que Python 3 es realmente un rediseño de ciertas partes del lenguaje, por lo que tiene un número de versión diferente y se considera una actualización importante. Esto explica por qué un paquete de sitio como arcpy para ArcMap (que usa Python 2.7) tiene comandos y herramientas idiosincrásicas, mientras que la nueva API de Python para ArcGIS está estructurada y se lee más como Python puro, haciéndola más accesible para los desarrolladores de Python nuevos en GIS. Usar y aprender. Esto puede ayudar a Esri a atraer a más desarrolladores no GIS a usar su plataforma y herramientas.

Añadir leyenda

5.   SIG ha abrazado Anaconda
Seamos realistas, el control de versiones y la gestión de paquetes con Python es un desafío. Tener una aplicación que hace esto por ti te hace la vida más fácil y esto es lo que Anaconda hace. También le proporciona la aplicación Jupyter Notebook que se mencionó anteriormente. Mientras que al principio era sospechoso sobre su utilidad, ahora prefiero Anaconda a usar virtualenv cuando uso Python 3. También trabaja muy bien si usted tiene Python 2.7 instalado, que viene con ArcMap (el truco no es agregar Anaconda a usted PythonPath como no Para mezclar ambos entornos en la misma máquina). Esri también se dio cuenta de que sería más conveniente ofrecer arcpy a través de Anaconda en lugar de enviarlo directamente con ArcMap o Pro. Esto se hará cargo de las próximas versiones de software.

6.   SIG ha abrazado la pila SciPy
Esri ahora envía la pila SciPy con cada instalación de ArcMap y Pro para que usted no tenga que instalar los paquetes usted mismo. Aunque los paquetes individuales de esta pila están orientados hacia la comunidad científica, los analistas y cartógrafos de SIG pueden utilizarlos para su ventaja, por ejemplo, incluyendo parcelas simples de la biblioteca matplotlib. Si usas pandas y Numpy juntos, tienes algo similar a lo que ofrece el lenguaje R. Nuevamente, los portátiles Jupyter son el entorno perfecto para escribir todo su código, comentarios, resultados y gráficos.
En cuanto al marketing, es interesante ver que Esri ha adoptado el eslogan "La Ciencia de Dónde" para enfatizar sus raíces científicas. La colaboración con Continuum Analytics (la compañía detrás de Anaconda), el envío de la pila SciPy con ArcGIS y la creación del puente ArcGIS para usuarios de R son iniciativas recientes para unir a diferentes comunidades a través de la tecnología.

7.   GeoPython es un emocionante nuevo evento para desarrolladores geoespaciales de Python

GeoPython es una conferencia europea de Python para la comunidad geoespacial. Aquí, los nuevos proyectos de la comunidad de código abierto, así como los vendedores propietarios, se comparten. Si no visitó el evento, podría estar interesado en el repositorio Github que muestra todos los proyectos que se presentaron en la última edición, a principios de este año en mayo.

Añadir leyenda

8.   La fuente abierta amplía el GIS propietario y viceversa

Esto es principalmente una repetición de los dos últimos puntos, pero hay más: hay muchas herramientas de código abierto que amplían GIS comerciales. Muchos de éstos se distribuyen a través del índice del paquete de PyPI Python, aunque Github está tomando su lugar lentamente. Le sorprenderá la cantidad de resultados de búsqueda cuando busque "GIS", "Geo" o "esri". Un paquete de código abierto llamado ArcREST ofrece un conjunto de herramientas python para ayudar a trabajar con ArcGIS REST API para ArcGIS Server (AGS), ArcGIS Online (AGOL) y ArcGIS WebMap JSON. Además, si utiliza Anaconda como gestor de paquetes, puede buscar paquetes utilizando la nube Anaconda o navegar a través de los paquetes predeterminados de 1000+ que vienen con su instalación.

9.   Python habilita GIS Web
Aunque el GIS de escritorio no desaparecerá pronto, el SIG web es la ola del futuro. Python es una gran herramienta para crear flujos de trabajo usando Web GIS, que es otra razón para comenzar a aprenderlo. En lugar de depender de las aplicaciones de escritorio, puede reproducir gran parte de su funcionalidad fuera de dicha aplicación, por ejemplo, utilizando un cuaderno de Jupyter y aprovechar las funcionalidades, datos y herramientas que necesita utilizando las API web, crear sus propios flujos de trabajo y compartirlos con otros. Python puede usarse como el lenguaje que une todo esto y las posibilidades son infinitas.

Los cuadernos de Jupyter II Parte

Comenzamos con la II parte de los Cuadernos de Jupyter aplicado a las Geo-Ciencias, si necesitas mayor información no dudes en escribir.

DATOS GEOGRÁFICOS CON BASEMAP

Un tipo común de visualización en la ciencia de datos es el de los datos geográficos. La herramienta principal de Matplotlib para este tipo de visualización es el kit de herramientas Basemap, que es uno de varios kits de herramientas Matplotlib que vive bajo el espacio de nombres mpl_toolkits. Es cierto que Basemap se siente un poco torpe de usar, y muchas veces incluso las visualizaciones simples tardan mucho más tiempo en renderizarse de lo que esperas. Las soluciones más modernas como el leaflet o API de Google Maps pueden ser una mejor opción para visualizaciones de mapas más intensivas. Aún así, Basemap es una herramienta útil para los usuarios de Python en sus herramientas virtuales. En esta sección, mostraremos varios ejemplos del tipo de visualización de mapas que es posible con este kit de herramientas.

La instalación de Basemap es simple; Si está utilizando conda, puede escribir esto y el paquete se descargará:

$ conda install basemap

Agregamos solo una nueva importación a nuestro estándar repetitivo:

In [1]:

%matplotlib inline

import numpy as np

import matplotlib.pyplot as plt

from mpl_toolkits.basemap import Basemap

Una vez que haya instalado e importado el juego de herramientas Basemap, las parcelas geográficas están a solo unas líneas de distancia (los gráficos a continuación también requieren el paquete PIL en Python 2, o el paquete pillow en Python 3):

In [2]:

plt.figure(figsize=(8, 8))

m = Basemap(projection='ortho', resolution=None, lat_0=0, lon_0=-73)

m.bluemarble(scale=0.6);

Lo útil es que el globo que se muestra aquí no es una mera imagen; ¡es un eje Matplotlib completamente funcional que comprende coordenadas esféricas y que nos permite sobrerreglar fácilmente datos en el mapa! Por ejemplo, podemos usar una proyección de mapa diferente, acercarnos a Suramérica y trazar la ubicación de Quito. Usaremos una imagen etopo (que muestra las características topográficas tanto en tierra como debajo del océano) como fondo del mapa:

In [13]:

fig = plt.figure(figsize=(8, 8))

m = Basemap(height=16700000,width=12000000,

            resolution='l',area_thresh=1000.,projection='omerc',\

            lon_0=-100,lat_0=15,lon_2=-120,lat_2=65,lon_1=-50,lat_1=-55)

m.drawcoastlines()

m.fillcontinents(color='coral',lake_color='aqua')

# draw parallels and meridians.

m.drawparallels(np.arange(-80.,81.,20.))

m.drawmeridians(np.arange(-180.,181.,20.))

m.drawmapboundary(fill_color='aqua')

# Map (long, lat) to (x, y) for plotting

x, y = m(-78.503722, -0.216471)

plt.plot(x, y, 'ok', markersize=5)

plt.text(x, y, ' Quito', fontsize=12);

Esto le da una breve visión del tipo de visualizaciones geográficas que son posibles con solo unas pocas líneas de Python. Ahora analizaremos las características del mapa base con más profundidad y proporcionaremos varios ejemplos de visualización de datos de mapas. Usando estos ejemplos breves como bloques de construcción, debería poder crear casi cualquier visualización de mapa que desee.

Proyecciones de mapas

Lo primero que debe decidir al usar mapas es qué proyección usar. Probablemente esté familiarizado con el hecho de que es imposible proyectar un mapa esférico, como el de la Tierra, sobre una superficie plana sin distorsionarlo o romper su continuidad. Estas proyecciones se han desarrollado en el transcurso de la historia humana, ¡y hay muchas opciones! Dependiendo del uso previsto de la proyección del mapa, existen ciertas características del mapa (por ejemplo, dirección, área, distancia, forma u otras consideraciones) que son útiles para mantener.

El paquete de mapa base implementa varias docenas de tales proyecciones, todas referenciadas por un código de formato corto. Aquí mostraremos brevemente algunos de los más comunes.

Comenzaremos por definir una rutina de conveniencia para dibujar nuestro mapa del mundo junto con las líneas de longitud y latitud:

In [30]:

from itertools import chain

​

def draw_map(m, scale=0.2):

    # draw a shaded-relief image

    m.shadedrelief(scale=scale)

   

    # lats and longs are returned as a dictionary

    lats = m.drawparallels(np.linspace(-90, 90, 13))

    lons = m.drawmeridians(np.linspace(-180, 180, 13))

​

    # keys contain the plt.Line2D instances

    lat_lines = chain(*(tup[1][0] for tup in lats.items()))

    lon_lines = chain(*(tup[1][0] for tup in lons.items()))

    all_lines = chain(lat_lines, lon_lines)

   

    # cycle through these lines and set the desired style

    for line in all_lines:

        line.set(linestyle='-', alpha=0.3, color='w')

Proyecciones cilíndricas

Las proyecciones de mapa más simples son proyecciones cilíndricas, en las que las líneas de latitud y longitud constantes se asignan a líneas horizontales y verticales, respectivamente. Este tipo de mapeo representa bastante bien las regiones ecuatoriales, pero produce distorsiones extremas cerca de los polos. El espaciado de las líneas de latitud varía entre las diferentes proyecciones cilíndricas, que conducen a diferentes propiedades de conservación y a una distorsión diferente cerca de los polos. En la siguiente figura mostramos un ejemplo de la proyección cilíndrica equidistante, que elige una escala de latitud que preserva las distancias a lo largo de los meridianos. Otras proyecciones cilíndricas son las proyecciones de Mercator (proyección = 'merc') y área cilíndrica igual (proyección = 'cea') .:

In [31]:

fig = plt.figure(figsize=(8, 6), edgecolor='w')

m = Basemap(projection='cyl', resolution=None,

            llcrnrlat=-90, urcrnrlat=90,

            llcrnrlon=-180, urcrnrlon=180, )

draw_map(m)

Los argumentos adicionales al mapa base para esta vista especifican la latitud (lat) y la longitud (lon) de la esquina inferior izquierda (llcrnr) y la esquina superior derecha (urcrnr) para el mapa deseado, en unidades de grados.

Proyecciones pseudo-cilíndricas

Las proyecciones pseudo-cilíndricas relajan el requisito de que los meridianos (líneas de longitud constante) permanezcan verticales; esto puede dar mejores propiedades cerca de los polos de la proyección. La proyección de Mollweide (proyección = 'moll') es un ejemplo común de esto, en el que todos los meridianos son arcos elípticos. Está construido para preservar el área a lo largo del mapa: aunque hay distorsiones cerca de los polos, el área de pequeños parches refleja el área verdadera. Otras proyecciones pseudo-cilíndricas son las proyecciones sinusoidales (proyección = 'sinu') y Robinson (proyección = 'robin').

In [32]:

fig = plt.figure(figsize=(8, 6), edgecolor='w')

m = Basemap(projection='moll', resolution=None,

            lat_0=0, lon_0=0)

draw_map(m)

Los argumentos adicionales al mapa base aquí se refieren a la latitud central (lat_0) y la longitud (lon_0) para el mapa deseado.

Proyecciones de perspectiva

Las proyecciones de perspectiva se construyen usando una opción particular de punto de perspectiva, similar a la Tierra desde un punto particular en el espacio (un punto que, para algunas proyecciones, se encuentra técnicamente dentro de la Tierra). Un ejemplo común es la proyección ortográfica (proyección = 'orto'), que muestra un lado del globo visto desde un espectador a una distancia muy larga. Como tal, puede mostrar solo la mitad del globo a la vez. Otras proyecciones basadas en perspectiva incluyen la proyección gnomónica (proyección = 'gnom') y la proyección estereográfica (proyección = 'estere'). Estos suelen ser los más útiles para mostrar pequeñas porciones del mapa.

Aquí hay un ejemplo de la proyección ortográfica:

In [33]:

fig = plt.figure(figsize=(8, 8))

m = Basemap(projection='ortho', resolution=None,

            lat_0=50, lon_0=0)

draw_map(m);

Proyecciones cónicas

Una proyección Cónica proyecta el mapa en un solo cono, que luego se desenrolla. Esto puede conducir a propiedades locales muy buenas, pero las regiones alejadas del punto de enfoque del cono pueden distorsionarse mucho. Un ejemplo de esto es la proyección Cónica Conforme de Lambert (proyección = 'lcc'). Proyecta el mapa en un cono dispuesto de tal manera que dos paralelos estándar (especificados en el mapa base por lat_1 y lat_2) tienen distancias bien representadas, con una escala decreciente entre ellos y un aumento fuera de ellos. Otras proyecciones cónicas útiles son la proyección cónica equidistante (proyección = 'eqdc') y la proyección de área igual de Albers (proyección = 'aea'). Las proyecciones cónicas, como las proyecciones de perspectiva, tienden a ser buenas opciones para representar parcelas pequeñas o medianas del globo.

In [34]:

fig = plt.figure(figsize=(8, 8))

m = Basemap(projection='lcc', resolution=None,

            lon_0=0, lat_0=50, lat_1=45, lat_2=55,

            width=1.6E7, height=1.2E7)

draw_map(m)

Otras proyecciones

Si vas a hacer mucho con visualizaciones basadas en mapas, te recomiendo que leas otras proyecciones disponibles, junto con sus propiedades, ventajas y desventajas. Lo más probable es que estén disponibles en el paquete de mapa base http://matplotlib.org/basemap/users/mapsetup.html. Si profundizas lo suficiente en este tema, encontrarás una increíble subcultura de geeks de geo-viz que estarán listos para argumentar fervientemente en apoyo de su proyección favorita para cualquier aplicación dada.

Dibujando un fondo de mapa

Anteriormente vimos los métodos bluemarble () y shadedrelief () para proyectar imágenes globales en el mapa, así como los métodos drawparallels () y drawmeridians () para dibujar líneas de latitud y longitud constantes. El paquete de mapa base contiene una gama de funciones útiles para dibujar bordes de características físicas como continentes, océanos, lagos y ríos, así como límites políticos, como países y estados y condados de EE. UU. Las siguientes son algunas de las funciones de dibujo disponibles que puede explorar con las funciones de ayuda de IPython:

Límites físicos y cuerpos de agua

drawcoastlines (): dibuja líneas de costa continental drawlsmask (): dibuja una máscara entre la tierra y el mar, para utilizar con la proyección de imágenes en uno u otro drawmapboundary (): dibuja el límite del mapa, incluido el color de relleno para los océanos. drawrivers (): Dibuja ríos en el mapa fillcontinents (): llena los continentes con un color dado; opcionalmente llenar lagos con otro color Límites políticos

drawcountries (): dibujar los límites del país drawstates (): dibuja los límites del estado de EE. UU. drawcounties (): Dibuja los límites del condado de EE. UU. Características del mapa

drawgreatcircle (): Dibuja un gran círculo entre dos puntos drawparallels (): Dibuja líneas de latitud constante drawmeridians (): Dibuja líneas de longitud constante drawmapscale (): Dibuja una escala lineal en el mapa

Imágenes de globo entero

bluemarble (): Proyecte la imagen de mármol azul de la NASA en el mapa shadedrelief (): Proyecte una imagen sombreada en relieve en el mapa etopo (): Dibuje una imagen de relieve etopo en el mapa warpimage (): proyecta una imagen proporcionada por el usuario en el mapa

Para las características basadas en límites, debe establecer la resolución deseada al crear una imagen de mapa base. El argumento de resolución de la clase Basemap establece el nivel de detalle en los límites, ya sea 'c' (crudo), 'l' (bajo), 'i' (intermedio), 'h' (alto), 'f' (completo) , o Ninguno si no se usarán límites. Esta elección es importante: establecer límites de alta resolución en un mapa global, por ejemplo, puede ser muy lento.

Aquí hay un ejemplo de cómo dibujar los límites tierra / mar, y el efecto del parámetro de resolución. Crearemos un mapa de baja y alta resolución de la bella isla de Skye de Escocia. Está ubicado a 57.3 ° N, 6.2 ° W, y un mapa de 90,000 × 120,000 kilómetros:

In [35]:

fig, ax = plt.subplots(1, 2, figsize=(12, 8))

​

for i, res in enumerate(['l', 'h']):

    m = Basemap(projection='gnom', lat_0=57.3, lon_0=-6.2,

                width=90000, height=120000, resolution=res, ax=ax[i])

    m.fillcontinents(color="#FFDDCC", lake_color='#DDEEFF')

    m.drawmapboundary(fill_color="#DDEEFF")

    m.drawcoastlines()

    ax[i].set_title("resolution='{0}'".format(res));

Tenga en cuenta que las líneas costeras de baja resolución no son adecuadas para este nivel de zoom, mientras que las de alta resolución funcionan bien. Sin embargo, el nivel bajo funcionaría bien para una vista global, ¡y sería mucho más rápido que cargar los datos de borde de alta resolución para todo el mundo! Puede requerir cierta experimentación encontrar el parámetro de resolución correcto para una vista determinada: la mejor ruta es comenzar con una gráfica rápida de baja resolución y aumentar la resolución según sea necesario.

Trazar datos en mapas

Tal vez la pieza más útil del conjunto de herramientas Basemap es la capacidad de sobreponer una variedad de datos en un fondo de mapa. Para gráficos y texto simples, cualquier función plt funciona en el mapa; puede usar la instancia del mapa base para proyectar coordenadas de latitud y longitud a coordenadas (x, y) para trazar con plt.

Además de esto, hay muchas funciones específicas del mapa disponibles como métodos de la instancia del mapa base. Estos funcionan de forma muy similar a sus contrapartidas Matplotlib estándar, pero tienen un argumento booleano latlon adicional, que si se establece en True le permite pasar latitudes y longitudes sin formato al método, en lugar de las coordenadas proyectadas (x, y).

Algunos de estos métodos específicos del mapa son:

contour () / contourf (): Dibuja líneas de contorno o contornos rellenos imshow (): Dibuja una imagen pcolor () / pcolormesh (): dibuja un gráfico de pseudocolor para mallas irregulares / regulares plot (): Dibuja líneas y / o marcadores. scatter (): Dibuja puntos con marcadores. carcaj (): dibujar vectores. púas (): dibuja barbas de viento. drawgreatcircle (): Dibuja un gran círculo.

Veremos algunos ejemplos de algunos de ellos mientras continuamos. Para obtener más información sobre estas funciones, incluidos varios diagramas de ejemplo, consulte la documentación en línea del mapa base http://matplotlib.org/basemap/

Ejemplo: Ciudades de California

Recordemos que en Customizing Plot Legends, demostramos el uso del tamaño y el color en un diagrama de dispersión para transmitir información sobre la ubicación, el tamaño y la población de las ciudades de California. Aquí, crearemos este gráfico de nuevo, pero usando Basemap para poner los datos en contexto.

Comenzamos cargando los datos, como lo hicimos antes:

In [48]:

import pandas as pd

cities = pd.read_csv('data/california_cities.csv')

​

# Extract the data we're interested in

lat = cities['latd'].values

lon = cities['longd'].values

population = cities['population_total'].values

area = cities['area_total_km2'].values

A continuación, configuramos la proyección del mapa, dispersamos los datos y luego creamos una barra de colores y una leyenda:

In [49]:

# 1. Draw the map background

fig = plt.figure(figsize=(8, 8))

m = Basemap(projection='lcc', resolution='h',

            lat_0=37.5, lon_0=-119,

            width=1E6, height=1.2E6)

m.shadedrelief()

m.drawcoastlines(color='gray')

m.drawcountries(color='gray')

m.drawstates(color='gray')

​

# 2. scatter city data, with color reflecting population

# and size reflecting area

m.scatter(lon, lat, latlon=True,

          c=np.log10(population), s=area,

          cmap='Reds', alpha=0.5)

​

# 3. create colorbar and legend

plt.colorbar(label=r'$\log_{10}({\rm population})$')

plt.clim(3, 7)

​

# make legend with dummy points

for a in [100, 300, 500]:

    plt.scatter([], [], c='k', alpha=0.5, s=a,

                label=str(a) + ' km$^2$')

plt.legend(scatterpoints=1, frameon=False,

           labelspacing=1, loc='lower left');

Esto nos muestra aproximadamente donde poblaciones más grandes de personas se han asentado en California: se agrupan cerca de la costa en las áreas de Los Ángeles y San Francisco, se extienden a lo largo de las carreteras en el plano valle central y evitan casi por completo las regiones montañosas a lo largo de las fronteras el estado.

Ejemplo: datos de temperatura superficial

Como ejemplo de visualizar algunos datos geográficos más continuos, consideremos el "vórtice polar" que golpeó la mitad oriental de los Estados Unidos en enero de 2014. Una gran fuente para cualquier tipo de datos climáticos es el Instituto Goddard de Estudios Espaciales de la NASA. Aquí usaremos los datos de temperatura de GIS 250, que podemos descargar usando comandos de shell (estos comandos pueden tener que ser modificados en máquinas con Windows). Los datos utilizados aquí se descargaron el 6/12/2016 y el tamaño del archivo es de aproximadamente 9 MB:

In [ ]:

# !curl -O http://data.giss.nasa.gov/pub/gistemp/gistemp250.nc.gz

# !gunzip gistemp250.nc.gz

Los datos vienen en formato NetCDF, que la biblioteca netCDF4 puede leer en Python. Puede instalar esta biblioteca como se muestra aquí

$ conda install netcdf4

Leemos los datos de la siguiente manera :

In [50]:

from netCDF4 import Dataset

data = Dataset('data/gistemp250.nc')

El archivo contiene muchas lecturas de temperatura global en una variedad de fechas; tenemos que seleccionar el índice de la fecha que nos interesa, en este caso, 15 de enero de 2014:

In [51]:

from netCDF4 import date2index

from datetime import datetime

timeindex = date2index(datetime(2014, 1, 15),

                       data.variables['time'])

Ahora podemos cargar los datos de latitud y longitud, así como la anomalía de temperatura para este índice:

In [52]:

lat = data.variables['lat'][:]

lon = data.variables['lon'][:]

lon, lat = np.meshgrid(lon, lat)

temp_anomaly = data.variables['tempanomaly'][timeindex]

Finalmente, usaremos el método pcolormesh () para dibujar una malla de color de los datos. Veremos América del Norte y utilizaremos un mapa en relieve sombreado en el fondo. Tenga en cuenta que para estos datos elegimos específicamente un mapa de color divergente, que tiene un color neutral en cero y dos colores que contrastan en valores negativos y positivos. También dibujaremos ligeramente las costas sobre los colores para referencia:

In [53]:

fig = plt.figure(figsize=(10, 8))

m = Basemap(projection='lcc', resolution='c',

            width=8E6, height=8E6,

            lat_0=45, lon_0=-100,)

m.shadedrelief(scale=0.5)

m.pcolormesh(lon, lat, temp_anomaly,

             latlon=True, cmap='RdBu_r')

plt.clim(-8, 8)

m.drawcoastlines(color='lightgray')

​

plt.title('January 2014 Temperature Anomaly')

plt.colorbar(label='temperature anomaly (°C)');

Los datos muestran una imagen de las anomalías localizadas de temperatura extrema que ocurrieron durante ese mes. La mitad oriental de los Estados Unidos era mucho más fría de lo normal, mientras que la mitad occidental y Alaska eran mucho más cálidas. Las regiones sin temperatura registrada muestran el fondo del mapa.