Análisis de la precipitación regional con Python

A continuación se analizarán los datos de las precipitaciones promedio anuales a escala regional para el norte de Perú. El conjunto de datos cuenta con 20 estaciones de la costa, los andes y la selva de Perú con las diferencias en la elevación y las tasas de precipitación.

El objetivo principal de esta es analizar los patrones de precipitación y sus tendencias utilizando Python.

Para esta tarea utilizamos IPython Notebook que usted puede instalar en su propio ordenador o ejecutarlo en la nube desde wakari.io. Además, puede descargar los datos de la estación para este ejercicio aquí.

 

La importación del archivo

Primero importamos los archivos en el Notebook IPython que se ejecuta en su navegador favorito.

In[]:%pylab inline
%cd C:\Users\Saul\Dropbox\Curso_19_Python_en_Hidrologia_v2\1_Doc\AuxData
import numpy as np
import matplotlib.pyplot as plt
datos = np.genfromtxt('DatosEstaciones.csv', delimiter=",")
datos

Out[]:
array([[ nan,nan,nan,nan,nan],
[-76.25, -10.68,4334.,nan,1052.04],
[-75.75,-9.9 ,1860.,nan, 378.],
[-76.78,-9.55,3430.,nan, 702.],
[-77.36, -10.06,3950.,nan, 755.04],
[-75.53, -10.03, 680.,nan,2424.96],
[-75.9 , -11.48,3632.,nan, 735.96],
[-77.43,-9.73,3561.,nan, 744.96],
[-76.21, -11.63,4193.,nan, 741.],
[-75.93, -11.55,3750.,nan, 567.],
[-75.13, -10.61,1050.,nan,3096.96],
[-75.3 , -11.13, 800.,nan,2007.96],
[-77.83, -10.68,46.,nan, 9.96],
[-76.38, -11.83,2379.,nan, 237.96],
[-75.95,-9.13, 669.,nan,3203.04],
[-77.51,-9.5 ,3050.,nan, 579.],
[-74.9 , -10.3 , 301.,nan,3309.],
[-77.6 ,-9.35,2760.,nan, 737.04],
[-78.05, -10.28, 100.,nan, 6.],
[-74.93,-9.86, 250.,nan,1701.96],
[-75.46, -11.78,3322.,nan, 714.]])

La primera fila son los nombres de columna que no podían ser importados con un arreglo numpy. Nombres de columnas son: longitud, latitud, elevación, nombre y precipitación anual.

No hay necesidad de fijar la matriz para eliminar la variable no-numérica indeseada "nan", ya que vamos a hacer arreglo para cortar el script.

 

Gráfica de los valores por separado

Se hará un par de gráficos para ver la distribución espacial de la precipitación y la elevación de las estaciones.

Distribución de las precipitaciones

En la primera gráfica se muestran las estaciones como círculos escalados con su precipitación anual promedio. El color es uniforme para todas las estaciones. El script también guarda el gráfico resultante como un archivo png.

In[]:plt.scatter(datos[:,0],datos[:,1], marker="o", s=datos[:,4])
plt.savefig('puntosconmagnitud.png')
Out[]:
Sin título.png

La precipitación aumenta hacia el este de Perú donde se encuentran los grandes círculos. En los andes las estaciones tienen círculos relativamente similares, mientras que en la costa de Perú la precipitación es mínima (dos pequeños puntos al lado derecho).

Distribución de Elevación

Hacemos un gráfico similar, pero esta vez con elevación de la estación trazó como hexágonos con su color y tamaño a escala de la elevación.

Sin título.png

Podemos ver los hexágonos más grandes y rojos que corresponden a las estaciones de mayor altura en los Andes. En ambos lados de los Andes (la costa del desierto y la selva tropical) elevaciones son más bajos.

 

Gráficos 3D para mostrar las precipitaciones y la elevación

Se ha analizado la elevación y precipitación por separado, ahora se necesitan algunas herramientas como gráficos o diagramas para representarlos en el mismo gráfico. En esta parte vamos a navegar entre algunas opciones gráficas para representar la precipitación y elevación en un mismo gráfico.

Elevación y precipitación en forma de puntos

La siguiente gráfica muestra la altura como puntos rojos y las precipitaciones como puntos azules. En el fondo de la caja podemos ver una ubicación de la estación como puntos grises.

In[]:from mpl_toolkits.mplot3d import Axes3D
import numpy as np
import matplotlib.pyplot as plt
fig = plt.figure()
fig.set_size_inches(10,10)
ax = fig.gca(projection='3d')
este = datos[:,0]
norte = datos[:,1]
elev = datos[:,2]
precip = datos[:,4]
ax.scatter(este, norte, elev, marker='o', c='r',s=(datos[:,4]/20))
ax.scatter(este, norte, precip, marker='o', c='b',s=(datos[:,4]/20))
ax.scatter(este, norte, 0, marker='o', c='#4b4b4b',s=(datos[:,4]/20))
ax.set_zlim3d(0, 5000)
plt.show()

Out[]:
Sin título.png

La representación de la elevación y la precipitación es posible en este gráfico 3D. En un enfoque "orientado al usuario" necesitamos una gráfica más "amigable" y más "habladora".

Se hace un gráfico similar, pero esta vez vamos a utilizar columnas en lugar de puntos. Tendremos dos columnas al lado del otro, una roja para la elevación y otra azul para la precipitación.

In[]:from mpl_toolkits.mplot3d import Axes3D
import numpy as np
import matplotlib.pyplot as plt
fig = plt.figure()
fig.set_size_inches(10,10)
ax = fig.gca(projection='3d')
este = datos[:,0]
norte = datos[:,1]
datum = np.zeros(este.size)
elev = datos[:,2]
precip = datos[:,4]
ax.bar3d(datos[:,0]-0.1,datos[:,1]-0.1,datum, .1, .1, datos[:,2], color='r', alpha=0.5)
ax.bar3d(datos[:,0],datos[:,1],datum, .1, .1, datos[:,4], color='b',alpha=0.5)
ax.set_xlabel('Este')
ax.set_ylabel('Norte')
ax.set_zlabel('Elevacion / Precipitacion')
ax.set_zlim3d(0, 5000)
plt.show()

Out[]:

Sin título.png

Esta gráfica es mejor para representar los valores de precipitación y elevación.

 

Interpolación y trazado de cuadrícula 3D

En lugar de utilizar valores discretos, vamos a interpolar la precipitación para tener una precipitación distribuida. Scipy tiene el módulo interpolate.griddata que es adecuado para este ejercicio. Al final de la secuencia de comandos se muestra la malla de interpolación.

In[]: from scipy import interpolate
%pylab inline
#generamos parametros para la grilla
delta=0.25
este_r = np.arange(-78.5,-74.5,delta)
norte_r = np.arange(-12.0,-9.0,delta)
este_g,norte_r = np.meshgrid(este_r,norte_r)
#generamos la grilla regular
X,Y =np.meshgrid(np.linspace(-78.5,-74.5,100), np.linspace(-12,-8,100))
#realizamos la interpolación
grid_datos = interpolate.griddata((datos[1:,0],datos[1:,1]),
datos[1:,4],(X,Y)
,method='cubic')
clf()
contourf(X,Y,grid_datos,8,cmap='Blues')
colorbar()
scatter(datos[:,0],datos[:,1])
xlim(np.min(datos[1:,0]),np.max(datos[1:,0]))
ylim(np.min(datos[1:,1]),np.max(datos[1:,1]))

Out[]:
Sin título.png

Es posible mostrar la precipitación interpolada como contornos con etiquetas.

In[]:clabel(contour(X,Y,grid_datos,8,cmap='Blues'),fmt='%r')
scatter(datos[:,0],datos[:,1])
xlim(np.min(datos[1:,0]),np.max(datos[1:,0]))
ylim(np.min(datos[1:,1]),np.max(datos[1:,1]))

Out[]:
Sin título.png

Finalmente se puede trazar una malla 3D de interpolación y tiene proyección cada eje.

In[]:from mpl_toolkits.mplot3d import Axes3D
%pylab inline
fig = figure()
ax = Axes3D(fig)
cset = ax.contourf(X, Y, grid_datos, zdir='z',
 offset=-500, cmap=cm.coolwarm)
cset = ax.contourf(X, Y, grid_datos, zdir='x',
 offset=-78.5, cmap=cm.coolwarm)
ax.plot_surface(X,Y,grid_datos,rstride=8,
cstride=8,alpha=0.3,
linewidth=0.1, antialiased=True)
Out[]:
Sin título.png
Footer_Suscribirse.png

Suscríbase a nuestro boletín gratuito para recibir noticias, datos interesantes y fechas de nuestros cursos en recursos hídricos.

Posted on March 5, 2014 and filed under Hidroinformática, Hidrología.