Water Productivity assessment

Download and process Level 3 data from WaPOR - Gezira Irrigation Scheme

Notebook prepared by Ben Maathuis, Bas Retsios and Willem Nieuwenhuis. ITC-University of Twente, Enschede. The Netherlands

Making use of WAPOR-DL (https://bitbucket.org/cioapps/wapordl/src/main/)

WaPOR data needed for Water Productivity assessment will be downloaded and processed for a part of the Gezira Scheme in Sudan for which Level 3 data is available. The crop selected is sorghum and the crop mask extraction is based on the information acquired during a field survey which was conducted in 2019. The sorghum crop mask will be extracted using some general assumptions based on the field experiences gained.

AETI, T and NPP data is downloaded from the WaPOR data repository (using WAPOR-DL) and is converted to the correct unit, aggregated over the crop growing season and Land (biomass and crop yield) and (gross and net) Water Productivity is assessed, followed by some other irrigation performance indicators like Uniformity, Beneficial fraction and Relative Water Deficit. The formula's applied are derived from the Water Productivity in Practice Project (Water-PIP).

All data is stored locally (under the notebook folder) and is available in Ilwis raster format. ILWIS386 can be used for visualization and the latest version can be downloaded from: https://filetransfer.itc.nl/pub/52n/ILWIS386/Software/. The installation manual is available from: https://filetransfer.itc.nl/pub/52n/ILWIS386/Tutorial/

In the Gezira Irrigation Scheme, the average sorghum yield is estimated to be around 1.212 tons per hectare (ha). Studies suggest that this yield can be significantly higher, reaching a potential yield of 2000-4000 kg/ha. Furthermore some reference figures from different sources:

• The Gezira scheme's sorghum crop experiences an average seasonal transpiration of 322 mm and an average seasonal evapotranspiration of 496 mm.
• The average seasonal gross crop water productivity for sorghum is 0.43 kg/m3.
• The scheme's sorghum yield can vary significantly across different locations, with El Naseeh block in the middle of the Managil major canal showing a yield of 1.212 tons/ha.

Challenges: The Gezira scheme has faced challenges like low productivity, poor water management, and irrigation mismanagement. A steady decline in sorghum yield was observed since the 1970s, with yields decreasing from 744.3 kg/ha to 476.6 kg/ha since 1982.

Some further information can be obtained from:

• https://www.researchgate.net/publication/299061128_Remote_Sensing_Derived_Crop_Coefficient_for_Estimating_Crop_Water_Requirements_for_Irrigated_Sorghum_in_the_Gezira_Scheme_Sudan
• https://www.mdpi.com/2073-4395/15/1/110
• https://github.com/eLEAF-Github/WAPORACT Water Productivity in Practice project (WaterPIP)

#if WaPOR-DL is not installed uncomment the line below and execute this field
#!pip install wapordl

#load libraries
import wapordl
import matplotlib.pyplot as plt
from matplotlib.colors import ListedColormap, BoundaryNorm
from IPython.display import HTML
from datetime import datetime, timedelta
import ilwis
import numpy as np
import sys
import os
import glob
import warnings

INFO: WaPORDL (`1.1.2`)

Get some general information on WaPOR level 3 areas and level 1, 2 and 3 variables

#uncomment the line below to get a listing of the level 3 areas
#here we will use the data from: 'GEZ': 'Gezira, Sudan'
wapordl.region_selector.l3_codes()

{'VTM': 'Tay Nguyen Pilot Study Region',
 'MAL': 'Malwathu Oya West Sub Catchment, Sri Lanka',
 'ODN': 'Office du Niger, Mali',
 'GEZ': 'Gezira, Sudan',
 'PAL': 'Jericho, Palestine',
 'SED': 'Senegal Delta, Senegal',
 'MIT': 'Mitidja, Algeria',
 'JVA': 'North Jordan Valley, Jordan',
 'JEN': 'Jendouba, Tunisia',
 'KTB': 'Tana and Bura irrigation schemes, Kenya',
 'LDA': 'Waddan region, Libya',
 'KMW': 'Mwea irrigation scheme, Kenya',
 'SAN': 'Sanaa basin, Yemen',
 'LOT': 'Tarhona region, Libya',
 'LOU': 'Moulay Bousselham, Morocco',
 'BKA': 'Bekaa, Lebanon',
 'ENO': 'Northern Egypt',
 'AWA': 'Awash, Ethiopia',
 'BUS': 'Busia, Kenya',
 'ERB': 'Erbil, Iraq',
 'GAR': 'West Gharraf, Iraq',
 'NAJ': 'Najaf, Iraq',
 'KOG': 'Koga, Ethiopia',
 'LCE': 'Fezzan region, Libya',
 'MBL': 'Baixo Limpopo, Mozambique',
 'LAK': 'Lower Akagera, Rwanda',
 'KWL': 'Khanewal, Pakistan',
 'YAN': 'Yanze catchment, Rwanda',
 'ZAN': 'Zankalon, Egypt',
 'JAF': 'Jafr-Shoubak, Jordan',
 'MAG': 'Magdalena, Colombia',
 'LAM': 'Lamego, Mozambique',
 'MUV': 'Muvumba catchment, Rwanda',
 'KAI': 'Kairouan, Tunisia',
 'SNG': 'Sanghar, Pakistan'}

#uncomment the line below to get a listing of the level 1, 2 and 3 variables
#here we use level 3
wapordl.variable_descriptions.WAPOR3_VARS

{'L1-AETI-A': {'long_name': 'Actual EvapoTranspiration and Interception',
  'units': 'mm/year'},
 'L1-AETI-D': {'long_name': 'Actual EvapoTranspiration and Interception',
  'units': 'mm/day'},
 'L1-AETI-M': {'long_name': 'Actual EvapoTranspiration and Interception',
  'units': 'mm/month'},
 'L1-E-A': {'long_name': 'Evaporation', 'units': 'mm/year'},
 'L1-E-D': {'long_name': 'Evaporation', 'units': 'mm/day'},
 'L1-GBWP-A': {'long_name': 'Gross Biomass Water Productivity',
  'units': 'kg/m³'},
 'L1-I-A': {'long_name': 'Interception', 'units': 'mm/year'},
 'L1-I-D': {'long_name': 'Interception', 'units': 'mm/day'},
 'L1-NBWP-A': {'long_name': 'Net Biomass Water Productivity',
  'units': 'kg/m³'},
 'L1-NPP-D': {'long_name': 'Net Primary Production', 'units': 'gC/m²/day'},
 'L1-NPP-M': {'long_name': 'Net Primary Production', 'units': 'gC/m²/month'},
 'L1-PCP-A': {'long_name': 'Precipitation', 'units': 'mm/year'},
 'L1-PCP-D': {'long_name': 'Precipitation', 'units': 'mm/day'},
 'L1-PCP-E': {'long_name': 'Precipitation', 'units': 'mm/day'},
 'L1-PCP-M': {'long_name': 'Precipitation', 'units': 'mm/month'},
 'L1-QUAL-LST-D': {'long_name': 'Quality Land Surface Temperature',
  'units': 'd'},
 'L1-QUAL-NDVI-D': {'long_name': 'Quality of Normalized Difference Vegetation Index',
  'units': 'd'},
 'L1-RET-A': {'long_name': 'Reference Evapotranspiration', 'units': 'mm/year'},
 'L1-RET-D': {'long_name': 'Reference Evapotranspiration', 'units': 'mm/day'},
 'L1-RET-E': {'long_name': 'Reference Evapotranspiration', 'units': 'mm/day'},
 'L1-RET-M': {'long_name': 'Reference Evapotranspiration',
  'units': 'mm/month'},
 'L1-RSM-D': {'long_name': 'Relative Soil Moisture', 'units': '%'},
 'L1-T-A': {'long_name': 'Transpiration', 'units': 'mm/year'},
 'L1-T-D': {'long_name': 'Transpiration', 'units': 'mm/day'},
 'L1-TBP-A': {'long_name': 'Total Biomass Production', 'units': 'kg/ha'},
 'L2-AETI-A': {'long_name': 'Actual EvapoTranspiration and Interception',
  'units': 'mm/year'},
 'L2-AETI-D': {'long_name': 'Actual EvapoTranspiration and Interception',
  'units': 'mm/day'},
 'L2-AETI-M': {'long_name': 'Actual EvapoTranspiration and Interception',
  'units': 'mm/month'},
 'L2-E-A': {'long_name': 'Evaporation', 'units': 'mm/year'},
 'L2-E-D': {'long_name': 'Evaporation', 'units': 'mm/day'},
 'L2-GBWP-A': {'long_name': 'Gross Biomass Water Productivity',
  'units': 'kg/m³'},
 'L2-I-A': {'long_name': 'Interception', 'units': 'mm/year'},
 'L2-I-D': {'long_name': 'Interception', 'units': 'mm/day'},
 'L2-NBWP-A': {'long_name': 'Net Biomass Water Productivity',
  'units': 'kg/m³'},
 'L2-NPP-D': {'long_name': 'Net Primary Production', 'units': 'gC/m²/day'},
 'L2-NPP-M': {'long_name': 'Net Primary Production', 'units': 'gC/m²/month'},
 'L2-QUAL-NDVI-D': {'long_name': 'Quality of Normalized Difference Vegetation Index',
  'units': 'd'},
 'L2-RSM-D': {'long_name': 'Relative Soil Moisture', 'units': '%'},
 'L2-T-A': {'long_name': 'Transpiration', 'units': 'mm/year'},
 'L2-T-D': {'long_name': 'Transpiration', 'units': 'mm/day'},
 'L2-TBP-A': {'long_name': 'Total Biomass Production', 'units': 'kg/ha'},
 'L3-AETI-A': {'long_name': 'Actual EvapoTranspiration and Interception',
  'units': 'mm/year'},
 'L3-AETI-D': {'long_name': 'Actual EvapoTranspiration and Interception',
  'units': 'mm/day'},
 'L3-AETI-M': {'long_name': 'Actual EvapoTranspiration and Interception',
  'units': 'mm/month'},
 'L3-E-A': {'long_name': 'Evaporation', 'units': 'mm/year'},
 'L3-E-D': {'long_name': 'Evaporation', 'units': 'mm/day'},
 'L3-GBWP-A': {'long_name': 'Gross Biomass Water Productivity',
  'units': 'kg/m³'},
 'L3-I-A': {'long_name': 'Interception', 'units': 'mm/year'},
 'L3-I-D': {'long_name': 'Interception', 'units': 'mm/day'},
 'L3-NBWP-A': {'long_name': 'Net Biomass Water Productivity',
  'units': ' kg/m³'},
 'L3-NPP-D': {'long_name': 'Net Primary Production', 'units': 'gC/m²/day'},
 'L3-NPP-M': {'long_name': 'Net Primary Production', 'units': 'gC/m²/month'},
 'L3-QUAL-NDVI-D': {'long_name': 'Quality of Normalized Difference Vegetation Index',
  'units': 'd'},
 'L3-RSM-D': {'long_name': 'Relative Soil Moisture', 'units': '%'},
 'L3-T-A': {'long_name': 'Transpiration', 'units': 'mm/year'},
 'L3-T-D': {'long_name': 'Transpiration', 'units': 'mm/day'},
 'L3-TBP-A': {'long_name': 'Total Biomass Production', 'units': 'kg/ha'}}

Season for Sorghum

• SOS = dekad 3 July
• EOS = dekad 3 November

Specify the year to be selected in the field below.

year = '2020'

#create a few folders to store the data downloaded / processed
folder = os.getcwd()+r'/waporl3_download'
os.chdir(".")
print("current dir is: %s" % (os.getcwd()))

if os.path.isdir(folder):
    print("Folder exists")
else:
    print("Folder doesn't exists")
    os.mkdir(folder)

current dir is: d:\jupyter\notebook_scripts\Special\WaPOR_Level3
Folder exists

#create a few folders to store the data downloaded / processed
folder = os.getcwd()+ r'/waporl3_'+year
os.chdir(".")
print("current dir is: %s" % (os.getcwd()))

if os.path.isdir(folder):
    print("Folder exists")
else:
    print("Folder doesn't exists")
    os.mkdir(folder)

current dir is: d:\jupyter\notebook_scripts\Special\WaPOR_Level3
Folder exists

#set sample data folder
#use the Python function os.getcwd() to retrieve the current folder used and add the apporriate sub-folder
download_dir = os.getcwd() + r'/waporl3_download'
working_dir = os.getcwd() +r'/waporl3_'+year

#set the working directory for ILWISPy
ilwis.setWorkingCatalog(working_dir)
print(ilwis.version())
print(download_dir)
print(working_dir)

1.0 build 20250513
d:\jupyter\notebook_scripts\Special\WaPOR_Level3/waporl3_download
d:\jupyter\notebook_scripts\Special\WaPOR_Level3/waporl3_2020

Retrieve time series data for AETI, NPP and T

Retrieve and process AETI

region = "GEZ"
folder = download_dir
variable = "L3-AETI-D"
period = [year+"-07-21", year+"-11-21"]
unit = "dekad" # or choose "day", "month", "year", "none" (default).
overview = -1 #or "None"

aeti_df = wapordl.wapor_ts(region, variable, period, overview=overview, unit_conversion=unit)
aeti_fp = wapordl.wapor_map(region, variable, period, folder, overview=overview, unit_conversion=unit)

INFO: Found 13 files for L3-AETI-D.
INFO: Converting units from [mm/day] to [mm/dekad] (use_xarray = False).                                               
INFO: Consider installing `xarray`, `rioxarray` and `dask` for faster unit conversions.
INFO: Found 13 files for L3-AETI-D.
INFO: Converting units from [mm/day] to [mm/dekad] (use_xarray = False).                                               
INFO: Consider installing `xarray`, `rioxarray` and `dask` for faster unit conversions.

aeti_df

	minimum	maximum	mean	start_date	end_date	number_of_days
0	0.0	63.8	2.2415	2020-07-21	2020-07-31	11 days
1	0.0	39.0	2.9257	2020-08-01	2020-08-10	10 days
2	0.0	45.0	7.0998	2020-08-11	2020-08-20	10 days
3	0.0	60.5	14.9136	2020-08-21	2020-08-31	11 days
4	0.0	68.0	19.3251	2020-09-01	2020-09-10	10 days
5	0.0	73.0	26.0379	2020-09-11	2020-09-20	10 days
6	0.0	77.0	30.2474	2020-09-21	2020-09-30	10 days
7	0.0	77.0	34.5848	2020-10-01	2020-10-10	10 days
8	0.0	80.0	36.4419	2020-10-11	2020-10-20	10 days
9	0.0	96.8	36.1667	2020-10-21	2020-10-31	11 days
10	0.0	96.0	29.7930	2020-11-01	2020-11-10	10 days
11	0.0	97.0	28.1959	2020-11-11	2020-11-20	10 days
12	0.0	90.0	27.8646	2020-11-21	2020-11-30	10 days

aeti_df.attrs

{'long_name': 'Actual EvapoTranspiration and Interception',
 'original_units': 'mm/day',
 'overview': '-1',
 'units': 'mm/dekad'}

aeti_fp

'd:\\jupyter\\notebook_scripts\\Special\\WaPOR_Level3/waporl3_download\\GEZ_L3-AETI-D_NONE_dekad_converted.tif'

# rename downloaded file
old_aeti = aeti_fp
new_aeti = download_dir+ '/GEZ_L3-AETI-D_NONE_dekad_converted_'+year+'.tif'

# Rename the file
os.rename(old_aeti, new_aeti)

ilwis.setWorkingCatalog(working_dir)
print(working_dir)

d:\jupyter\notebook_scripts\Special\WaPOR_Level3/waporl3_2020

aeti_ds = ilwis.RasterCoverage(new_aeti)

#check meta data
print(aeti_ds.size())
print(aeti_ds.envelope())
coordSys = aeti_ds.coordinateSystem()
coordSys.toWKT()

Size(2318, 1523, 13)
516820.000000 1567400.000000 563180.000000 1597860.000000

'PROJCS["gez_l3-aeti-d_none_dekad_converted_2020.tif",GEOCS["gez_l3-aeti-d_none_dekad_converted_2020.tif",DATUM[" WGS 84",[DWGS84],ELLIPSOID["WGS 84",6378137.000000000000,298.257223563000],PRIMEM["Greenwich",0, AUTHORITY["EPSG",8901"]]],PROJECTION["Transverse_Mercator"],PARAMETER["scale_factor",0.9996],PARAMETER["false_easting",500000],PARAMETER["false_northing",0],PARAMETER["scale",0.9996],PARAMETER["latitude_of_origin",0],PARAMETER["central_meridian",33],UNIT[meter,1.0]]'

#Note pandas df for date of rasterband selected (2021-09-11	2021-09-20)
ts1 = ilwis.do('selection',aeti_ds,"rasterbands(5)")
ts1 = ilwis.do('mapcalc', 'iff(@1==?, -10, @1*0.1)',ts1) #*apply scaling, see also max of layer 05!

ts_2np = np.fromiter(iter(ts1), np.float64, ts1.size().linearSize()) 
ts_2np = ts_2np.reshape((ts1.size().ysize, ts1.size().xsize))

min_val = np.nanmin(ts_2np)
max_val = np.nanmax(ts_2np)
print(min_val)
print(max_val)

-10.0
73.0

# Create custom colormap

cmap = plt.cm.jet
cmap_colors = cmap(np.linspace(0, 1, 256))

# Add white at the beginning for the -1 value
new_colors = np.vstack(([1, 1, 1, 1], cmap_colors))
new_cmap = ListedColormap(new_colors)

# Create bounds and norm
max_val = np.nanmax(ts_2np[ts_2np >= 0])  # Ignore -1 in max calc
bounds = np.linspace(0, max_val, num=258)  # One extra bin for -1
norm = BoundaryNorm(boundaries=bounds, ncolors=new_cmap.N)

# Plot
plt.imshow(ts_2np, cmap=new_cmap, norm=norm)
plt.colorbar(label='Data Value', shrink=0.5)
#plt.title('Custom Colormap with -1 as White')
plt.title('AETI-ts05 (mm/dek)')
plt.show()

ts_aeti = ilwis.do('mapcalc', 'iff(@1==?,-10, @1*0.1)',aeti_ds) #set background to value of -10 and apply scaling factor to all layers
ts_aeti.store('gez_aeti.mpl')

Retrieve and process NPP

region = "GEZ"
folder = download_dir
variable = "L3-NPP-D"
period = [year+"-07-21", year+"-11-21"]
overview = -1 #or "None"
unit = "dekad" # or choose "day", "month", "year", "none" (default).

npp_df = wapordl.wapor_ts(region, variable, period, overview=overview, unit_conversion=unit)
npp_fp = wapordl.wapor_map(region, variable, period, folder, unit_conversion=unit)

INFO: Found 13 files for L3-NPP-D.
INFO: Converting units from [gC/m²/day] to [gC/m²/dekad] (use_xarray = False).                                         
INFO: Consider installing `xarray`, `rioxarray` and `dask` for faster unit conversions.
INFO: Found 13 files for L3-NPP-D.
INFO: Converting units from [gC/m²/day] to [gC/m²/dekad] (use_xarray = False).                                         
INFO: Consider installing `xarray`, `rioxarray` and `dask` for faster unit conversions.

npp_df.attrs

{'long_name': 'Net Primary Production',
 'original_units': 'gC/m²/day',
 'overview': '-1',
 'units': 'gC/m²/dekad'}

npp_df

	minimum	maximum	mean	start_date	end_date	number_of_days
0	0.0	23.166	0.459571	2020-07-21	2020-07-31	11 days
1	0.0	37.010	2.002801	2020-08-01	2020-08-10	10 days
2	0.0	34.370	3.760595	2020-08-11	2020-08-20	10 days
3	0.0	32.153	6.243658	2020-08-21	2020-08-31	11 days
4	0.0	32.120	5.937101	2020-09-01	2020-09-10	10 days
5	0.0	25.940	5.771858	2020-09-11	2020-09-20	10 days
6	0.0	24.830	6.425725	2020-09-21	2020-09-30	10 days
7	0.0	24.730	7.377189	2020-10-01	2020-10-10	10 days
8	0.0	24.300	7.529264	2020-10-11	2020-10-20	10 days
9	0.0	21.967	5.716068	2020-10-21	2020-10-31	11 days
10	0.0	28.500	6.213743	2020-11-01	2020-11-10	10 days
11	0.0	29.300	6.094072	2020-11-11	2020-11-20	10 days
12	0.0	37.410	8.114781	2020-11-21	2020-11-30	10 days

npp_fp

'd:\\jupyter\\notebook_scripts\\Special\\WaPOR_Level3/waporl3_download\\GEZ_L3-NPP-D_NONE_dekad_converted.tif'

# rename downloaded file
old_npp = npp_fp
new_npp = download_dir+ '/GEZ_L3-NPP-D_NONE_dekad_converted_'+year+'.tif'

# Rename the file
os.rename(old_npp, new_npp)

npp_ds = ilwis.RasterCoverage(new_npp)

#Note rasterband time step selected (2021-09-11	2021-09-20)
ts2 = ilwis.do('selection',npp_ds,"rasterbands(5)")
ts2 = ilwis.do('mapcalc', 'iff(@1==?, -10, @1*0.001)',ts2) #apply scaling factor, see max of layer 05!

ts_2np = np.fromiter(iter(ts2), np.float64, ts2.size().linearSize()) 
ts_2np = ts_2np.reshape((ts2.size().ysize, ts2.size().xsize))

min_val = np.nanmin(ts_2np)
max_val = np.nanmax(ts_2np)
print(min_val)
print(max_val)

-10.0
25.94

# Create custom colormap

cmap = plt.cm.jet
cmap_colors = cmap(np.linspace(0, 1, 256))

# Add white at the beginning for the -1 value
new_colors = np.vstack(([1, 1, 1, 1], cmap_colors))
new_cmap = ListedColormap(new_colors)

# Create bounds and norm
max_val = np.nanmax(ts_2np[ts_2np >= 0])  # Ignore -1 in max calc
bounds = np.linspace(0, max_val, num=258)  # One extra bin for -1
norm = BoundaryNorm(boundaries=bounds, ncolors=new_cmap.N)

# Plot
plt.imshow(ts_2np, cmap=new_cmap, norm=norm)
plt.colorbar(label='Data Value', shrink=0.5)
#plt.title('Custom Colormap with -1 as White')
plt.title('NPP-ts05 (gC/m2/dek)')
plt.show()

ts_npp = ilwis.do('mapcalc', 'iff(@1==?, -10, @1*0.001)',npp_ds) #set background to value of -10 and apply scaling factor to all layers
ts_npp.store('gez_npp.mpl')

Retrieve and process T

region = "GEZ"
folder = download_dir
variable = "L3-T-D"
period = [year+"-07-21", year+"-11-21"]
overview = -1 #or "None"
unit = "dekad" # or choose "day", "month", "year", "none" (default).

t_df = wapordl.wapor_ts(region, variable, period, overview=overview, unit_conversion=unit)
t_fp = wapordl.wapor_map(region, variable, period, folder, overview=overview, unit_conversion=unit)

INFO: Found 13 files for L3-T-D.
INFO: Converting units from [mm/day] to [mm/dekad] (use_xarray = False).                                               
INFO: Consider installing `xarray`, `rioxarray` and `dask` for faster unit conversions.
INFO: Found 13 files for L3-T-D.
INFO: Converting units from [mm/day] to [mm/dekad] (use_xarray = False).                                               
INFO: Consider installing `xarray`, `rioxarray` and `dask` for faster unit conversions.

t_df.attrs

{'long_name': 'Transpiration',
 'original_units': 'mm/day',
 'overview': '-1',
 'units': 'mm/dekad'}

t_df

	minimum	maximum	mean	start_date	end_date	number_of_days
0	0.0	58.3	0.6784	2020-07-21	2020-07-31	11 days
1	0.0	35.0	1.2575	2020-08-01	2020-08-10	10 days
2	0.0	38.0	4.1427	2020-08-11	2020-08-20	10 days
3	0.0	53.9	10.8664	2020-08-21	2020-08-31	11 days
4	0.0	61.0	15.2986	2020-09-01	2020-09-10	10 days
5	0.0	64.0	20.7341	2020-09-11	2020-09-20	10 days
6	0.0	68.0	24.1808	2020-09-21	2020-09-30	10 days
7	0.0	68.0	26.5088	2020-10-01	2020-10-10	10 days
8	0.0	71.0	26.9833	2020-10-11	2020-10-20	10 days
9	0.0	81.4	26.5512	2020-10-21	2020-10-31	11 days
10	0.0	81.0	22.5422	2020-11-01	2020-11-10	10 days
11	0.0	82.0	21.0567	2020-11-11	2020-11-20	10 days
12	0.0	77.0	20.3343	2020-11-21	2020-11-30	10 days

t_fp

'd:\\jupyter\\notebook_scripts\\Special\\WaPOR_Level3/waporl3_download\\GEZ_L3-T-D_NONE_dekad_converted.tif'

# rename downloaded file
old_t = t_fp
new_t = download_dir+ '/GEZ_L3-T-D_NONE_dekad_converted_'+year+'.tif'

# Rename the file
os.rename(old_t, new_t)

t_ds = ilwis.RasterCoverage(new_t)

#Note rasterband time step selected (2021-09-01	2021-09-10)
ts7 = ilwis.do('selection',t_ds,"rasterbands(7)")
ts7 = ilwis.do('mapcalc', 'iff(@1==?, -10, @1*0.1)',ts7) #apply scaling factor, see max of layer 05!

ts_2np = np.fromiter(iter(ts7), np.float64, ts7.size().linearSize()) 
ts_2np = ts_2np.reshape((ts7.size().ysize, ts7.size().xsize))

min_val = np.nanmin(ts_2np)
max_val = np.nanmax(ts_2np)
print(min_val)
print(max_val)

-10.0
68.0

# Create custom colormap

cmap = plt.cm.jet
cmap_colors = cmap(np.linspace(0, 1, 256))

# Add white at the beginning for the -1 value
new_colors = np.vstack(([1, 1, 1, 1], cmap_colors))
new_cmap = ListedColormap(new_colors)

# Create bounds and norm
max_val = np.nanmax(ts_2np[ts_2np >= 0])  # Ignore -1 in max calc
bounds = np.linspace(0, max_val, num=258)  # One extra bin for -1
norm = BoundaryNorm(boundaries=bounds, ncolors=new_cmap.N)

# Plot
plt.imshow(ts_2np, cmap=new_cmap, norm=norm)
plt.colorbar(label='Data Value', shrink=0.5)
#plt.title('Custom Colormap with -1 as White')
plt.title('T-ts07 (mm/dek)')
plt.show()

ts_t = ilwis.do('mapcalc', 'iff(@1==?, -10, @1*0.1)',t_ds) #set background to value of -10 and apply scaling factor to all layers
ts_t.store('gez_t.mpl')

Aggregate timeseries of AETI, NPP and T to season sum

AETI_Season = ilwis.do('aggregaterasterstatistics', ts_aeti, 'sum')
AETI_Season.store('aeti_sum.mpr') # mm/season

NPP_Season = ilwis.do('aggregaterasterstatistics', ts_npp, 'sum')
NPP_Season.store('npp_sum.mpr') #gC/m2/season

T_Season = ilwis.do('aggregaterasterstatistics', ts_t, 'sum')
T_Season.store('t_sum.mpr') # mm/season

Create a Sorghum crop mask

Next to Sorghum other main crops are cotton, groundnut, pigeonpea and onion. Also there is wide range of other crops cultivated but in small areas such as okra, cowpea, sesame and white bean among others. Sorghum is planted in (mid-end) July and harvested in (end) November, which is destinctively different from the other main crops. To approximate the sorghum crop cycle from the Transpiration time series from the mid-season a minimum and maximum transpiration threshold (20 mm/dek / 70 mm/dek) is expected, at the start of the season - initial stage an increase in transpiration, from mid-season towards end of season a decrease in Transpiration and during the late season Transpiration is less than 20 mm/dek.

#selection of timesteps to be used in the sorghum crop mask creation (from Transpiration cycle)
t1= ilwis.do('selection',ts_t,"rasterbands(0")
t2= ilwis.do('selection',ts_t,"rasterbands(1)")
t3= ilwis.do('selection',ts_t,"rasterbands(2)")
t4= ilwis.do('selection',ts_t,"rasterbands(3)")
t5= ilwis.do('selection',ts_t,"rasterbands(4)")
t6= ilwis.do('selection',ts_t,"rasterbands(5)")
t7= ilwis.do('selection',ts_t,"rasterbands(6)")
t8= ilwis.do('selection',ts_t,"rasterbands(7)")
t9= ilwis.do('selection',ts_t,"rasterbands(8)")
t10= ilwis.do('selection',ts_t,"rasterbands(9)")
t11= ilwis.do('selection',ts_t,"rasterbands(10)")
t12= ilwis.do('selection',ts_t,"rasterbands(11)")
t13= ilwis.do('selection',ts_t,"rasterbands(12)")

Creating a Sorghum crop mask using the individual Transpiration layers selected

mask_season = ilwis.do('mapcalc', 'iff((@1>250)and(@1<400),1,0)',T_Season)# using season mask to differentiate from cotton which is having overall 
#seasonal higher Transpiration rates, and from groundnut, pigeonpea or onion which are having lower Transpiration rates

mask_th = ilwis.do('mapcalc', 'iff((@1>25)and(@1<70),1,0)',t7) #threshold on maximum mid-season T, to differentiate from cotton which is having higher 
#Transpiration rates during mid-season, and from groundnut, pigeonpea or onion which are having lower mid-season Transpiration rates

mask_up = ilwis.do('mapcalc', 'iff((@1<@2)and(@3<@4),1,0)',t2,t3,t6,t7) #increase in T from initial through development stage

mask_down = ilwis.do('mapcalc', 'iff((@2<@1),1,0)',t10,t11) #decrease in T in late season stage

mask_end = ilwis.do('mapcalc', 'iff(@1<20,1,0)',t13) #hardly any T during late - end of season stage

#extract the final sorghum mask - should meet all above consitions
mask = ilwis.do('mapcalc', 'iff(@1+@2+@3+@4+@5==5,1,0)',mask_season,mask_th,mask_up,mask_down,mask_end)
mask.store('sorghum_mask_'+year+'.mpr')

Calculate land productivity (LP)

Following crop parameters for Sorghum are used:

• Harvest Index (HI) = 0.25
• Above ground over total biomass (AOT) = 0.8 # above ground over total biomass production ratio
• Moisture Content (MC) = 0.2 # moisture content, dry matter over freshbiomass
• Light use efficiency correction factor (fc) = 1.6 # Light use efficiency correction factor for C4 crops

AGBM = (AOT * fc * (NPP * 22.222 / (1 - MC))) / 1000 # factor of 1000 is to convert from kg to ton

where AGMB is the above ground biomass, AOT is the above ground over total biomass production fraction, fc is the light use efficiency correction factor, NPP is seasonal net primary production in gC/m²/season, MC is moisture content in the fresh biomass

to calculate crop yield: Crop yield = HI*AGBM

where HI is harvest index, and AGBM is above ground biomass

AGBM = ilwis.do('mapcalc', '(0.8*1.6*(@1*22.222/(1-0.2)))/1000',NPP_Season) #convert from kg to ton
AGBM_sorghum = ilwis.do('mapcalc', 'iff(@1=1,@2,?)',mask, AGBM) 

Yield = ilwis.do('mapcalc', '(0.25*@1)',AGBM) 
Yield_sorghum = ilwis.do('mapcalc', 'iff(@1=1,@2,?)',mask, Yield) 

AGBM_sorghum.store('agbm_sorghum.mpr') #tons/ha/season
Yield_sorghum.store('yield_sorghum_'+year+'.mpr') #tons/ha/season

#Get statistics and plot graph
def descriptive_statistics(raster_input):
    stats_input = raster_input.statistics(ilwis.PropertySets.pHISTOGRAM)
    print()
    print('Minimum: ',stats_input[ilwis.PropertySets.pMIN]) 
    print('Maximum: ',stats_input[ilwis.PropertySets.pMAX]) 
    print('Mean: ',stats_input[ilwis.PropertySets.pMEAN]) 
    print()
    x=[a for (a,b) in stats_input.histogram()][:-1] 
    y=[b for (a,b) in stats_input.histogram()][:-1]
    plt.figure(figsize=(4, 4))
    plt.plot(x,y,label='Map values')
    plt.xlabel('Data Range')
    plt.ylabel('Data Frequency')
    plt.title('Map Histogram')
    plt.legend()

descriptive_statistics(AGBM_sorghum)

Minimum:  0.0
Maximum:  6.136969740800002
Mean:  3.5321959069872015

descriptive_statistics(Yield_sorghum)

Minimum:  0.0
Maximum:  1.5342424352000006
Mean:  0.8830489767468004

Calculate gross and net water productivity (WP)

Calculate gross , net biomass and crop yield water productivity

To calculate Gross Biomass Water Productivity: WPg_biomass = Biomass / AETI *100

where Biomass = seasonal AGBM and AETI = seasonal AETI

To calculate Net Biomass Water Productivity: WPn_biomass = Biomass / T *100

where Biomass = seasonal AGBM and T = seasonal Transpitation #fraction beneficially consumed by the crop

WPg_biomass = ilwis.do('mapcalc', '(@1/@2*100)',AGBM, AETI_Season) # for kg/m3
WPg_biomass = ilwis.do('mapcalc', 'iff(@1=1,@2,?)',mask, WPg_biomass) 

WPg_biomass.store('WPg_biomass_sorghum.mpr')

descriptive_statistics(WPg_biomass)

Minimum:  0.0
Maximum:  1.3081625578006872
Mean:  0.8797796511394146

WPn_biomass = ilwis.do('mapcalc', '(@1/@2*100)',AGBM, T_Season) # for kg/m3
WPn_biomass = ilwis.do('mapcalc', 'iff(@1=1,@2,?)',mask, WPn_biomass) 

WPn_biomass.store('WPn_biomass_sorghum_'+year+'.mpr')

descriptive_statistics(WPn_biomass)

Minimum:  0.0
Maximum:  1.7165088163609028
Mean:  1.0680250834134348

To calculate Crop Yield Water Productivity: Yield WP = Yield / AETI *100
where Biomass = seasonal Yield and AETI = seasonal AETI

WPyield = ilwis.do('mapcalc', '(@1/@2*100)',Yield, AETI_Season) # for kg/m3
WPyield = ilwis.do('mapcalc', 'iff(@1=1,@2,?)',mask, WPyield) 

WPyield.store('WPyield_sorghum_'+year+'.mpr')

descriptive_statistics(WPyield)

Minimum:  0.0
Maximum:  0.3270406394501718
Mean:  0.21994491278485365

Calculate some other performance indicators

Uniformity of water consumption

Equity is defined as the coefficients of variation (CV) of seasonal AETI in the area of interest. It measures the evenness of the water supply in an irrigation scheme. It is defined as: CV_AETI = (AETIsd / AETIm) * 100

where:

• AETIsd = standard deviation
• AETIm = mean

CV classification, according to Bastiaanssen et al., 1996):

• CV of 0 to 10% is good
• CV of 10 to 25% is fair
• CV > 25% is poor uniformity

#transform ilwis raster to numpy array 
AETI_np = np.fromiter(iter(AETI_Season), np.float64, AETI_Season.size().linearSize()) #transform to 1D array
AETI_np = AETI_np.reshape((AETI_Season.size().ysize, AETI_Season.size().xsize))

#assign background value of -130 to np.nan
AETI_np[AETI_np == -130] = np.nan 

pmean_a = np.nanmean(AETI_np)
pstdev_a = np.nanstd(AETI_np)

print('mean =', pmean_a) # mean value
print('std =', pstdev_a) #std value

mean = 295.83773173010417
std = 176.36385829131305

CV_AETI_wholescheme = (pstdev_a/pmean_a) * 100
print('CV of the whole scheme (%) = ', CV_AETI_wholescheme)

CV of the whole scheme (%) =  59.61506575240092

#set background to -130 and transform ilwis raster to numpy array for sorghum only
CV_AETI_sorghum_2np = ilwis.do('mapcalc', 'iff(@1==1,@2,-130)', mask, AETI_Season) 
#CV_AETI_sorghum_2np.store('CV_AETI_sorghum_2np.mpr')

AETI_np = np.fromiter(iter(CV_AETI_sorghum_2np), np.float64, CV_AETI_sorghum_2np.size().linearSize()) #transform to 1D array
AETI_np = AETI_np.reshape((CV_AETI_sorghum_2np.size().ysize, CV_AETI_sorghum_2np.size().xsize))

#assign background value of -130 to np.nan
AETI_np[AETI_np == -130] = np.nan 

pmean_s = np.nanmean(AETI_np)
pstdev_s = np.nanstd(AETI_np)

print('mean =', pmean_s) # mean value
print('std =', pstdev_s) #std value

mean = 398.8134665368804
std = 45.566986769397786

CV_AETI_sorghum = (pstdev_s/pmean_s) * 100
print('CV of sorghum in the scheme (%) = ',CV_AETI_sorghum)

CV of sorghum in the scheme (%) =  11.425638949727883

Calculate Beneficial Fraction (BF)

Beneficial fraction is the ratio of the water that is consumed as transpiration (T) compared to overall field water consumption (AETI).
It is a measure of the efficiency of on farm water and agronomic practices in use of water for crop growth.

It is defined as: BF = T/AETI

BF_scheme = ilwis.do('mapcalc', 'iff(@1>=0,@1/@2,?)',T_Season, AETI_Season)
BF_scheme.store('BF_scheme.mpr')

BF_sorghum = ilwis.do('mapcalc', 'iff(@1==1,@2/@3,?)',mask,T_Season, AETI_Season)
BF_sorghum.store('BF_sorghum.mpr')

descriptive_statistics(BF_sorghum)

Minimum:  0.39272197962154287
Maximum:  0.904914004914005
Mean:  0.8261628844977315

Calculate Relative Water Deficit (RWD)

Relative water deficit = ETa / ETx

ETx = 99 percentile of the actual evapotranspiration

#transform from ilwis raster to numpy array whole scheme
AETI_np = np.fromiter(iter(AETI_Season), np.float64, AETI_Season.size().linearSize()) #transform to 1D array

#assign -130 to np.nan
AETI_np[AETI_np==-130]=np.nan

ETx = np.nanpercentile(AETI_np, 99)
AETI_mean = pmean_a
RWD = 1-(AETI_mean/ETx)

print('ETx = ',ETx)
print('AETI_mean = ',AETI_mean)
print('RWD = ',RWD)

ETx =  665.4000000000001
AETI_mean =  295.83773173010417
RWD =  0.5553986598585752

#transform from ilwis raster to numpy array sorghum only
AETI_np = np.fromiter(iter(CV_AETI_sorghum_2np), np.float64, CV_AETI_sorghum_2np.size().linearSize()) #transform to 1D array

#assign -130 to np.nan
AETI_np[AETI_np==-130]=np.nan

ETx = np.nanpercentile(AETI_np, 99)
AETI_mean = pmean_s
RWD = 1-(AETI_mean/ETx)

print('ETx = ',ETx)
print('AETI_mean = ',AETI_mean)
print('RWD = ',RWD)

ETx =  482.3
AETI_mean =  398.8134665368804
RWD =  0.17310083653974617

Eventually process another year and compare it with the results obtained for 2020