# ---
# jupyter:
#   jupytext:
#     text_representation:
#       extension: .py
#       format_name: light
#       format_version: '1.5'
#       jupytext_version: 1.14.4
#   kernelspec:
#     display_name: Python 3
#     language: python
#     name: python3
# ---

# +
import numpy as np

import pandas as pd
from pandas.plotting import register_matplotlib_converters
register_matplotlib_converters()

from matplotlib import pyplot as plt
import seaborn as sns
sns.set()
sns.set_context('notebook')
# -

# ## 1. Hypsometric curve
# ### Import DEM

# %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)

# ### Derive the hypsometric curve

# 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

# Example: 100 m elevation

# set a threshold
Z = 100

# n. of cells below the threshold
n = np.sum(dem < Z)
n

# value in the hypsometric 
n / ncells

# Loop for all elevation thresholds

# 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 (-)');

# ### Calculate the weights

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 (-)');

# ## Precipitation vs elevation
# ### Import data

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

# ### Mean annual precipitaion at the stations
# #### Example: station 1125

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

# Loop for all the stations

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)');

# Linear regression

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)');

# Estimate elevation for each elevation band

Zbands = np.array(Zbands)
Pz = m * Zbands + n
Pz = pd.Series(data=Pz, index=Zbands)
Pz

# ## 3. Areal precipitation

Pareal = np.sum(weights * Pz)
Pareal

Pareal2 = stns.Pan.mean()
Pareal2



StarcoderdataJupyterScriptsDedupFiltered