Regional precipitation analysis with Python

We want to analyze average annual precipitation data on a regional scale for the North of Peru. The dataset has 20 stations from the Peruvian desert coast, Andes and rainforest with differences in elevation and precipitation rates.

The main objective is to analyze precipitation patterns and trends by plotting them in different ways using Python.

For this task we use IPython Notebook, which you can install on your own computer or run it on the cloud from wakari.io. You can download the station data for this exercise here.

Importing the file

First, we import the files on the IPython Notebook that runs on your favorite browser. Variables are written in Spanish to differenciate them from the rest of the python commands.

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.]])

The first row are column names that could not be imported with a numpy array. Column names are: Longitude, Latitude, Elevation, Name, Annual Precipitation.

There is no need to fix the array to delete the unwated not-a-number “nan” since we will do array slicing on our script.

Plotting values separately

We will do a couple of plots to see the spatial distribution of station precipitation and station elevation.

Precipitation distribution

The first plot shows stations as circles scaled with their average annual precipitation. Color is uniform for all stations. The script also saves the resulting plot as a .png file.

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

Precipitation increases towards East of Peru where the big circles are located. In the Andes, stations have relative similar circles, while in the Peruvian coast, precipitation is minimal (two small dots on the right side).

Elevation distribution

We do a similar plot, but this time with station elevation plotted as hexagons with a colour and size scaled to the elevation.

We can see the biggest and red hexagons corresponding to stations high in the Andes. On both sides of the Andes (desert coast and rainforest) elevations are lower.

3D plotting to show precipitation and elevation

We have analyzed elevation and precipitation separately. We need some graphic tools or plots to have them on the same graphic. In this part, we will surf among some plot options to analyze both precipitation and elevetation on a same graphic.

Elevation and precipitation as dots

The next plot shows elevation as red dots and precipitation as blue dots. On the box bottom we can see a station location as grey dots.

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[]:

Representation of elevation and precipitation is possible on this 3D plot. In a “user oriented” approach we need a more “friendly” and “talkative” plot.

We will do a similar plot, but this time we will use columns instead of dots. We will have two columns next to each other, one red for elevation and another blue for precipitation.

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[]:

Interpolation and 3D grid plotting

Instead of using discrete values, we will interpolate precipitation to have a distributed precipitation. Scipy has the interpolate.griddata module that is suitable for this exercise. At the end of the script, we show the interpolated grid.

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[]:

We can show the interpolated precipitation as contours with labels.

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[]:

Finally, we can 3D plot the interpolated grid and have a fancy projection on each axis.

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[]: