#!/usr/bin/python3
# -*- coding: utf-8 -*-
########################################################################
#
# This file is part of python module <pyspc>.
# Copyright (C) 2013-2021 R. Marty
# (renaud.marty@developpement-durable.gouv.fr)
#
# This program is free software: you can redistribute it and/or modify
# it under the terms of the GNU General Public License as published by
# the Free Software Foundation, either version 3 of the License, or
# (at your option) any later version.
#
# This program is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
# GNU General Public License for more details.
#
# You should have received a copy of the GNU General Public License
# along with this program (see COPYING.txt).
# If not, see <http://www.gnu.org/licenses/>.
#
########################################################################
"""
Evapo-transpiration par Oudin
"""
from datetime import timedelta as td
import math
import numpy as np
import pandas as pnd
import pyspc.core.exception as _exception
[docs]
def daily2hourly(data_ej=None):
"""
Ventilation de l'ETP journalière au pas de temps horaire.
Parameters
----------
data_ej : pnd.DataFrame
Dataframe des données ETP journalières
Returns
-------
data_eh : pnd.DataFrame
Dataframe des données ETP horaires
Examples
--------
>>> import pyspc.model.oudin as _oudin
>>> data_ej = pnd.DataFrame({
... 'Kxxx_EJ': [0.9, 0.520, 0.376]} ,
... index=pnd.date_range(
... dt(2014, 11, 3),
... dt(2014, 11, 5),
... freq='D'
... )
... )
>>> data_ej
Kxxx_EJ
2014-11-03 0.900
2014-11-04 0.520
2014-11-05 0.376
>>> data_eh = _oudin.daily2hourly(data_ej)
>>> data_eh
Kxxx_EJ
2014-11-03 00:00:00 0.000000
2014-11-03 01:00:00 0.000000
2014-11-03 02:00:00 0.000000
2014-11-03 03:00:00 0.000000
2014-11-03 04:00:00 0.000000
2014-11-03 05:00:00 0.000000
2014-11-03 06:00:00 0.000000
2014-11-03 07:00:00 0.031500
2014-11-03 08:00:00 0.055800
2014-11-03 09:00:00 0.071100
2014-11-03 10:00:00 0.087300
2014-11-03 11:00:00 0.099000
2014-11-03 12:00:00 0.105300
... ...
2014-11-05 12:00:00 0.043992
2014-11-05 13:00:00 0.043992
2014-11-05 14:00:00 0.041360
2014-11-05 15:00:00 0.036472
2014-11-05 16:00:00 0.029704
2014-11-05 17:00:00 0.023312
2014-11-05 18:00:00 0.013160
2014-11-05 19:00:00 0.000000
2014-11-05 20:00:00 0.000000
2014-11-05 21:00:00 0.000000
2014-11-05 22:00:00 0.000000
2014-11-05 23:00:00 0.000000
Notes
-----
Voir le site internet du projet [1]_, et les publications de L. Oudin
[2]_ et [3]_.
.. [1] http://webgr.irstea.fr/modeles/modele-devapotranspiration/
.. [2] Oudin, L.; Hervieu, F.; Michel, C.; Perrin, C.; Andreassian, V.;
Anctil, F. & Loumagne, C. Which potential evapotranspiration input
for a lumped rainfall-runoff model?: Part 2 Towards a simple and
efficient potential evapotranspiration model for rainfall-runoff
modelling Journal of Hydrology , 2005, 303, 290 - 306
.. [3] Oudin, L. Recherche d'un modèle d'évapotranspiration potentielle
pertinent comme entrée d'un modèle pluie-débit global ENGREF
(AgroParisTech), 2004
"""
# Contrôles
_exception.check_dataframe(data_ej)
# Pré-traitement et contrôle des données journalières
data_ej = data_ej.sort_index()
ts = list(set(data_ej.index.shift(1) - data_ej.index))
_exception.raise_valueerror(
len(ts) != 1,
'Données à pas de temps multiples'
)
_exception.raise_valueerror(
ts[0] != td(days=1),
'Données non journalières'
)
# Création du dataframe final
data_eh = pnd.DataFrame(
index=pnd.date_range(
start=data_ej.index[0].replace(hour=0),
end=data_ej.index[-1].replace(hour=23),
freq='h'
),
columns=list(data_ej.columns)
)
# Répartition horaire de l'ETP journalière (de 0TU à 23TU)
ratios = np.array([
0, 0, 0, 0, 0, 0,
0, 0.035, 0.062, 0.079, 0.097, 0.11,
0.117, 0.117, 0.11, 0.097, 0.079, 0.062,
0.035, 0, 0, 0, 0, 0
])
# Remplissage
for c in data_ej.columns:
eh = []
ej = list(data_ej[c].values)
for x in ej:
eh.extend(ratios * x)
data_eh[c] = eh
return data_eh
[docs]
def daily_ETP(data_tj=None, latitude=None):
"""
Calculer l'évapo-transpiration potentielle journalière.
A partir des températures journalières et de la latitude
Parameters
----------
data_tj : pnd.DataFrame
Dataframe des données TA journalières
latitude : float
Valeur de la latitude où l'utilisateur souhaite calculer l'ETP
Returns
-------
data_ej : pnd.DataFrame
Dataframe des données ETP journalières
Examples
--------
>>> import pyspc.model.oudin as _oudin
>>> data_tj = pnd.DataFrame(
... {'Kxxx_EJ': [11.625, 9.929, 9.267]},
... index=pnd.date_range(
... dt(2014, 11, 1),
... dt(2014, 11, 3),
... freq='D'
... )
... )
>>> data_tj
Kxxx_EJ
2014-11-01 11.625
2014-11-02 9.929
2014-11-03 9.267
>>> data_ej = _oudin.daily_ETP(data_tj)
>>> data_ej
Kxxx_EJ
2014-11-01 1.10308
2014-11-02 0.97756
2014-11-03 0.92200
Notes
-----
Voir le site internet du projet [1]_, et les publications de L. Oudin
[2]_ et [3]_.
.. [1] http://webgr.irstea.fr/modeles/modele-devapotranspiration/
.. [2] Oudin, L.; Hervieu, F.; Michel, C.; Perrin, C.; Andreassian, V.;
Anctil, F. & Loumagne, C. Which potential evapotranspiration input
for a lumped rainfall-runoff model?: Part 2 Towards a simple and
efficient potential evapotranspiration model for rainfall-runoff
modelling Journal of Hydrology , 2005, 303, 290 - 306
.. [3] Oudin, L. Recherche d'un modèle d'évapotranspiration potentielle
pertinent comme entrée d'un modèle pluie-débit global ENGREF
(AgroParisTech), 2004
"""
# Contrôles
_exception.check_dataframe(data_tj)
_exception.check_numeric(latitude)
# Pré-traitement et contrôle des données journalières
data_tj = data_tj.sort_index()
ts = list(set(data_tj.index.shift(1) - data_tj.index))
_exception.raise_valueerror(
len(ts) != 1,
'Données à pas de temps multiples'
)
_exception.raise_valueerror(
ts[0] != td(days=1),
'Données non journalières'
)
# Création du dataframe final
data_ej = pnd.DataFrame(
index=data_tj.index,
columns=list(data_tj.columns)
)
# Remplissage du dataframe final
for c in data_tj.columns:
data_ej[c] = _compute_ETP(
data_ej.index.dayofyear,
data_tj[c].values,
latitude
)
return data_ej
def _compute_ETP(doy, tmp, lat):
"""
Calcul de l'ETP d'Oudin
Parameters
----------
doy : list
Jours de l'année
tmp : list
Températures
lat : float
Latitude
Returns
-------
etp : list
ETP
Notes
-----
Voir le site internet du projet [1]_, et les publications de L. Oudin
[2]_ et [3]_.
.. [1] http://webgr.irstea.fr/modeles/modele-devapotranspiration/
.. [2] Oudin, L.; Hervieu, F.; Michel, C.; Perrin, C.; Andreassian, V.;
Anctil, F. & Loumagne, C. Which potential evapotranspiration input
for a lumped rainfall-runoff model?: Part 2 Towards a simple and
efficient potential evapotranspiration model for rainfall-runoff
modelling Journal of Hydrology , 2005, 303, 290 - 306
.. [3] Oudin, L. Recherche d'un modèle d'évapotranspiration potentielle
pertinent comme entrée d'un modèle pluie-débit global ENGREF
(AgroParisTech), 2004
"""
# Contrôles
_exception.check_numeric(lat)
_exception.raise_valueerror(
len(doy) != len(tmp),
'Données incohérentes'
)
# Calcul
etp = []
for d, t in zip(doy, tmp):
# Date en jour julien
JD = d
# Application de l'algorithme d'Oudin
RD = 57.3
# Calculation of extra-atmospheric global radiation
# (Appendix C in Morton (1983), Eq. C-6 to C-11, p.60-61)
# Converts latitude in radians
FI = lat / RD
COSFI = math.cos(FI)
# AFI = abs(latitude / 42.0)
# TETA: Declination of the sun in radians
TETA = 0.4093 * math.sin(JD / 58.1 - 1.405)
COSTETA = math.cos(TETA)
COSGZ = max([0.001, math.cos(FI - TETA)])
# GZ: Noon angular zenith distance of the sun
# GZ = math.acos(COSGZ)
# COSGZ2 = COSGZ * COSGZ
# if COSGZ2 >= 1:
# SINGZ = 0.
# else:
# SINGZ = math.sqrt(1.0 - COSGZ2)
COSOM = 1.0 - COSGZ / COSFI / COSTETA
COSOM = min(max(COSOM, -1.0), 1.0)
# if COSOM < -1.0:
# COSOM = -1.0
# if COSOM > 1.0:
# COSOM = 1.0
COSOM2 = COSOM * COSOM
if COSOM2 > 1.0:
SINOM = 0.
else:
SINOM = math.sqrt(1.0 - COSOM2)
OM = math.acos(COSOM)
# PZ: Average angular zenith distance of the sun
COSPZ = COSGZ \
+ COSFI * COSTETA * (SINOM / OM - 1.)
COSPZ = max(COSPZ, 0.001)
# if COSPZ < 0.001:
# COSPZ = 0.001
# ETA: Radius vector of the sun
ETA = 1. + math.cos(JD / 58.1) / 30.
# GE: extra-atmospheric global radiation
GE = 446. * OM * COSPZ * ETA
# Daily PE by Oudin et al. (2006) formula
DPE = max(0., GE * (t + 5.) / 100. / 28.5)
etp.append(DPE)
return etp