Import GLEAM E data

• E - Actual evaporation [mm/day]

• by definition: E = Et + Eb + Ew + Ei + Es + Ec

#import required libraries
import os
import ilwis
import netCDF4 as nc
import numpy as np
import numpy.ma as ma
#from osgeo import gdal, osr
import matplotlib.pyplot as plt

data_dir = os.getcwd() + '/data/'

work_dir = os.getcwd() + '/results'

ilwis.setWorkingCatalog(work_dir)
print(work_dir)

d:\jupyter\notebook_scripts\Special\GLEAM/results

Provide the location of your input data file, note that the folder specifications provided below expects that the data is within a specific product folder within the notebook folder. Here 2010 is processed

fn = data_dir+'/E_2008_GLEAM_v4.3a.nc'
ds = nc.Dataset(fn)

print(ds)

<class 'netCDF4.Dataset'>
root group (NETCDF4 data model, file format HDF5):
    Dataset: Global Land Evaporation Amsterdam Model
    Version: 4.3a
    Authors: Hydro-Climate Extremes Lab (H-CEL)
    Institution: Ghent University
    Contact: info@gleam.eu
    Reference1: Miralles, D.G. et al. 2025: global land evaporation and soil moisture dataset at 0.1° resolution from 1980 to near present, Scientific Data, 12, 416, doi: 10.1038/s41597-025-04610-y
    Reference2: Miralles, D.G. et al. 2011: Global land-surface evaporation estimated from satellite-based observations, Hydrology and Earth System Sciences, 15, 453-469, doi: 10.5194/hess-15-453-2011
    dimensions(sizes): time(366), lat(1800), lon(3600)
    variables(dimensions): float32 E(time, lat, lon), float64 lon(lon), int64 time(time), float64 lat(lat)
    groups: 

for dim in ds.dimensions.values():
    print(dim)

"<class 'netCDF4.Dimension'>": name = 'time', size = 366
"<class 'netCDF4.Dimension'>": name = 'lat', size = 1800
"<class 'netCDF4.Dimension'>": name = 'lon', size = 3600

for var in ds.variables.values():
    print(var)

<class 'netCDF4.Variable'>
float32 E(time, lat, lon)
    _FillValue: -999.0
    standard_name: Actual evaporation
    long_name: Actual evaporation from GLEAM 4.3a
    units: mm.day-1
unlimited dimensions: 
current shape = (366, 1800, 3600)
filling on
<class 'netCDF4.Variable'>
float64 lon(lon)
    _FillValue: nan
    standard_name: longitude
    long_name: longitude
    units: degrees_east
unlimited dimensions: 
current shape = (3600,)
filling on
<class 'netCDF4.Variable'>
int64 time(time)
    standard_name: time
    long_name: time
    units: days since 1900-01-01
    calendar: proleptic_gregorian
unlimited dimensions: 
current shape = (366,)
filling off
<class 'netCDF4.Variable'>
float64 lat(lat)
    _FillValue: nan
    standard_name: latitude
    long_name: latitude
    units: degrees_north
unlimited dimensions: 
current shape = (1800,)
filling on

et = ds.variables["E"][:] 
et.shape

(366, 1800, 3600)

print(type(et))

<class 'numpy.ma.MaskedArray'>

et_day1 = et[0]

plt.figure(figsize=(12, 6))

im = plt.imshow(
    et_day1,
    cmap="viridis",
    vmin = 0, 
    vmax = 10
)

plt.colorbar(im, label="ET (mm/day)",  shrink = 0.5, extend = 'max')
plt.title("Evapotranspiration - Day 1")
plt.xlabel("Longitude index")
plt.ylabel("Latitude index")

plt.show()

grf_LL= ilwis.GeoReference('code=georef:type=corners, csy=epsg:4326, envelope=-180 90 180 -90 , gridsize=3600 1800, cornerofcorners=yes')

dfNum = ilwis.DataDefinition(ilwis.NumericDomain('code=value'), ilwis.NumericRange(0.0, 10000.0, 0))
GLEAMNew = ilwis.RasterCoverage()
GLEAMNew.setSize(ilwis.Size(3600, 1800,366))
GLEAMNew.setGeoReference(grf_LL)
GLEAMNew.setDataDef(dfNum)

# 1 D array
E = np.array([et]).flatten()
print(E.shape)
print(E)

(2371680000,)
[-9.9900000e+02 -9.9900000e+02 -9.9900000e+02 ...  1.3336545e-01
  1.3336235e-01 -9.9900000e+02]

GLEAMNew.array2raster(E)
print(GLEAMNew.size())

Size(3600, 1800, 366)

grf_LLsub = ilwis.GeoReference('code=georef:type=corners, csy=epsg:4326, envelope=105 -4.58333333 115.833333 -10 , gridsize=1300 650, cornerofcorners=yes')

et_res = ilwis.do('resample', GLEAMNew, grf_LLsub, 'nearestneighbour')
print(et_res.size())

Size(1300, 650, 366)

year = os.path.basename(fn).split("_")[1]
print(year)

2008

et_res1 = ilwis.do('mapcalc', 'iff(@1==-999, ?, @1)',et_res) 
et_res1.store('Gleam_E_sub_'+year+'.mpl')

select = ilwis.do('selection',et_res1,"rasterbands(0)") 

#Just for visualization, as data is already stored in ILWIS format and map can be displayed using ILWIS386
ETsub_2np = np.fromiter(iter(select), np.float64, select.size().linearSize()) 
ETsub_2np = ETsub_2np.reshape((select.size().ysize, select.size().xsize))

#display the numpy array created in the previous step using matplotlib 
fig = plt.figure(figsize =(12, 7))

maximum = 8
plt.imshow(ETsub_2np, interpolation='none', vmin=0, vmax=maximum, cmap='jet')

plt.axis('on')

cbar = plt.colorbar(shrink=0.45, extend='max')
cbar.set_label("ETIa (mm/day)", fontsize=10)
plt.title('GLEAM E of Java');