In [1]:
#load library
import os
import ilwis
from osgeo import gdal, osr
import numpy as np
import matplotlib.pyplot as plt
import pandas as pd
from datetime import datetime, timedelta
import warnings
warnings.filterwarnings('ignore')
gdal connector: ilwis3 connector: Organizing data: 

Retrieve data from the Copernicus climate data store¶

install the cdsapi: The CDS API client is a python based library. It provides support for both Python 2.7.x and Python 3.

login to the climate data store, register and if you click on the name you get the details of your user profile. Copy and paste the UID and API Key to a text file and change it according to the two line example below:

    url: https://cds.climate.copernicus.eu/api/v2
    key: 'Your_UID':'Your_API-KEY'

Upload the file created onto the root of your /surface_flood folder and use as filename: cdsapirc.txt

Execute the code field below. This will transfer the 'cdsapirc.txt' file to the root of the sytem (here /home/eoafrica) and will rename the file to '.cdsapirc'. You can check your result using in terminal. Navigate to /home/eoafrica and type the followning command: ls -a Executing this command will also show the hidden files.

In [ ]:
#handle credentials (using user defined Copernicus Climate Data Store credentials)
#create a text file with filename 'cdsapirc.txt' in the current folder and add the 2 lines as indicated above

from pathlib import Path
print(Path.home())
!pwd
!mv cdsapirc.txt .cdsapirc
!cp .cdsapirc /home/eoafrica #check if the file .cdsapirc is appearing in the listing
!ls -a /home/eoafrica
!mv .cdsapirc cdsapirc.txt


#check if the file .cdsapirc is appearing in the listing
#!ls -a /home/eoafrica

You can Install the CDS API client via the package management system pip, type the following command promt window in the python folder: !pip install cdsapi - see code field below

Now you are ready to start retrieving data from the CDS

In [ ]:
!pip install cdsapi

The gridded river discharge data used below will be downloaded for a selected area, the Buzi catchment in Mozambique, for the years 2021 and 2022. The data represents the discharge in m3/sec (per pixel) generated by the GloFAS system (see: https://www.globalfloods.eu/). The data will be retrieved in grib2 format and will be dowloaded from: https://cds.climate.copernicus.eu/cdsapp#!/dataset/cems-glofas-historical?tab=form. The data is going to be assigned as an ILWISPy raster coverage in the cells below, the dimensions are checked and a plot of the merged dataset is created.

Link to data portal https://cds.climate.copernicus.eu/cdsapp#!/dataset/cems-glofas-historical?tab=overview

Limitation of data download is that only up to 500 layers (days) can be retrieved given a download request, so here a yearly data download of the area of interest is used

In [2]:
#specify your output folder
data_dir = r'./glofas_data/'
In [ ]:
# Import of version 4 glofas river discharge Consolidated!

import cdsapi
year = '2022'
c = cdsapi.Client()

c.retrieve(
    'cems-glofas-historical',
    {
        'system_version': 'version_4_0',
        'variable': 'river_discharge_in_the_last_24_hours',
        'format': 'grib',
        'hydrological_model': 'lisflood',
        'product_type': 'consolidated',
        'hyear': (year),
        'hday': [
            '01', '02', '03',
            '04', '05', '06',
            '07', '08', '09',
            '10', '11', '12',
            '13', '14', '15',
            '16', '17', '18',
            '19', '20', '21',
            '22', '23', '24',
            '25', '26', '27',
            '28', '29', '30',
            '31',
        ],
        'hmonth': [
            'april', 'august', 'december',
            'february', 'january', 'july',
            'june', 'march', 'may',
            'november', 'october', 'september',
        ],
        'area': [
            -18.5, 32.25, -21.0,
            35.0,
        ],
    },
    (data_dir) +'/' +'historical_'+(year)+'.grib')

Processing the GLOFAS Historical Discharge data, note that GloFAS version 4 is used here, resolution is 0.05 degree, see 'system_version' above

In [3]:
work_dir = os.getcwd() + '/glofas_data'

#set the working directory for ILWISPy
ilwis.setWorkingCatalog(work_dir)
print(work_dir)
ilwis.version()
/home/eoafrica/surface_flood/glofas_data
gdal connector: ilwis3 connector: Organizing data: 

Out[3]:
'1.0 build 20230929'
In [6]:
# Load retrieved data - ignore grib Warning
glofas_2021 = ilwis.RasterCoverage('historical_2021.grib')
In [7]:
# Load retrieved data - ignore grib Warning 
glofas_2022 = ilwis.RasterCoverage('historical_2022.grib')
In [8]:
print(glofas_2021.size())
print(glofas_2022.size())
Size(57, 52, 365)
Size(57, 52, 365)

Check some of the other meta data

In [9]:
coordSys = glofas_2021.coordinateSystem()
print(coordSys.toWKT())
print()
print(coordSys.toProj4())
print()
print(glofas_2021.envelope())
print()
stats2021 = glofas_2021.statistics(ilwis.PropertySets.pHISTOGRAM)
#print(stats.histogram())
print(stats2021[ilwis.PropertySets.pMIN]) # minimum value on the map
print(stats2021[ilwis.PropertySets.pMAX]) # maximum value on the map
print()
stats2022 = glofas_2022.statistics(ilwis.PropertySets.pHISTOGRAM)
#print(stats.histogram())
print(stats2022[ilwis.PropertySets.pMIN]) # minimum value on the map
print(stats2022[ilwis.PropertySets.pMAX]) # maximum value on the map
PROJCS["historical_2021.grib",GEOCS["historical_2021.grib",DATUM[" unnamed",[?],ELLIPSOID["Normal Sphere (r=6370997)",6370997.000000000000,0.000000000000],PRIMEM["Greenwich",0, AUTHORITY["EPSG",8901"]]],PROJECTION[""],,UNIT[degrees,1.0]]

+proj=longlat +a=6370997.000000000000 +b=6370997.000000000000

32.200000 -21.050000 35.050000 -18.450000

0.0
7799.625

0.0
2311.28125

Transform the stack (in ilwis format) to a 1D numpy array (using np.fromiter(iter()) for the two years under consideration, combine the two 1D arrays (using append) and transform back to a 3D array (using reshape)

In [10]:
#create 1D arrays
glofas_2np_2021 = np.fromiter(iter(glofas_2021), np.float64, glofas_2021.size().linearSize())
glofas_2np_2022 = np.fromiter(iter(glofas_2022), np.float64, glofas_2022.size().linearSize())

#combine the two - using append
array_both = np.append(glofas_2np_2021, glofas_2np_2022)

#create a 3D stack of the tow years combined
glofas_2np = array_both.reshape((glofas_2021.size().zsize+glofas_2022.size().zsize, glofas_2021.size().ysize, glofas_2022.size().xsize))
In [11]:
#check the dimensions obtained
print(glofas_2np.shape)
(730, 52, 57)

Create a simple interactive plot using MatPlotLib Imshow of the first layer

In [13]:
# create a plot of the numpy array using MatPlotLib - Imshow
fig1 = plt.figure(figsize =(10, 7))
plt.imshow(glofas_2np[0], interpolation='none', vmin=0, vmax=500, cmap='jet')

plt.axis('on')
plt.colorbar(shrink=0.6)
plt.title('Discharge (m3-sec)');

In the cell below the data for the outlet location of the Buzi River is retrieved for 2020, this is column, row location X=50 and y=28

In [14]:
l = 28
c = 50 
z = (glofas_2021.size().zsize)
print(z)

Q_2021 = []
for n in range(0,(z)):
    point_value = (glofas_2021.pix2value(ilwis.Pixel(c,l,n)))
    Q_2021.append(point_value)

print('Values extracted for selected location:', Q_2021)
365
Values extracted for selected location: [791.15625, 780.796875, 739.453125, 711.84375, 680.453125, 662.1875, 659.875, 686.21875, 730.328125, 795.375, 836.4375, 833.859375, 814.40625, 809.59375, 840.25, 896.75, 951.71875, 1277.3125, 2162.0, 2325.15625, 1961.109375, 1730.96875, 2065.4375, 3331.03125, 2911.03125, 2335.71875, 2016.5625, 1814.15625, 1693.859375, 1624.109375, 1531.90625, 1427.3125, 1366.5, 1436.625, 1428.78125, 1335.34375, 1240.3125, 1163.0625, 1107.4375, 1076.3125, 1026.90625, 987.96875, 962.3125, 930.84375, 999.40625, 1170.0, 1684.03125, 3062.3125, 3439.125, 2875.4375, 2359.5625, 2030.03125, 1798.8125, 1652.3125, 1555.71875, 1455.5, 1357.75, 1275.40625, 1261.375, 1270.1875, 1217.78125, 1141.28125, 1067.1875, 1008.03125, 973.0625, 931.125, 889.15625, 847.6875, 803.40625, 758.59375, 717.9375, 680.15625, 647.96875, 617.375, 587.84375, 558.5625, 532.125, 509.6875, 487.28125, 466.75, 447.4375, 431.25, 434.5, 455.3125, 457.21875, 441.28125, 435.21875, 430.5625, 416.625, 399.5625, 381.3125, 363.40625, 346.5625, 330.84375, 316.40625, 303.3125, 291.65625, 282.5625, 274.90625, 265.28125, 255.15625, 245.96875, 237.40625, 229.0625, 220.90625, 213.0625, 205.96875, 203.21875, 206.25, 204.78125, 198.40625, 191.28125, 184.78125, 179.28125, 175.3125, 171.5625, 166.40625, 160.375, 154.4375, 149.28125, 145.875, 146.03125, 154.75, 170.15625, 175.09375, 169.34375, 161.34375, 153.5625, 146.25, 139.5, 133.3125, 127.8125, 122.75, 118.03125, 113.3125, 108.71875, 104.40625, 100.34375, 96.6875, 93.65625, 91.3125, 89.25, 87.03125, 85.0625, 83.9375, 83.6875, 82.90625, 80.84375, 77.96875, 74.9375, 72.03125, 69.34375, 67.0, 65.21875, 65.4375, 68.71875, 72.75, 73.875, 72.59375, 71.1875, 71.09375, 71.65625, 70.75, 68.65625, 66.21875, 63.65625, 60.90625, 58.125, 55.4375, 52.9375, 50.78125, 49.28125, 49.40625, 51.78125, 54.5625, 55.3125, 54.21875, 52.34375, 50.234375, 48.03125, 45.8125, 43.640625, 41.609375, 39.703125, 37.9375, 36.265625, 34.703125, 33.28125, 32.09375, 31.34375, 31.578125, 32.40625, 32.640625, 31.9375, 30.75, 29.609375, 29.09375, 29.703125, 32.703125, 38.703125, 42.953125, 43.234375, 41.65625, 39.703125, 37.734375, 35.828125, 34.109375, 32.71875, 31.703125, 31.09375, 30.859375, 30.609375, 30.015625, 29.171875, 28.203125, 27.171875, 26.171875, 25.234375, 24.375, 23.5625, 22.75, 21.96875, 21.265625, 20.6875, 20.28125, 20.40625, 21.765625, 25.5, 38.546875, 60.171875, 68.5625, 66.421875, 62.0, 57.671875, 53.796875, 50.328125, 47.515625, 45.40625, 43.421875, 41.1953125, 38.890625, 36.796875, 35.0546875, 34.3125, 35.6328125, 37.3515625, 37.34375, 35.9921875, 34.2578125, 32.515625, 30.859375, 29.34375, 28.0546875, 27.125, 26.90625, 28.2890625, 31.9375, 35.234375, 35.7421875, 34.421875, 32.7265625, 31.3046875, 30.65625, 32.046875, 35.1875, 37.1796875, 37.0, 36.3046875, 36.7265625, 37.296875, 36.59375, 34.9609375, 33.1015625, 31.5, 30.4765625, 29.9453125, 29.703125, 29.4921875, 29.3125, 30.171875, 33.1796875, 37.3828125, 41.078125, 43.6640625, 44.5, 45.390625, 48.140625, 49.7421875, 48.5625, 46.0703125, 43.7265625, 42.3828125, 42.4453125, 45.1640625, 52.8046875, 61.1796875, 63.234375, 60.921875, 57.9140625, 56.6015625, 57.25, 57.46875, 57.0625, 57.3359375, 57.3359375, 55.75, 53.078125, 50.2421875, 47.5703125, 45.1328125, 42.984375, 41.09375, 39.2421875, 37.4140625, 35.734375, 34.3203125, 33.453125, 31.5859375, 29.9921875, 29.453125, 28.9609375, 27.96875, 26.6640625, 25.2890625, 24.1953125, 24.1015625, 27.7421875, 40.984375, 57.9140625, 68.6640625, 76.96875, 85.59375, 90.5625, 89.890625, 85.3125, 79.375, 73.5, 69.078125, 65.15625, 61.0625, 57.25, 53.828125, 50.84375, 48.203125, 45.75, 43.34375, 41.0, 38.828125, 36.84375, 35.21875, 34.8125, 41.0, 65.6875, 115.15625, 146.75, 149.359375, 147.609375, 155.546875, 178.265625, 200.390625, 212.84375, 222.921875, 223.5625, 211.75, 200.609375]

Repeat the procedure for 2022

In [15]:
l = 28
c = 50 
z = (glofas_2021.size().zsize)
print(z)

Q_2022 = []
for n in range(0,(z)):
    point_value = (glofas_2022.pix2value(ilwis.Pixel(c,l,n)))
    Q_2022.append(point_value)

print('Values extracted for selected location:', Q_2022)
365
Values extracted for selected location: [194.78125, 187.890625, 180.328125, 177.90625, 180.296875, 179.4375, 183.5, 200.4375, 209.875, 207.59375, 201.125, 190.3125, 180.125, 173.0625, 170.90625, 172.671875, 172.46875, 166.90625, 159.53125, 152.53125, 146.09375, 139.5, 134.1875, 131.390625, 152.015625, 433.0, 758.890625, 773.015625, 687.515625, 620.828125, 563.8125, 512.421875, 465.734375, 425.796875, 391.234375, 363.703125, 359.34375, 367.203125, 364.203125, 348.015625, 326.09375, 304.46875, 284.09375, 265.765625, 249.203125, 233.703125, 219.71875, 207.90625, 198.765625, 193.25, 192.5, 191.375, 186.65625, 179.0625, 170.84375, 166.4375, 168.8125, 174.5625, 174.9375, 169.53125, 162.8125, 155.59375, 148.09375, 141.375, 134.96875, 128.90625, 125.90625, 128.625, 131.84375, 131.6875, 129.59375, 125.40625, 120.0, 114.40625, 109.09375, 107.4375, 156.625, 272.125, 342.40625, 341.9375, 313.125, 284.53125, 260.34375, 240.5, 223.59375, 209.65625, 220.96875, 261.84375, 275.21875, 263.5, 245.53125, 228.40625, 213.9375, 201.53125, 190.6875, 183.9375, 191.03125, 211.28125, 229.6875, 299.1875, 461.375, 569.3125, 550.46875, 505.0, 469.75, 437.875, 413.03125, 402.53125, 393.46875, 380.09375, 376.65625, 409.84375, 442.1875, 421.625, 386.65625, 355.1875, 329.34375, 312.875, 300.9375, 286.15625, 269.40625, 253.21875, 239.4375, 228.53125, 219.15625, 210.53125, 202.25, 193.78125, 184.375, 175.03125, 166.65625, 159.5, 153.0625, 146.59375, 140.0, 133.6875, 127.96875, 123.28125, 119.9375, 116.96875, 114.1875, 111.53125, 112.5, 134.6875, 173.5, 181.875, 173.375, 163.34375, 154.21875, 146.34375, 139.375, 132.59375, 127.0625, 123.5625, 120.09375, 116.4375, 114.6875, 114.5, 114.625, 117.03125, 116.6875, 112.3125, 107.15625, 103.34375, 102.59375, 101.40625, 97.625, 92.96875, 88.5, 84.40625, 80.71875, 77.28125, 73.96875, 71.125, 69.0, 67.53125, 66.671875, 66.203125, 65.984375, 69.265625, 79.6875, 96.90625, 110.171875, 115.0, 114.484375, 110.671875, 105.109375, 98.96875, 93.0, 87.4375, 82.296875, 77.546875, 73.1875, 69.234375, 65.90625, 63.40625, 61.28125, 58.921875, 56.5625, 54.640625, 53.625, 53.546875, 53.484375, 53.34375, 53.90625, 54.84375, 54.703125, 53.453125, 51.859375, 50.40625, 48.765625, 46.6875, 44.4375, 42.234375, 40.265625, 38.671875, 38.078125, 38.90625, 39.671875, 39.09375, 37.640625, 35.96875, 34.375, 32.921875, 31.546875, 30.25, 29.03125, 27.890625, 26.859375, 25.921875, 25.1875, 24.6875, 24.34375, 23.96875, 23.53125, 23.171875, 23.5, 24.90625, 26.3125, 26.7421875, 26.3203125, 25.5234375, 24.625, 23.8125, 23.3359375, 24.40625, 28.203125, 32.0078125, 33.1328125, 32.4140625, 31.0703125, 29.6796875, 28.421875, 27.3046875, 26.25, 25.171875, 24.1171875, 23.2109375, 22.640625, 23.4765625, 27.6171875, 32.65625, 34.328125, 33.390625, 31.6953125, 30.0, 28.4375, 27.078125, 25.8984375, 24.8203125, 23.8046875, 22.90625, 22.140625, 21.4375, 20.8046875, 20.25, 19.703125, 19.1328125, 18.5546875, 17.984375, 17.4609375, 16.984375, 16.578125, 16.234375, 15.921875, 15.6484375, 15.375, 15.1171875, 14.875, 14.703125, 14.5859375, 14.5234375, 14.4921875, 14.484375, 14.4921875, 14.484375, 14.4609375, 14.421875, 14.421875, 14.9140625, 16.890625, 20.6171875, 24.2421875, 26.1015625, 26.8125, 28.1796875, 32.9140625, 56.5546875, 105.4765625, 132.484375, 134.734375, 130.984375, 127.8515625, 125.7734375, 131.2421875, 147.96875, 172.53125, 180.265625, 181.6484375, 208.5, 229.859375, 220.0078125, 202.453125, 186.2109375, 172.0859375, 160.7109375, 151.8359375, 145.3359375, 137.8984375, 129.9375, 123.671875, 130.1796875, 336.2421875, 561.8828125, 539.9140625, 477.71875, 418.53125, 385.7421875, 357.46875, 327.3984375, 301.640625, 288.96875, 299.453125, 315.109375, 312.4296875, 296.5, 279.03125, 262.21875, 247.109375, 237.7109375, 244.5625, 271.296875, 308.375, 371.3515625, 444.859375, 499.3984375, 538.1875, 535.640625, 500.265625, 458.7734375, 421.7890625, 389.265625, 360.5078125, 335.265625, 313.1875]

The two lists obtained above are now merged into a single list

In [16]:
#concatonate 2 lists
Q_all = Q_2021 + Q_2022
print(len(Q_all))
730
In [17]:
# initializing date range for X-Axis
plot_date = datetime.strptime("01-1-2021", "%d-%m-%Y")

z = (glofas_2021.size().zsize+glofas_2022.size().zsize)

date_generated = pd.date_range(plot_date, periods=z)
print(date_generated.strftime("%d-%m-%Y"))
Index(['01-01-2021', '02-01-2021', '03-01-2021', '04-01-2021', '05-01-2021',
       '06-01-2021', '07-01-2021', '08-01-2021', '09-01-2021', '10-01-2021',
       ...
       '22-12-2022', '23-12-2022', '24-12-2022', '25-12-2022', '26-12-2022',
       '27-12-2022', '28-12-2022', '29-12-2022', '30-12-2022', '31-12-2022'],
      dtype='object', length=730)
In [18]:
fig2 = plt.figure(figsize =(14, 7))
plt.plot(date_generated, Q_all, label='Discharge')
plt.xlabel('Time')
plt.ylabel('Discharge (m3/sec)')
plt.title('Discharge extracted for a given pixel')
plt.legend();

Exercise:¶

Create a new plot, including also the data from 2019 and 2020, note that 2020 is a leap year!

In [ ]: