<jupyter_start><jupyter_text>1. Hypsometric curve Import DEM<jupyter_code>%run functions_precipitation.ipynb 
dem, attributes = read_ascii('../data/dem_pas.asc')
attributes
dem.shape
dem
# plot the DEM
im = plt.imshow(dem, cmap='pink')
cb = plt.colorbar(im)
cb.set_label('elevation (masl)')
plt.axis('off');
np.nanmin(dem), np.nanmax(dem)<jupyter_output><empty_output><jupyter_text>Derive the hypsometric curve<jupyter_code># define elevation thresholds 
Zs = np.arange(start=0, stop=1701, step=100)
Zs
# n. of cells in the catchment
ncells = np.sum(~np.isnan(dem))
ncells<jupyter_output><empty_output><jupyter_text>Example: 100 m elevation<jupyter_code># set a threshold
Z = 100
# n. of cells below the threshold
n = np.sum(dem < Z)
n
# value in the hypsometric 
n / ncells<jupyter_output><empty_output><jupyter_text>Loop for all elevation thresholds<jupyter_code># table where to save the data
hypso = pd.Series(index=Zs)
hypso.head()
# compute the hypsometric curve
for Z in Zs:
    hypso[Z] = np.sum(dem < Z) / ncells
hypso.head()
# plot the hypsometric curve
plt.plot(Zs, hypso)
plt.title('Hypsometric curve')
plt.xlabel('elevation (masl)')
plt.ylabel('area (-)');<jupyter_output><empty_output><jupyter_text>Calculate the weights<jupyter_code>Zbands = []
weights = []
for Zi, Zj in zip(Zs[:-1], Zs[1:]):
    Zbands.append(np.mean([Zi, Zj]))
    weights.append(hypso[Zj] - hypso[Zi])
weights = pd.Series(data=weights, index=Zbands)
weights
plt.plot(weights)
plt.xlabel('elevation (masl)')
plt.ylabel('weight (-)');<jupyter_output><empty_output><jupyter_text>Precipitation vs elevation Import data<jupyter_code>pcp_d = pd.read_csv('../data/daily_precipitation_Pas.csv', parse_dates=True, index_col=0)
pcp_d.tail()
stns = pd.read_csv('../data/stations_pas.csv', index_col='CODE')
stns<jupyter_output><empty_output><jupyter_text>Mean annual precipitaion at the stations Example: station 1125<jupyter_code>stn = '1115'
plt.figure(figsize=(15, 3.5))
plt.plot(pcp_d.index, pcp_d[stn], linewidth=.1)
plt.ylabel('P (mm)');
# annual series
pcp_an = pcp_d[stn].groupby(by=pcp_d[stn].index.year).mean()
plt.plot(pcp_an);
daysYear = pcp_d[stn].groupby(by=pcp_d[stn].index.year).count()
plt.plot(daysYear);
thr = 330
pcp_an.loc[daysYear > 330].mean() * 365.25<jupyter_output><empty_output><jupyter_text>Loop for all the stations<jupyter_code>pcp_an = pcp_d.groupby(by=pcp_d.index.year).mean()
daysYear = pcp_d.groupby(by=pcp_d.index.year).count()
# series where to save the results
Pan = pd.Series(index=stns.index)
Pan.head()
for stn in stns.index:
    Pan[stn] = pcp_an.loc[daysYear[stn] > 330, stn].mean() * 365.25
Pan
stns['Pan'] = Pan
stns
plt.scatter(stns.Z, stns.Pan)
plt.ylim(0, 2500)
plt.xlabel('elevation (masl)')
plt.ylabel('Pannual (mm)');<jupyter_output><empty_output><jupyter_text>Linear regression<jupyter_code>from scipy.stats import linregress
m, n, *perf = linregress(stns.Z, stns.Pan)
m, n
perf
plt.scatter(stns.Z, stns.Pan)
x = np.array([0, 2000])
y = m * x + n
plt.plot(x, y, '--k')
plt.ylim(0, 2500)
plt.xlabel('elevation (masl)')
plt.ylabel('Pannual (mm)');<jupyter_output><empty_output><jupyter_text>Estimate elevation for each elevation band<jupyter_code>Zbands = np.array(Zbands)
Pz = m * Zbands + n
Pz = pd.Series(data=Pz, index=Zbands)
Pz<jupyter_output><empty_output><jupyter_text>3. Areal precipitation<jupyter_code>Pareal = np.sum(weights * Pz)
Pareal
Pareal2 = stns.Pan.mean()
Pareal2<jupyter_output><empty_output>
StarcoderdataJupyterStructuredCleanDedup