⚠️ This is the first operational version of the handbook, but it is still a work in progress and will be heavily updated during 2024! ⚠️

Agricultural drought Hazard Assessment#

Click Binder to launch the workflow on Binder!

In this workflow we will assess the impact of water deficit on crop yields.

The aim is to estimate the potential loss in yield for a given crop in the absence of an artificial irrigation system compensating for precipitation scarcity. This is particularly relevant for semi-arid regions which are increasingly prone to prolonged drought periods making artificial irrigation unfeasible, as well as historically wet regions that have not yet implemented artificial irrigation at large-scale but might experience a significant decline in precipitation rates with future climate change. For more details, please refer to the description document.

Hazard Asssessment methodology#

The starting point of the hazard assessment is the calculation of the soil standard evapotranspiration (ET0) using the Penman-Monteith equation described by FAO (Allen et al., 1998). The Penman-Monteith equation allows to calculate the reference evapotranspiration demand (i.e., how much water does a standard plant requires to live) combining thermal radiation and wind data. ET0 does not hold any plant-specific information, but represents the climate-driven evapotranspiration demand for a reference crop.

The standard evapotranspiration is then combined with a series of crop-specific information to derive the crop standard evapotranspiration (ETc). The crop-specific parameters used to calculate ETc vary depending on the thermal climate zone considered. Thus, ETc represents the maximum evapotranspiration potential of each crop for a given thermal climate zone. From ETc, it is possible to derive the actual crop evapotranspiration (ETa) using precipitation data to estimate the local water availability for the plant. Hence, ETa represents a fraction of ETc. Using the equations from the FAO I&D 33 paper (Doorenbos et al., 1979), is then possible to relate the ratio between the rainfed and the maximum evapotranspiration potential to the crop yield loss (%) in rainfed-only conditions.

The hazard assessment can currently be performed for the 14 crops parameterised in the crop table, but this can be modfied to add any crop the user might be interested in. The assessment is performed using climate data averaged inter-annually to show the impact of precipitation scarcity on yield for an average growing season in the selected period.


One limitation of this approach to calculate yield losses derive from the paramterization of the crop-specific indicators. Indeed, the thermal climate zone division is rather coarse and does not allow to accurately capture all sub-regional climate characteristics. As a consequence, the definition of the season start and end days, as well as the length of the growing period might not accurately represent the real conditions in all studied regions. However, the user can overcome these limitations by manually edit the crop table csv file and insert the parameters that best represent the local conditions they want to assess.

Another limitation derives from the use of the Doorenbos et al., (1979) equation to relate water scarcity to yield losses. On one hand, this approach has been widely used in the literature as it is extremely straight-forward but still offers a good estimate of yield reductions. On the other, more sophisticated approaches have now been developed to get more accurate quantification of aridity-related yield losses taking into account additional processes and parameters (Steduto et al., 2012). For the purpose of this workflow, the choice of the Doorenbos method seemed the most approapriate due to its long-standing solidity and computational efficiency. The user is invited to cross-check the results with those produced by more sophisticated softwares like CropWat to assess the robustness of the workflow.

Load libraries#

In this notebook we will use the following Python libraries:

import pandas as pd
import numpy as np
import rasterio
import geopandas as gpd
import pooch
import cdsapi
import os
import xarray as xr
import json
import urllib
import pyproj
import cartopy.crs as ccrs
import zipfile
import matplotlib.pyplot as plt
from mpl_toolkits.basemap import Basemap
import re

Create the directory structure#

First, we need to set up the directory structure to make the workflow work. The next cell will create the directory called ‘agriculture_workflow’ in the same directory where this notebook is saved. A directory for data and one for results will also be created inside the main workflow directory to store the downloaded data and the final plots.

workflow_dir = 'agriculture_workflow'

# Define directories for data and results within the previously defined workflow directory
data_dir = os.path.join(workflow_dir, 'data')
results_dir = os.path.join(workflow_dir, 'results')

# Check if the workflow directory exists, if not, create it along with subdirectories for data and results
if not os.path.exists(workflow_dir):

Define the studied area#

The cells below allow to download the boundaries of any NUTS2 region in the EU as a GeoJson file given the region code (in this case ES51 for Catalunya). You can look up the NUTS2 code for all EU regions here by simply searching the document for the region name.

The coordinates of the selected regions are extracted and saved in an array. Finally, the geometry of the GeoJson file is saved as a shapefile to be used in the plotting phase.

region = ['ES51'] #Replace the code in [''] with that of your region 
#auxiliary function to load region GeoJson file.
def load_nuts_json(json_path):
    # dependencies: json, urllib, geopandas, 
    while True:
        uh = urllib.request.urlopen(json_path)
        data = uh.read()
    gdf = gpd.GeoDataFrame.from_features(json.loads(data)["features"])
    gdf['Location'] = gdf['CNTR_CODE'] 
    gdf = gdf.set_index('Location')
    #gdf.to_crs(pyproj.CRS.from_epsg(4326), inplace=True)
    return gdf

# load nuts2 spatial data
json_nuts_path = 'https://gisco-services.ec.europa.eu/distribution/v2/nuts/geojson/NUTS_RG_10M_2021_4326_LEVL_2.geojson'
nuts = load_nuts_json(json_nuts_path)
nuts = nuts.loc[nuts['NUTS_ID'].isin(region)]

#extract coordinates

#save geometry as shapefile

The code below creates the study area bounding box using the coordinates of the region GeoJson file. The coordinates are then reprojected from the source projection to the climate data projection.

In some cases, it might be needed to expand the selected area through the ‘scale’ parameter to avoid the corners of the region being left out from the data extraction. The units of the ‘scale’ parameter are degrees, so setting scale=1 will increase the extraction area by approximately 100 km. A scale of 0-0.5 should be sufficient to fully cover most regions.


The larger the scale parameter, the larger the extracted area, the longer the workflow will run for. Thus, the user is invited to have a first run of the workflow with scale=0.5, then increase it if not satisfied with the data coverage of the final map.


#defining region bounding box

#setting up the projection transformation tool
crs = ccrs.RotatedPole(pole_latitude=39.25, pole_longitude=-162)
transformer = pyproj.Transformer.from_crs('epsg:4326',crs)

# New bbox coordinates matching EURO-CORDEX projection.
RLON_MIN, RLAT_MIN = transformer.transform(bbox[1], bbox[0])
RLON_MAX, RLAT_MAX = transformer.transform(bbox[3], bbox[2])

Download the climate datasets#

Climate data is sourced from the Copenicus Climate Data Store. For this analysis, EURO-CORDEX climate projections for precipitation flux, maximum temperature, minimum temperature, relative humidity, wind speed and solar downward radiation at 12 km spatial resolution and daily temporal resolution have been employed. These projections are readily accessible to the public through the Climate Data Store (CDS) portal. The EURO-CORDEX data can be downloaded through the CDS API reported below.

For this example, we will guide you through downloading 5-year timeframes for each needed variable from the EURO-CORDEX repository. In this example the selected timeframe is 2046-2050 and the emission scenario is RCP2.6. We will also select a specific combination of General Circulation Model (GCM) and Regional Climate model (RCM). To download the data through the API below you will need to register an account on the Copernicus CDS. A complete guide on how to set up the API can be found here. Once your API is set-up, insert your own ‘KEY’ below and run the cell to download the data.

In the cell below, specify the starting and end year, as well as the RCP scenario of the period you are interested in. Climate data are stored in the CDS as separate 5 years simulations, so start and end years need to be specified for all of them. For simplicity, we suggest to run simulations for 2026-2046 (near future), 2036-2065 (mid-century), or 2071-2090 (end-of-century). Below, an example for the mid-century scenario is displayed.


The cell below creates a list of start and end years based on your selection. These lists are needed to download the data since daily climate data from the CDS are stored as separate 5 years simulations.

list_start=[str(x) for x in arr_start]
list_end=[str(x) for x in arr_end]


The download could take a few minutes (depending on your internet speed) due to the large dataset to download and the availability of the Copernicus repository. Once it starts, the status bar below the cell will keep you updated on the status of the download.


Only some GCM-RCM combinations will allow you to download all the variables needed in this workflow. Extracting all the variables from the same GCM-RCM combination is important for data consistency. In this example we used mpi_m_mpi_esm_lr as the GCM and smhi_rca4 as the RCM. The table below shows the 5 GCM-RCM combinations that allow to download all needed variables for all EURO-CORDEX emission scenarios (RCP2.6, RCP4.5 and RCP8.5). Although all climate projections are equally plausible, it is recommended to re-run the workflow using data from different models to test the results uncertainty.

Global Climate Model

Regional Climate Model











#remove previously downloaded files
if os.path.exists(f'{data_dir}/EUR-11*{rcp}*day*.nc'):

#start new download
zip_path_cordex = os.path.join(data_dir, 'cordex_data.zip')

KEY= None #INSERT your own

c = cdsapi.Client(url=URL,key=KEY)

        'domain': 'europe',
        'experiment': rcp_download,
        'horizontal_resolution': '0_11_degree_x_0_11_degree',
        'temporal_resolution': 'daily_mean',
        'variable': [
            '10m_wind_speed', '2m_relative_humidity', 'maximum_2m_temperature_in_the_last_24_hours',
            'mean_precipitation_flux', 'minimum_2m_temperature_in_the_last_24_hours', 'surface_solar_radiation_downwards'
        'gcm_model': 'mpi_m_mpi_esm_lr',
        'rcm_model': 'smhi_rca4',
        'ensemble_member': 'r1i1p1',
        'start_year': list_start,
        'end_year': list_end,
        'format': 'zip',

with zipfile.ZipFile(zip_path_cordex, 'r') as zObject:


You can download data for different emission scenarios, as well as start and end years by modifying the experiment, start_year and end_year parameters in the code above.

Explore the data#

The downloaded files from CDS have a filename structure that describes the exact dataset contained. The general structure for daily files is the following:


Load one of the files and explore the content and structure of the dataset. Notice the dimensions, coordinates, data variables, indexes and attributes as well as the description of the spatial reference system in rotated_pole().

# Open the netCDF file as a dataset using xarray
dataset_precipitation_example = xr.open_mfdataset(
# Display said dataset
<xarray.Dataset> Size: 1GB
Dimensions:       (rlat: 412, rlon: 424, time: 1826, bnds: 2)
    lat           (rlat, rlon) float64 1MB dask.array<chunksize=(412, 424), meta=np.ndarray>
    lon           (rlat, rlon) float64 1MB dask.array<chunksize=(412, 424), meta=np.ndarray>
  * rlat          (rlat) float64 3kB -23.38 -23.27 -23.16 ... 21.62 21.73 21.84
  * rlon          (rlon) float64 3kB -28.38 -28.27 -28.16 ... 17.94 18.05 18.16
    rotated_pole  |S1 1B ...
  * time          (time) datetime64[ns] 15kB 2046-01-01T12:00:00 ... 2050-12-...
    time_bnds     (time, bnds) datetime64[ns] 29kB dask.array<chunksize=(1, 2), meta=np.ndarray>
Dimensions without coordinates: bnds
Data variables:
    pr            (time, rlat, rlon) float32 1GB dask.array<chunksize=(1, 412, 424), meta=np.ndarray>
Attributes: (12/23)
    Conventions:                    CF-1.4
    contact:                        rossby.cordex@smhi.se
    creation_date:                  2017-02-17-T01:34:52Z
    experiment:                     RCP26
    experiment_id:                  rcp26
    driving_experiment:             MPI-M-MPI-ESM-LR, rcp26, r1i1p1
    ...                             ...
    rossby_run_id:                  201622
    rossby_grib_path:               /home/rossby/prod/201622/raw/
    rossby_comment:                 201622: CORDEX Europe 0.11 deg | RCA4 v1a...
    rcm_version_id:                 v1a
    c3s_disclaimer:                 This data has been produced in the contex...
    tracking_id:                    hdl:21.14103/3dcc2803-8bd0-4a29-b9e6-661c...

Extract the climate data for the studied region#

In the cell below the climate datasets are cut for the studied region. Data is also averaged so that the final datasets are composed of 365 values representing the inter-annual mean of the studied period. Thus, the assessment refers to a ‘typical’ annual growing season in the selected timeframe.


If you run the workflow using a different future period, remember to change the start and end years (2046 and 2050) in the file names below.

# Auxiliary function to slice each dataset to a particular region with rotated coordinates.
def cut_to_region(ds):
    ds = ds.sel(rlat = slice(RLAT_MIN, RLAT_MAX), rlon = slice(RLON_MIN, RLON_MAX))
    return ds

#Extraction of each variable and averaging across years

ds_tasmax = xr.open_mfdataset(f'{data_dir}/tasmax_EUR-11*{rcp}*day*.nc', decode_coords='all',

ds_tasmin = xr.open_mfdataset(f'{data_dir}/tasmin_EUR-11*{rcp}*day*.nc', decode_coords='all',


ds_precipitation=xr.open_mfdataset(f'{data_dir}/pr_EUR-11*{rcp}*day*.nc', decode_coords='all',
precipitation=ds_precipitation['pr'].convert_calendar('noleap').groupby('time.dayofyear').mean().to_numpy()*86400 #conversion from kg/m2/s to mm/day

#Wind speed
ds_ws=xr.open_mfdataset(f'{data_dir}/sfcWind_EUR-11*{rcp}*day*.nc', decode_coords='all',
ws=ds_ws['sfcWind'].convert_calendar('noleap').groupby('time.dayofyear').mean().to_numpy()*(4.87/np.log(67.8*10**-5.42)) #conversion from 10m wind speed to 2m wind speed

#Relative humidity
ds_hurs = xr.open_mfdataset(f'{data_dir}/hurs_EUR-11*{rcp}*day*.nc',decode_coords='all',

#Solar downward surface radiation
rad=ds_rsds['rsds'].convert_calendar('noleap').groupby('time.dayofyear').mean().to_numpy()*0.0864 #conversion from Watts to MegaJoules

Download the Available Water Capacity, Elevation and Thermal Climate Zones datasets#

To calculate the crop evapotranspiration potential (ETc) we need information about the local Available Water Capacity, elevation and thermal climate zone. The cells below show how to download these datasets and extract information for the studied region.

1. Available Water Capacity#

The available water capacity (AWC) represents the maximum amount of water the soil can store within a certain depth. In this workflow we are using data for soil depths of 0-200 cm. The data are sourced from Hengl and Gupta (2019). The cell below downloads the data as a .tif file (~5 GB size) from the zenodo repository and saves it in the data directory.

url = 'doi:10.5281/zenodo.2629149/sol_available.water.capacity_usda.mm_m_250m_0..200cm_1950..2017_v0.1.tif'
filename = 'sol_available.water.capacity_usda.mm_m_250m_0..200cm_1950..2017_v0.1.tif'


2. Elevation#

Elevation data are sourced from the USGS GDTEM 2010 digital elevation model (Danielson & Gesch, 2011). This dataset has relatively coarse resolution compared to more recent DEM, but it is adequate for the resolution of this assessment. The cell below downloads the DEM at 30 arc-second resolution as a .zip folder, then extracts it in the data directory.

url = 'https://edcintl.cr.usgs.gov/downloads/sciweb1/shared/topo/downloads/GMTED/Grid_ZipFiles/mn30_grd.zip'
filename = 'mn30_grd.zip'
with zipfile.ZipFile(f'{data_dir}/mn30_grd.zip', 'r') as zObject:

3. Thermal Climate Zones#

The thermal climate zones represent a classification of global climates in 8 categories characterised by different annual temperature ranges and rain seasons (Van Velthuizen et al., 2007). This classification is helpful for agricultural modelling as crops have different growing calendars and evapotranspirative responses depending on the temperature/precipitation regime. The cell below downloads the thermal climate zones raster from the FAO repository as a .zip folder and extracts it in the data directory.

url = 'https://storage.googleapis.com/fao-maps-catalog-data/uuid/68790fd0-690c-11db-a5a5-000d939bc5d8/resources/Map4_1.zip'
filename = 'Map4_1.zip'

with zipfile.ZipFile(f'{data_dir}/Map4_1.zip', 'r') as zObject:

Extract the AWC, elevation and thermal climate zones data for the studied region#

#Select the lon and lat coordinates from one of the climate datasets extracted before and create an array of points

#empty dataframe to store extracted data temporarily

#Sample each raster file at the coordinates found above and save the data

raster_path = f'{data_dir}/sol_available.water.capacity_usda.mm_m_250m_0..200cm_1950..2017_v0.1.tif'
with rasterio.open(raster_path) as src:
    dff['AWC'] = [x[0] for x in src.sample(coords_11)]  

with rasterio.open(raster_path) as src:
    dff['elev'] = [x[0] for x in src.sample(coords_11)]    

#Climate zones
raster_path = f'{data_dir}/therm_clim/w001001.adf'
with rasterio.open(raster_path) as src:
    dff['clim'] = [x[0] for x in src.sample(coords_11)] 
climate_zone[climate_zone<0]=dff['clim'].mode() #where climate zone gives error, assign the most commmon value.

Explore the data#

Now that we have extracted the data, you can simply call the AWC, elev and climate_zone variables to visualise them. All variables are stored as a numpy array and have the same shape (i.e., lon x lat dimension) of the climate datasets.
To see the shape of the extracted files:

(24, 34)

This tells us that in this case, the extracted area for Catalunya is formed by a lon-lat grid of 24x34 12.5 km pixels.

Import crop-specific information#

Crop specific paramters needed for the hazard assessment are imported from the crop table. The cell below reads the crop table file and prints it.


A complete description of the crop table is available here.

FAO_Code Crop Clim Kc_in Kc_mid Kc_end lgp_f1 lgp_f2 lgp_f3 lgp_f4 Season start Season End RD1 RD2 DF Type Ky
0 111 wheat 2 0.7 1.15 0.30 0.125000 0.208333 0.416667 0.250000 135.0 255 0.2 1.25 0.55 1 1.0
1 111 wheat 3 0.7 1.15 0.30 0.111111 0.333333 0.388889 0.166667 319.0 134 0.2 1.25 0.55 1 1.0
2 111 wheat 4 0.7 1.15 0.30 0.477612 0.223881 0.223881 0.074627 289.0 259 0.2 1.25 0.55 1 1.0
3 111 wheat 5 0.7 1.15 0.30 0.111111 0.333333 0.388889 0.166667 336.0 151 0.2 1.25 0.55 1 1.0
4 111 wheat 7 0.4 1.15 0.30 0.477612 0.223881 0.223881 0.074627 320.0 290 0.2 1.25 0.55 1 1.0
... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ...
79 122 sweet_potato 3 0.5 1.15 0.65 0.133333 0.200000 0.400000 0.266667 91.0 241 0.2 1.25 0.65 1 1.0
80 122 sweet_potato 4 0.5 1.15 0.65 0.133333 0.200000 0.400000 0.266667 121.0 271 0.2 1.25 0.65 1 1.0
81 122 sweet_potato 5 0.5 1.15 0.65 0.133333 0.200000 0.400000 0.266667 121.0 271 0.2 1.25 0.65 1 1.0
82 122 sweet_potato 7 0.5 1.15 0.65 0.133333 0.200000 0.400000 0.266667 121.0 271 0.2 1.25 0.65 1 1.0
83 122 sweet_potato 8 0.5 1.15 0.65 0.133333 0.200000 0.400000 0.266667 121.0 271 0.2 1.25 0.65 1 1.0

84 rows × 17 columns

The resulting table is a DataFrame that you can explore using the column keys. For instance, to check information for all available crops in the temperate oceanic climate zone you can use:

crop_table[(crop_table["Clim"] == 4)]
FAO_Code Crop Clim Kc_in Kc_mid Kc_end lgp_f1 lgp_f2 lgp_f3 lgp_f4 Season start Season End RD1 RD2 DF Type Ky
2 111 wheat 4 0.70 1.15 0.30 0.477612 0.223881 0.223881 0.074627 289.0 259 0.2 1.25 0.55 1 1.00
8 113 maize 4 0.30 1.20 0.50 0.200000 0.200000 0.266667 0.333333 106.0 256 0.2 1.00 0.55 1 1.50
14 112 rice 4 1.05 1.20 0.90 0.200000 0.200000 0.400000 0.200000 122.0 302 0.2 0.50 0.20 1 1.35
20 114 barley 4 0.30 1.15 0.25 0.477612 0.223881 0.223881 0.074627 289.0 259 0.2 1.00 0.55 1 1.00
26 152 sorghum 4 0.30 1.00 0.55 0.160000 0.280000 0.320000 0.240000 92.0 217 0.2 1.00 0.55 1 0.90
32 202 soybean 4 0.40 1.15 0.50 0.133333 0.166667 0.500000 0.200000 122.0 272 0.2 0.60 0.50 1 0.85
38 204 sunflower 4 0.35 1.15 0.35 0.192308 0.269231 0.346154 0.192308 122.0 252 0.2 0.80 0.45 1 0.95
44 121 potato 4 0.50 1.15 0.75 0.206897 0.241379 0.344828 0.206897 92.0 237 0.2 0.40 0.35 1 1.10
50 157 sugar_beet 4 0.35 1.20 0.70 0.277778 0.222222 0.277778 0.222222 136.0 331 0.2 1.00 0.55 1 0.90
56 176 beans 4 0.40 1.15 0.35 0.333333 0.333333 0.400000 0.266667 166.0 241 0.2 0.80 0.45 1 1.15
62 270 rapeseed 4 0.35 1.15 0.35 0.200000 0.400000 0.200000 0.200000 122.0 302 0.2 1.00 0.60 1 1.00
68 191 chickpea 4 1.00 0.35 0.40 0.318182 0.227273 0.272727 0.181818 105.0 215 0.2 0.80 0.50 1 0.60
74 201 lentil 4 0.40 1.10 0.30 0.133333 0.200000 0.400000 0.266667 121.0 271 0.2 0.70 0.50 1 1.00
80 122 sweet_potato 4 0.50 1.15 0.65 0.133333 0.200000 0.400000 0.266667 121.0 271 0.2 1.25 0.65 1 1.00

Hazard calculation#

Now that we have imported all the data we needed, we can start with the hazard quantification. First, you will have to select the crops you want to assess. In this example we are using wheat and maize.

#Select crops to study
crop_list = ['maize','wheat']


Write the crop names inside the crop_list exactly how they are written in the crop table.


You can select as many crops as you want, but the more the longer the workflow will take to complete.

The cell below creates empty arrays to be filled with the outputs of the calculations. Some of them are only required for debugging and validation purposes, but advanced users might want to explore them.

The fundamental parameters are:

  • irrigation requirement (irr_req), representing the amount of water needed to reach the maximum evapotranspiration potential

  • actual evapotranspiration in rainfed conditions (ETa_rf)

  • actual evapotranspiration in irrigated conditions (ETa_irr)

  • percentage yield reduction (yield_loss_perc)

The validation pameters are:

  • the daily crop coefficient (kc_daily), used to relate the crop growing stage to the water demand

  • the daily water stress coefficient (ks_daily), expressing the percentage reduction in evapotranspiration due to water shortage

  • the daily water deficit due to precipitation shortage (deficit_meteo_daily)

  • the daily rooting depth of the studied crops (daily_rooting_depth_all)

  • the effective precipitation (prec_eff_daily), representing the amount of rainfall reaching the crop roots

  • the daily standard evapotranspiration potential (ET0 daily).

# Create empty arrays for different parameters

Irr_req = np.zeros([len(Tmax[0,:,0]),len(Tmax[0,0,:]),len(crop_list),365])
ETa_rf = np.zeros([len(Tmax[0,:,0]),len(Tmax[0,0,:]),len(crop_list),365])
ETa_irr= np.zeros([len(Tmax[0,:,0]),len(Tmax[0,0,:]),len(crop_list),365])
yield_loss_perc = np.zeros([len(Tmax[0,:,0]),len(Tmax[0,0,:]),len(crop_list)])

kc_daily = np.zeros([len(Tmax[0,:,0]),len(Tmax[0,0,:]),len(crop_list),365])
ks_daily = np.zeros([len(Tmax[0,:,0]),len(Tmax[0,0,:]),len(crop_list),365])
daily_rooting_depth_all= np.zeros([len(Tmax[0,:,0]),len(Tmax[0,0,:]),len(crop_list),365])
prec_eff_daily= np.zeros([len(Tmax[0,:,0]),len(Tmax[0,0,:]),len(crop_list),365])

Reference soil evapotranspiration potential (ET0)#

The cell below defines an auxiliary function to compute the daily reference evaportanspiration potential (ET0) in every grid-cell of the selected region following the FAO56 methodology (Allen et al., 1998). The user is referred to Zotarelli et al., (2009) for a synthetic, step-by-step guide on how the calculation is performed.

def FAO_56_ET0(Tmax,Tmin,Tmean,relhum,ws,lat_radians,elev,rad,albedo_crop):
    #wind component
    es_Tx = 0.6108 * np.exp((17.27 * Tmax[day,lat,lon]) / (Tmax[day,lat,lon] + 237.3)) 
    es_Tn = 0.6108 * np.exp((17.27 * Tmin[day,lat,lon]) / (Tmin[day,lat,lon] + 237.3))
    es = ((es_Tx + es_Tn) / 2) #mean sat vapour pressure
    ea = relhum[day,lat,lon]/100 * es  # rel_hum as fraction
    P_barometrica = 101.3 * ((293 - 0.0065 * elev[lat,lon]) / 293)**5.26# Psychrometric constant, [kPa C-1]
    gamma = (P_barometrica * 0.001013) / (0.622 * 2.45)
    Delta = 4098 * (0.6108 * np.exp((17.27 * Tmean[day,lat,lon]) / (Tmean[day,lat,lon] + 237.3))) / ((Tmean[day,lat,lon] + 237.3)**2) # Slope of Saturation Vapor Pressure [kPa C-1]
    #thermal radiation component
    #standard evapotranspiration potential
    ET0 = 0.408*Rnet+DT + PT*TT*(es-ea) 
    return ET0

Yield percentage reduction#

The cells below combine the ET0 calculation results with the precipitation data and the crop specific parameters to compute the percentage yield loss from irrigation deficit. The main calculation steps are described within the cell.


The calculation might take a few minutes, especially if a large region is selected.

First we need to create a dataframe to store the yield loss results. This will contain information about the grid cell coordinates, the selected RCP scenario and the start and end years of the studied period. This dataframe will be saved as a .csv spreadsheet for you to explore and to serve as input for the risk assessment workflow.

#create dataframe where results will be stored

initial_SM = 0   # Initial Soil Moisture, default is 0 but can be edited by expert user 
sowing_SM = 0   # Sowing Soil Moisture, default is 0 but can be edited by expert user 
deficit_start= 0  # Water deficit at the beginning first day, default is 0 but can be edited by expert user
deficit_final= 0  # Water deficit at the end of first day, default is 0 but can be edited by expert user

zm = 10  # height of wind measurements [m]
zh = 2  # height of humidity measurements [m]
h = 0.12  # crop height [m]
albedo_crop = 0.23  #constant
G=0 #for daily time-steps in Et0 calculation

#Latitude in radians

#The for loop allows to perform the yield loss calculation for all grid-cells for 365 days
for crop in crop_list:
    for lat in np.arange(len(Tmax[0,:,0])):
        for lon in np.arange(len(Tmin[0,0,:])):
            for day in range(1, 365):
                i = crop_list.index(crop)
                row = crop_table[(crop_table["Crop"] == crop) & (crop_table["Clim"] == climate_zone[lat,lon])]

                # Defining Length of Growing Periods, Kc values, and Rooting Depth
                season_start = row.iloc[0, 10]
                season_end = row.iloc[0, 11]

                if season_start < season_end:
                    LGP = season_end - season_start + 1
                    LGP = season_end - season_start + 366

                # Definition of kc
                kc_in = row.iloc[0, 2]
                kc_mid = row.iloc[0, 3]
                kc_end = row.iloc[0, 4]

                # Definition of growing period phases
                L1 = round(LGP * row.iloc[0, 5])
                L2 = round(LGP * row.iloc[0, 6])
                L3 = round(LGP * row.iloc[0, 7])
                L4 = LGP - (L1 + L2 + L3)

                L12 = L1 + L2
                L123 = L1 + L2 + L3
                L2_den = L2
                L4_den = L4
                L12_den = L12

                s = season_start
                h = season_end

                # Definition of rooting depth
                RD = row.iloc[0, [12, 13]]
                # Calculation of ET0
            # Determination of effective rainfall (prec_eff) folowing Mancosu et al., (2015)
                if (h >= s) & (day >= s):
                    prec_eff=np.where(precipitation[day,lat,lon]  < (ETa_irr[lat,lon,i,day-1]+deficit_final),
                elif (h >= s) & (day < s):
                    ETa_irr[lat,lon,i,day] = 0 # if day=1 there is no ETa  
                elif (h<s) & (day<s) & (day>h):
                    ETa_irr[lat,lon,i,day] = 0 # if day=1 there is no ETa  
                    prec_eff=np.where(precipitation[day,lat,lon]  < (ETa_irr[lat,lon,i,day-1]+deficit_final),

                #1st step: Calculation of crop coefficient Kc   
                kc = np.where((h > s) & (day >= s) & (day < (s + L1 - 1)), kc_in, 0) + np.where(
                    (h > s) & (day >= (s + L1 - 1)) & (day < (s + L12 - 1)), kc_in + ((day - (s + L1) + 1) / L2_den) * (kc_mid - kc_in), 0) + np.where(
                    (h > s) & (day >= (s + L12 - 1)) & (day < (s + L123 - 1)), kc_mid, 0) + np.where(
                    (h > s) & (day >= (s + L123 - 1)) & (day <= h),kc_mid + ((day - (s + L123) + 1) / L4_den) * (kc_end - kc_mid), 0) + np.where(
                    (h < s) & (day >= s) & (day < (s + L1 - 1)), kc_in, 0) + np.where(
                    (h < s) & (day < (h - (LGP - L1))), kc_in, 0) + np.where(
                    (h < s) & (day >= (s + L1 - 1)) & (day < (s + L12 - 1)),kc_in + ((day - (s + L1) + 1) / L2_den) * (kc_mid - kc_in), 0)+ np.where(
                    (h < s) & (day >= (h - (LGP - L1))) & (day < (h - (LGP - L12))), kc_in + ((366 - s - L1 + day) / L2_den) * (kc_mid - kc_in), 0) + np.where(
                    (h < s) & (day >= (s + L12 - 1)) & (day < (s + L123 - 1)), kc_mid, 0) + np.where(
                    (h < s) & (day >= (h - (LGP - L12))) & (day < (h - (LGP - L123))), kc_mid, 0) + np.where(
                    (h < s) & (day >= (s + L123 - 1)) & (day <= (s + LGP - 1)),kc_mid + ((day - (s + L123) + 1) / L4_den) * (kc_end - kc_mid), 0) +np.where(
                    (h < s) & (day >= (h - (LGP - L123))) & (day <= h),kc_mid + ((366 - s - L123 + day) / L4_den) * (kc_end - kc_mid), 0)


                # 2nd step: Increment of Soil Moisture due to Precipitation input
                if row.iloc[0, 15] ==1: # Annual crops
                    daily_rooting_depth = (np.where((s < h) & (s <= day) & (s + L12 - 1 > day),(RD[1] - RD[0]) / L12_den * (day - s + 1) + RD[0],0)+
                                           np.where((s < h) and (s + L12 - 1 <= day) and (day <= h),RD[1],0)+
                                           np.where((s > h) and (s <= day) and (s + L12 - 1 > day),((RD[1] - RD[0]) / L12_den) * (day - s + 1) + RD[0],0)+
                                           np.where((s > h) and (h - LGP - L12>day),((RD[1] - RD[0]) / L12_den) * (day - s + 1) + RD[0],0)+
                                           np.where((s > h) and (s + L12 - 1<= day),RD[1],0)+
                                           np.where((s > h) and (h >= day) and (h - LGP - L12 <= day), RD[1],0))
                    deficit_meteo = np.where(day == s, (AWC[lat,lon] - sowing_SM) * 1000 * daily_rooting_depth + 
                                             (AWC[lat,lon] - initial_SM) * 1000 * daily_rooting_depth,(AWC[lat,lon] - initial_SM) * 1000 * daily_rooting_depth - prec_eff)
                else: # Perennial crops
                    daily_rooting_depth = RD[1]
                    deficit_meteo = (AWC[lat,lon] - initial_SM) * 1000 * daily_rooting_depth - prec_eff 

            # 3rd step: Calculation of daily maximum water capacity in the rooting zone (TAW) and amount of water available until water stress occurs (RAW)
                TAW = AWC[lat,lon] * daily_rooting_depth 
                RAW = AWC[lat,lon] * daily_rooting_depth * row.iloc[0, 14] 

            # 4th step: Calculation of water-stress coefficient
                ks = np.where((TAW>0) & (TAW - deficit_meteo > 0), (TAW - deficit_meteo) / (TAW - RAW), 0)
                ks[ks > 1] = 1

            # 5th step: Calculation of crop Standard (ETc) and Actual (ETa) Evapotranspiration potential
                ETc = ET0*kc 
                ETa = ET0* kc * ks
                ETa_rf[lat,lon,i,day] = ETa
                ETa_irr[lat,lon,i,day]= ETc 

            # 6th step: Reduction of Soil Moisture due to Actual Evapotranspiration
                deficit_final = np.where((deficit_meteo >= 0) & (deficit_meteo <= TAW),deficit_meteo + ETa,

                deficit_final_limit = TAW - ks * (TAW - RAW)  # Deficit_max = RAW

            # 7th step: irrigation requirement, expressed as difference in maximum vs actual evapotranspiration
                irrigation = np.where(ETa_rf[lat,lon,i,day] <= ETa_irr[lat,lon,i,day],
                                  ETa_irr[lat,lon,i,day] - ETa_rf[lat,lon,i,day], 0)

            # wat_soil= additional water required to bring soil moisture at required amount
                wat_soil = np.where(((deficit_final - deficit_final_limit) >= 0)
                                & ((deficit_final - deficit_final_limit) >= ETa_irr[lat,lon,i,day]),
                                deficit_final - deficit_final_limit - irrigation, 0)

            # 8° step: total irrigation requirement
                Irr_req[lat,lon,i,day] = irrigation + wat_soil

            # 9° step: calculation of final soil moisture at the end of the day to restart the loop
                final_SM = AWC[lat,lon] - ((deficit_final - Irr_req[lat,lon,i,day]) / (1000 * daily_rooting_depth))
                initial_SM = final_SM

    # Yield water response calculation based on FAO I&D n.33 (Doorenbos et al., 1979). 
            ky = row.iloc[0, 16]
            #Yield loss expressed as percentage loss from maximum yield 
            yield_loss_perc= (ky * (1 - annual_ETa_rf / annual_ETa_irr))*100 

# Save the output as a spreadsheet for visualization and a numpy file to be used in the risk assessment

Plotting the results#

The cells below allow to plot the results of the hazard assessment. First, the AWC, precipitation and ET0 are plotted to help the user contestualise the yield reduction results. Precipitation and ET0 are summed across the selected growing season so that the results displayed show the total rainfall availability and standard evapotranspiration demand. The plotting procedure steps are described within the cell.


Areas showing higher cumulate rainfall and lower cumulate standard evapotranspiration will be less likely to experience severe yield losses. However, it is important to remember tha the standard evapotranspiration is only based on climate conditions and does not take into account crop-specific paramters (introduced in ETa) that might significantly change the final yield change results.

rcp=ds_tasmax.attrs['experiment_id'] #identify the climate projection
ystart=ds_tasmax['time.year'][0].to_numpy() #identify the projection start year
yend=ds_tasmax['time.year'][-1].to_numpy() #identify the projection end year
#zoom parameter

#define the longitude and latitude coordinates

#plot soil available water capacity
map_AWC =Basemap(
               llcrnrlat=bbox[1]-zoom, urcrnrlat=bbox[3]+zoom,
               llcrnrlon=bbox[0]-zoom, urcrnrlon=bbox[2]+zoom,

map_AWC.colorbar().set_label('Soil Available Water Capacity (mm)',fontsize=13)
map_AWC.readshapefile(f'{data_dir}/{nuts_name}',nuts_name,color='b',linewidth=2) #draw region boundaries from shapefile
plt.title('AWC '+str(nuts.iloc[0,4]),fontweight='bold') #title
plt.savefig(f'{results_dir}/'+str(nuts.iloc[0,4])+'_AWC.png') #save plots in results directory

#plot cumulate precipitation through the selected growing season
map_pr =Basemap(
               llcrnrlat=bbox[1]-zoom, urcrnrlat=bbox[3]+zoom,
               llcrnrlon=bbox[0]-zoom, urcrnrlon=bbox[2]+zoom,

map_pr.colorbar().set_label('Precipitation intensity (mm)',fontsize=13)
map_pr.readshapefile(f'{data_dir}/{nuts_name}',nuts_name,color='b',linewidth=2) #draw region boundaries from shapefile
plt.suptitle('Cumulate precipitation intensity through the growing season',fontsize=13,fontweight='bold')
plt.title(str(nuts.iloc[0,4]+' '+rcp+' '+str(ystart)+'-'+str(yend))) #title
plt.savefig(f'{results_dir}/'+str(nuts.iloc[0,4])+'_AWC.png') #save plots in results directory

#plot cumulate ET0 through the selected growing season
map_et0 =Basemap(
               llcrnrlat=bbox[1]-zoom, urcrnrlat=bbox[3]+zoom,
               llcrnrlon=bbox[0]-zoom, urcrnrlon=bbox[2]+zoom,

map_et0.colorbar().set_label('Standard evapotranspiration (mm)',fontsize=13)
map_et0.readshapefile(f'{data_dir}/{nuts_name}',nuts_name,color='b',linewidth=2) #draw region boundaries from shapefile
plt.suptitle('Cumulate standard evapotranspiration (ET0) through the growing season',fontsize=13,fontweight='bold')
plt.title(str(nuts.iloc[0,4]+' '+rcp+' '+str(ystart)+'-'+str(yend))) #title
plt.savefig(f'{results_dir}/'+str(nuts.iloc[0,4])+'_AWC.png') #save plots in results directory
../../../../_images/0cde99ba4069c6e664217a4be7c65ee43d39940aa852d30f6009714d883f046a.png ../../../../_images/acf063d81baa810cb046ab8b7d979a85bf5c257d7de97f3607b5b9555c2011f2.png ../../../../_images/edd480299561d2104088073f928fec475d2654cc767281133fda964128652504.png

The yield loss results for the studied region are displayed using the cell below. The resulting plot states the crop, RCP scenario and reference period used in the assessment.


Use the ‘zoom’ paramerter to set how much you would like the final plot to be zoomed-out from the region boundaries (0=no zoom out,1=100 km). Here a zoom of 0.5 degrees (50 km) is used.

#plot the yield loss
for a in np.arange(len(crop_list)):
    #create the basemap around the region boundaries
    map_yield =Basemap(
               llcrnrlat=bbox[1]-zoom, urcrnrlat=bbox[3]+zoom,
               llcrnrlon=bbox[0]-zoom, urcrnrlon=bbox[2]+zoom,

    #plot the yield loss
    map_yield.colorbar().set_label('Yield loss (%)',fontsize=13)
    map_yield.readshapefile(f'{data_dir}/{nuts_name}',nuts_name,color='b',linewidth=2) #draw region boundaries from shapefile
    plt.suptitle(crop_list[a].title()+' yield loss from precipitation deficit',fontsize=13,fontweight='bold') #title
    plt.title(str(nuts.iloc[0,4])+' '+rcp+' '+str(ystart)+'-'+str(yend)) #subtitle
    plt.savefig(f'{results_dir}/'+str(nuts.iloc[0,4])+'_'+str(crop_list[a])+'_yield_loss.png') #save plots in results directory
../../../../_images/a8f0bca2f1bc1b2ddec3ca7a302e5fcb3293a8f2af51765b96b0deda19e64968.png ../../../../_images/4118264699677706e361d5d157c595f86d718fca43194bfbbc5db8bdfdd86f62.png

The maps show the potential yield loss from irrigation deficit in the studied region for the selected crops (here maize and wheat), emission scenario (here RCP2.6) and period (here 2046-2050). Yield loss is expressed as the percentage reduction in yield if crops are grown under rainfed only conditions rather than fully irrigated.

These maps can be used by demonstrators to understand in which areas of their region cropland will potentially experience the highest hydro-climatic stress. This allows them to target adaptation efforts in the most affected areas favouring a cost-effective use of resources. At the same time, the map provides a snap-shot of a potential future growing season that can be used to guide cropland expansion towards areas less affected by water stress.


In this workflow we learned:

  • how to download and explore high-resolution climate, elevation and available water capacity data for your region.

  • how to use and combine these datasets to calculate the actual evapotranspiration potential of different crops using the FAO56 methodology.

  • how to derive potential yield losses of different crops without full irrigation from the actual evapotranspiration deficit.

The yield loss results produced in this workflow have been stored in .csv and .npy formats to be explored and later used in the risk assessment workflow. The yield loss results can also be visualised in the maps saved locally.


Euro-Mediterranean Center on Climate Change (CMCC), Italy.

Author of the workflow: Andrea Rivosecchi


  • Allen, R.G., Pereira, L.S., Raes, D. & Smith, M. (1998) Crop evapotranspiration-Guidelines for computing crop water requirements-FAO Irrigation and drainage paper 56. Fao, Rome. 300 (9), D05109.

  • Danielson, J.J. & Gesch, D.B. (2011) Global multi-resolution terrain elevation data 2010 (GMTED2010): U.S. Geological Survey Open-File Report 2011–1073, p.26.

  • Doorenbos, J. & Kassam, A.H. (1979) Yield response to water. Irrigation and drainage paper. 33, 257.

  • Hengl, T. & Gupta S. (2019). Soil available water capacity in mm derived for 5 standard layers (0-10, 10-30, 30-60, 60-100 and 100-200 cm) at 250 m resolution (v0.1). Zenodo. https://doi.org/10.5281/zenodo.2629149

  • Mancosu, N., Spano, D., Orang, M., Sarreshteh, S. & Snyder, R.L. (2016) SIMETAW# - a Model for Agricultural Water Demand Planning. Water Resources Management. 30 (2), 541–557. doi:10.1007/s11269-015-1176-7.

  • Steduto, P., Hsiao, T.C., Fereres, E. & Raes, D. (2012) Crop yield response to water. FAO Rome.

  • Van Velthuizen, H., Huddleston, B., Fischer, G., Salvatore, M., Ataman, E., Nachtergaele, F. O., Zanetti, M., Bloise, M., Antonicelli, A., Bel, J., De Liddo, A., De Salvo, P., & Franceschini, G. (2007). Mapping biophysical factors that influence agricultural production and rural vulnerability (Issue 11). Food & Agriculture Org.

  • Zotarelli, L., Dukes, M. D., Romero, C. C., Migliaccio, K. W., & Morgan, K. T. (2009). Step by step calculation of the Penman-Monteith Evapotranspiration (FAO-56 Method).