Nota
To download the notebook, go to the following link
Daily reference evapotranspiration field
(Última actualización 1 sep 2023)
Ejemplo para calcular la evapotranspiración de referencia (ET0).
Example to calculate the reference evapotranspiration(ET0).
[ ]:
# En caso de utilizar Google Colab, descomentar las siguientes líneas
# In case of using Google Colab, uncomment the following lines
#!pip install --no-binary shapely shapely --force
#!pip install h5netcdf
#!pip install cartopy
#!pip install metpy
[ ]:
# Importamos las librerías necesarias
# We import the necessary libraries
import xarray as xr
import h5netcdf
import datetime
import numpy as np
import metpy.calc as mpcalc
import metpy.constants as constants
import matplotlib.pyplot as plt
Definimos la fecha y hora de inicialización del pronóstico:
We define the forecast initialization date:
[ ]:
init_year = 2023
init_month = 2
init_day = 1
init_hour = 0
INIT_DATE = datetime.datetime(init_year, init_month, init_day, init_hour)
Leemos los archivos para todos los plazos de pronóstico:
We read the files for all forecast lead times:
[ ]:
start_lead_time = 0
end_lead_time = 72
# Descomentar la opción elegida:
# --------
# Opción 1: Para acceder al archivo online
# Option 1: To access files online
#!pip install s3fs
#import s3fs
#fs = s3fs.S3FileSystem(anon=True)
#files = [f'smn-ar-wrf/DATA/WRF/DET/{INIT_DATE:%Y/%m/%d/%H}/WRFDETAR_01H_{INIT_DATE:%Y%m%d_%H}_{fhr:03d}.nc' for fhr in range(start_lead_time, end_lead_time)]
#ds_list = []
#for s3_file in files:
# print(s3_file)
# if fs.exists(s3_file):
# f = fs.open(s3_file)
# ds_tmp = xr.open_dataset(f, decode_coords = 'all', engine = 'h5netcdf')
# ds_list.append(ds_tmp)
# else:
# print('The file {} does not exist'.format(s3_file))
# --------
# --------
# Opción 2: Para abrir los archivos ya descargados
# Option 2: To open the already downloaded files
#files = ['WRFDETAR_01H_{:%Y%m%d_%H}_{:03d}.nc'.format(INIT_DATE, lead_time) for lead_time in range(start_lead_time, end_lead_time + 1)]
#ds_list = []
#for filename in files:
# print(filename)
# ds_tmp = xr.open_dataset(filename, decode_coords = 'all', engine = 'h5netcdf')
# ds_list.append(ds_tmp)
# Combinamos los archivos en un unico dataset
# We combine all the files in one dataset
ds = xr.combine_by_coords(ds_list, combine_attrs = 'drop_conflicts')
Definimos una función para calcular la evapotranspiracion de referencia diaria usando la ecuación FAO Penman-Monteith:
We define a function to calculate the daily reference evapotranspiration using theFAO Penman-Monteith equation:
[ ]:
def calc_ET0(ds):
ds = ds.metpy.quantify()
days = np.unique(ds['time'].astype('datetime64[D]'))
if ds['time'][0].dt.hour != 0:
days = days[1:]
ET0_list = []
for n, day in enumerate(days[:-1]):
#Calculamos valores medios diarios de las variables meteorologicas
#We calculate the daily mean meteorological variables
ds_mean = ds[['magViento10', 'PSFC', 'T2', 'HR2']].sel(time = slice(days[n], days[n+1])).mean(dim = 'time')
#Radiacion acumulada durante el dia
#Accumulated radiation during the day
lwd_day = (ds['ACLWDNB'].sel(time = days[n + 1]) - ds['ACLWDNB'].sel(time = days[n])) # Onda larga entrante / Longwave downward
lwu_day = (ds['ACLWUPB'].sel(time = days[n + 1]) - ds['ACLWUPB'].sel(time = days[n])) # Onda larga saliente / Longwave upward
swd_day = (ds['ACSWDNB'].sel(time = days[n + 1]) - ds['ACSWDNB'].sel(time = days[n])) # Onda corta entrante / Shortwave downward
swu_day = 0.23*swd_day # Onda corta saliente / Shortwave upward
#Radiacion neta [MJ/m**2]
#Net radiation [MJ/m**2]
rn_day = ((swd_day - swu_day) + (lwd_day - lwu_day)).metpy.convert_units('MJ/m**2').metpy.dequantify()
#Calculamos la velocidad del viento a 2m [m/s]
#We calculate 2m wind speed [m/s]
wind2m_day = (0.745 * ds_mean['magViento10']).metpy.dequantify()
#Calculamos la constante psicrometrica [kPa/°C]
#We calculate the psychrometric constant [kPa/°C]
gamma = ((constants.dry_air_spec_heat_press * ds_mean['PSFC'])/(constants.molecular_weight_ratio * constants.water_heat_vaporization)).metpy.convert_units('kPa/degC').metpy.dequantify()
#Calculamos la presion parcial de vapor de saturación [kPa]
#We calculate the vapour pressure of the air at saturation [kPa]
es2m_day = (mpcalc.saturation_vapor_pressure(ds_mean['T2'])).metpy.convert_units('kPa').metpy.dequantify()
ds_mean = ds_mean.metpy.dequantify()
#Calculamos la pendiente de la curva de presión de vapor:
#We calculate the slope of the vapour pressure curve
delta_num = 4098*0.6108*np.exp(17.27*ds_mean['T2']/(ds_mean['T2'] + 237.3))
delta_den = (ds_mean['T2'] + 237.3)**2
delta = delta_num/delta_den
#Calculamos la ET0 [mm/dia]
#We calculate de ET0 [mm/day]
ET0_num1 = 0.408*delta*(rn_day)
ET0_num2 = gamma * (900/(ds_mean['T2'] + 273))*wind2m_day*es2m_day*(1 - ds_mean['HR2']/100)
ET0_den = delta + gamma*(1 + 0.34*wind2m_day)
ET0_day = (ET0_num1 + ET0_num2)/ET0_den
#Donde ET0 resulta negativa, reemplazamos por cero
#Where ET0 is negative, we change it by zero
ET0_day = ET0_day.where(ET0_day >= 0, 0)
ET0_day.name = 'ET0'
ET0_day = ET0_day.expand_dims({'time':[day]})
ET0_list.append(ET0_day)
ET0_total = xr.concat(ET0_list, dim = 'time')
return ET0_total
Calculamos la ET0 diaria y graficamos el campo del primer dia:
We calculate the daily ET0 and plot the first day field:
[ ]:
ET0 = calc_ET0(ds)
print(ET0)
plt.pcolormesh(ET0['lon'], ET0['lat'], ET0.isel(time = 0))
plt.colorbar()
<xarray.DataArray 'ET0' (time: 3, y: 1249, x: 999)>
array([[[1.5088633, 1.5132143, 1.4095869, ..., 2.132155 , 2.1326497,
2.148784 ],
[1.8586245, 1.7998089, 1.7808553, ..., 2.2562513, 2.2513156,
2.039669 ],
[1.8637698, 1.8004587, 1.7870052, ..., 2.1955142, 2.1815789,
2.0429382],
...,
[4.83387 , 4.8193235, 4.8054843, ..., 5.779634 , 5.7052574,
5.652413 ],
[4.842905 , 4.840371 , 4.833255 , ..., 5.796604 , 5.7327604,
5.6763706],
[4.84976 , 4.8477564, 4.8447657, ..., 5.7560043, 5.7638516,
5.7496104]],
[[1.6991397, 1.6982177, 1.6491075, ..., 1.4275619, 1.4659595,
1.4102427],
[1.7688186, 1.7902625, 1.7953228, ..., 1.7655959, 1.7632679,
1.3559755],
[1.7708207, 1.792065 , 1.8047507, ..., 1.7875934, 1.7759477,
1.4457896],
...
[4.264168 , 4.466234 , 4.4394145, ..., 5.4323554, 5.3156705,
5.3448086],
[4.263845 , 4.497948 , 4.515789 , ..., 5.4381967, 5.327081 ,
5.328658 ],
[4.401255 , 4.380332 , 4.3858075, ..., 5.2771287, 5.335723 ,
5.370829 ]],
[[1.8218285, 1.8230447, 1.8138897, ..., 2.2263365, 2.199836 ,
2.1875033],
[1.8340664, 1.8581399, 1.8674679, ..., 1.7753291, 1.754859 ,
0.7104781],
[1.8350779, 1.8582087, 1.868867 , ..., 1.5891689, 1.5744501,
0.5838893],
...,
[4.4109006, 4.706371 , 4.688888 , ..., 5.3870797, 5.261983 ,
5.2725234],
[4.5294986, 4.7443376, 4.7301965, ..., 5.386678 , 5.266292 ,
5.2651258],
[4.823922 , 4.8188915, 4.818117 , ..., 5.27438 , 5.3473697,
5.372849 ]]], dtype=float32)
Coordinates:
* time (time) datetime64[ns] 2023-02-01 2023-02-02 2023-02-03
lon (y, x) float32 -94.33 -94.28 -94.22 ... -48.0 -47.97
lat (y, x) float32 -54.39 -54.4 -54.41 ... -11.65 -11.65
* x (x) float32 -1.996e+06 -1.992e+06 ... 1.992e+06 1.996e+06
* y (y) float32 -2.496e+06 -2.492e+06 ... 2.492e+06 2.496e+06
Lambert_Conformal <U1 ''
<ipython-input-14-3b377cddb4be>:4: UserWarning: The input coordinates to pcolormesh are interpreted as cell centers, but are not monotonically increasing or decreasing. This may lead to incorrectly calculated cell edges, in which case, please supply explicit cell edges to pcolormesh.
plt.pcolormesh(ET0['lon'], ET0['lat'], ET0.isel(time = 0))
<matplotlib.colorbar.Colorbar at 0x7a0136f85720>
