OpenEO Platform logo

openEO Platform - Analysis Ready Data¶

Get started with openEO Platform: Basics ** Review article "The openEO API–Harmonising the Use of Earth Observation Cloud Services Using Virtual Data Cube Functionalities"** https://www.researchgate.net/publication/350109855_The_openEO_API-Harmonising_the_Use_of_Earth_Observation_Cloud_Services_Using_Virtual_Data_Cube_Functionalities

Review: https://hub.openeo.org/

Review: https://editor.openeo.cloud/

Review: https://openeo.org/documentation/1.0/datacubes.html

In [1]:
# import libraries
import os
import sys
import ilwis
import numpy as np
from netCDF4 import Dataset
from osgeo import gdal, osr
from osgeo.gdalconst import *
import matplotlib.pyplot as plt
gdal connector: ilwis3 connector: Organizing data: 

In [2]:
import openeo
openeo.client_version()
Out[2]:
'0.25.0'
In [3]:
!rm -rf /home/eoafrica/.local/share/openeo-python-client/refresh-tokens.json

Connect to openEO Platform and print list of available collections¶

In [4]:
con = openeo.connect(url="openeo.vito.be")
In [5]:
con.authenticate_oidc()
Visit https://aai.egi.eu/device?user_code=NBMA-YNSB 📋 to authenticate.
[###################################--] ✅ Authorized successfully
Authenticated using device code flow.
Out[5]:
<Connection to 'https://openeo.vito.be/openeo/1.2/' with OidcBearerAuth>
In [6]:
print(con.list_collection_ids())
['MAPEO_WATER_TUR_V1', 'COP_DEM_EU_25M', 'ESA_WORLDCEREAL_ACTIVECROPLAND', 'ESA_WORLDCEREAL_IRRIGATION', 'ESA_WORLDCEREAL_TEMPORARYCROPS', 'ESA_WORLDCEREAL_WINTERCEREALS', 'ESA_WORLDCEREAL_MAIZE', 'ESA_WORLDCEREAL_SPRINGCEREALS', 'SENTINEL1_GRD_SIGMA0', 'S1_GRD_SIGMA0_ASCENDING', 'S1_GRD_SIGMA0_DESCENDING', 'SENTINEL3_SYNERGY_VG1', 'SENTINEL3_SYNERGY_VG10', 'TERRASCOPE_S2_FAPAR_V2', 'TERRASCOPE_S2_NDVI_V2', 'TERRASCOPE_S2_LAI_V2', 'TERRASCOPE_S2_FCOVER_V2', 'TERRASCOPE_S2_TOC_V2', 'TERRASCOPE_S1_SLC_COHERENCE_V1', 'SENTINEL1_GAMMA0_SENTINELHUB', 'SENTINEL1_GRD', 'SENTINEL2_L1C_SENTINELHUB', 'SENTINEL2_L2A_SENTINELHUB', 'SENTINEL2_L2A_MOSAIC_120', 'PROBAV_L3_S10_TOC_333M', 'PROBAV_L3_S5_TOC_100M', 'PROBAV_L3_S1_TOC_100M', 'PROBAV_L3_S1_TOC_333M', 'TERRASCOPE_S5P_L3_CO_TD_V1', 'TERRASCOPE_S5P_L3_CO_TM_V1', 'TERRASCOPE_S5P_L3_CO_TY_V1', 'LANDSAT8-9_L1', 'LANDSAT8-9_L2', 'LANDSAT4-5_TM_L1', 'LANDSAT4-5_TM_L2', 'LANDSAT7_ETM_L1', 'LANDSAT7_ETM_L2', 'MODIS', 'SENTINEL3_OLCI_L1B', 'SENTINEL3_SLSTR', 'COPERNICUS_30', 'COPERNICUS_90', 'TERRASCOPE_S2_RHOW_V1', 'TERRASCOPE_S2_CHL_V1', 'TERRASCOPE_S2_SPM_V1', 'TERRASCOPE_S2_TUR_V1', 'ESA_WORLDCOVER_10M_2020_V1', 'ESA_WORLDCOVER_10M_2021_V2', 'CORINE_LAND_COVER', 'CORINE_LAND_COVER_ACCOUNTING_LAYERS', 'GHS_BUILT_S2', 'GLOBAL_LAND_COVER', 'GLOBAL_SURFACE_WATER', 'SEA_ICE_INDEX', 'SENTINEL_5P_L2', 'MAPZEN_DEM', 'LANDSAT1-5_MSS_L1', 'WATER_BODIES', 'SENTINEL1_CARD4L', 'VEGETATION_INDICES', 'SEASONAL_TRAJECTORIES', 'VEGETATION_PHENOLOGY_AND_PRODUCTIVITY_PARAMETERS_SEASON_1', 'VEGETATION_PHENOLOGY_AND_PRODUCTIVITY_PARAMETERS_SEASON_2', 'FRACTIONAL_SNOW_COVER', 'OPENEO_CROPTYPE_2021', 'EGM2008', 'CGLS_NDVI_LTS_V2_GLOBAL', 'CGLS_GDMP300_V1_GLOBAL', 'CGLS_GDMP_V2_GLOBAL', 'CGLS_SSM_V1_EUROPE', 'CGLS_FAPAR_V2_GLOBAL', 'CGLS_LAI_V2_GLOBAL', 'CGLS_FCOVER_V2_GLOBAL', 'CGLS_LAI300_V1_GLOBAL', 'CGLS_NDVI300_V1_GLOBAL', 'CGLS_FCOVER300_V1_GLOBAL', 'CGLS_FAPAR300_V1_GLOBAL', 'CGLS_NDVI300_V2_GLOBAL', 'CGLS_NDVI_V3_GLOBAL', 'CGLS_NDVI_V2_GLOBAL', 'CGLS_SWI_V1_EUROPE', 'TERRASCOPE_S1_GAMMA0_V1', 'SENTINEL2_L1C_INCD', 'SENTINEL2_L2A', 'TERRASCOPE_S5P_L3_NO2_TD_V1', 'TERRASCOPE_S5P_L3_NO2_TM_V1', 'TERRASCOPE_S5P_L3_NO2_TY_V1', 'AGERA5', 'PLANETSCOPE', 'POPULATION_DENSITY', 'SENTINEL2_L1C']

Get some more details about the first collection / first couple of collections

In [7]:
print(con.list_collections()[0])
{'description': 'Optical data collected with airborne drone platform, and processed with MapEO Water software v1 at VITO into turbidity. Turbidity indicates the relative opacity of the water column. It is an optical water property and a measure for the amount of light scattered by constituents in the water column. The higher the scattered light intensity, the higher the turbidity. Constituents that causes water to be turbid include clay, silt, very tiny inorganic and organic matter, algae, dissolved coloured organic compounds, plankton, and other microscopic organisms. Turbidity is expressed in Formazin Nephelometric Units (FNU, according to the ISO 7027 standard).', 'extent': {'spatial': {'bbox': [[-180.0, -84.0, 180.0, 84.0]]}, 'temporal': {'interval': [['2018-02-01T00:00:00Z', None]]}}, 'id': 'MAPEO_WATER_TUR_V1', 'keywords': ['Orthoimagery', 'Water quality', 'Turbidity', 'Airborne drone data', 'UAV', 'RPAS', 'Drone', 'MicaSense', 'DJI', 'cm resolution', 'super high resultion', 'VITO'], 'license': 'proprietary', 'links': [{'href': 'https://services.terrascope.be/collectioncatalogue/srv/eng/catalog.search#/metadata/urn:eop:VITO:MAPEO_WATER_TUR_V1', 'rel': 'alternate', 'title': 'Collection Catalogue Entry'}, {'href': 'https://services.terrascope.be/catalogue/catalogue/description.geojson?collection=urn:eop:VITO:MAPEO_WATER_TUR_V1', 'rel': 'alternate', 'title': 'OpenSearch entry point'}, {'href': 'https://openeo.vito.be/openeo/1.2/collections', 'rel': 'root'}, {'href': 'https://openeo.vito.be/openeo/1.2/collections', 'rel': 'parent'}, {'href': 'https://openeo.vito.be/openeo/1.2/collections/MAPEO_WATER_TUR_V1', 'rel': 'self'}], 'stac_extensions': ['https://stac-extensions.github.io/datacube/v2.2.0/schema.json', 'https://stac-extensions.github.io/eo/v1.1.0/schema.json'], 'stac_version': '0.9.0', 'title': 'MapEO Water Turbidity (TUR) Products - V1'}
In [ ]:
print(con.list_collections()[0:3])

Get detailed information about a selected collection¶

In [16]:
con.describe_collection('GLOBAL_SURFACE_WATER')#web listing
Out[16]:
In [8]:
collection_id = 'GLOBAL_SURFACE_WATER'
In [9]:
temporal_extent = ["2019-01-01", "2020-01-01"]
spatial_extent = {"west": 32.37, "south": -20.941, "east": 34.758, "north": -18.669}

Run code field below 2 times, first time you get a warning message. We did not specify a 'band', as we would like to retrieve all subdatasets of the global surface water coverage

In [11]:
cube = con.load_collection(
    collection_id, 
    temporal_extent=temporal_extent,
    spatial_extent=spatial_extent,
)

Verify if output folder exists, else create a new folder, here 'data'

In [12]:
work_dir = os.getcwd()+'/data'

print("current dir is: %s" % (os.getcwd()))
print("current working directory is:",work_dir) 

if os.path.isdir(work_dir):
    print("Folder exists")
else:
    print("Folder doesn't exists")
    os.mkdir(work_dir)
current dir is: /home/eoafrica/surface_flood
current working directory is: /home/eoafrica/surface_flood/data
Folder exists
In [13]:
%%time 
cube.download(work_dir + '/GSW.nc')
CPU times: user 41.8 ms, sys: 54.1 ms, total: 95.9 ms
Wall time: 1min 2s
In [17]:
input_file = (work_dir + '/GSW.nc')
print(input_file)
/home/eoafrica/surface_flood/data/GSW.nc
In [18]:
# read the meta data of the downloaded NC file - raster data
dataset = gdal.Open(input_file)
nc_info = gdal.Info(dataset)
print (nc_info)
Driver: netCDF/Network Common Data Format
Files: /home/eoafrica/surface_flood/data/GSW.nc
Size is 512, 512
Metadata:
  NC_GLOBAL#Conventions=CF-1.9
  NC_GLOBAL#institution=openEO platform
Subdatasets:
  SUBDATASET_1_NAME=NETCDF:"/home/eoafrica/surface_flood/data/GSW.nc":change
  SUBDATASET_1_DESC=[1x9088x9552] change (16-bit integer)
  SUBDATASET_2_NAME=NETCDF:"/home/eoafrica/surface_flood/data/GSW.nc":extent
  SUBDATASET_2_DESC=[1x9088x9552] extent (16-bit integer)
  SUBDATASET_3_NAME=NETCDF:"/home/eoafrica/surface_flood/data/GSW.nc":occurrence
  SUBDATASET_3_DESC=[1x9088x9552] occurrence (16-bit integer)
  SUBDATASET_4_NAME=NETCDF:"/home/eoafrica/surface_flood/data/GSW.nc":recurrence
  SUBDATASET_4_DESC=[1x9088x9552] recurrence (16-bit integer)
  SUBDATASET_5_NAME=NETCDF:"/home/eoafrica/surface_flood/data/GSW.nc":seasonality
  SUBDATASET_5_DESC=[1x9088x9552] seasonality (16-bit integer)
  SUBDATASET_6_NAME=NETCDF:"/home/eoafrica/surface_flood/data/GSW.nc":transitions
  SUBDATASET_6_DESC=[1x9088x9552] transitions (16-bit integer)
Corner Coordinates:
Upper Left  (    0.0,    0.0)
Lower Left  (    0.0,  512.0)
Upper Right (  512.0,    0.0)
Lower Right (  512.0,  512.0)
Center      (  256.0,  256.0)

In [19]:
subdatasets = dataset.GetSubDatasets()
print(subdatasets)
[('NETCDF:"/home/eoafrica/surface_flood/data/GSW.nc":change', '[1x9088x9552] change (16-bit integer)'), ('NETCDF:"/home/eoafrica/surface_flood/data/GSW.nc":extent', '[1x9088x9552] extent (16-bit integer)'), ('NETCDF:"/home/eoafrica/surface_flood/data/GSW.nc":occurrence', '[1x9088x9552] occurrence (16-bit integer)'), ('NETCDF:"/home/eoafrica/surface_flood/data/GSW.nc":recurrence', '[1x9088x9552] recurrence (16-bit integer)'), ('NETCDF:"/home/eoafrica/surface_flood/data/GSW.nc":seasonality', '[1x9088x9552] seasonality (16-bit integer)'), ('NETCDF:"/home/eoafrica/surface_flood/data/GSW.nc":transitions', '[1x9088x9552] transitions (16-bit integer)')]

Read a data subset, here 'occurrence'

In [20]:
#read second subdataset, retrieve the name as given in the first portion of the tuple, here 'NETCDF:"data/GSW.nc":occurrence'
subdataset = gdal.Open(subdatasets[2][0])

nc_info = gdal.Info(subdataset)
print (nc_info)
Driver: netCDF/Network Common Data Format
Files: /home/eoafrica/surface_flood/data/GSW.nc
Size is 9552, 9088
Coordinate System is:
GEOGCRS["WGS 84",
    DATUM["World Geodetic System 1984",
        ELLIPSOID["WGS 84",6378137,298.257223563,
            LENGTHUNIT["metre",1]],
        ID["EPSG",6326]],
    PRIMEM["Greenwich",0,
        ANGLEUNIT["degree",0.0174532925199433],
        ID["EPSG",8901]],
    CS[ellipsoidal,2],
        AXIS["longitude",east,
            ORDER[1],
            ANGLEUNIT["degree",0.0174532925199433]],
        AXIS["latitude",north,
            ORDER[2],
            ANGLEUNIT["degree",0.0174532925199433]]]
Data axis to CRS axis mapping: 1,2
Origin = (32.369999999999997,-18.669000000000000)
Pixel Size = (0.000250000000000,-0.000250000000000)
Metadata:
  crs#crs_wkt=GEOGCS["WGS 84", DATUM["World Geodetic System 1984", SPHEROID["WGS 84", 6378137.0, 298.257223563, AUTHORITY["EPSG","7030"]], AUTHORITY["EPSG","6326"]], PRIMEM["Greenwich", 0.0, AUTHORITY["EPSG","8901"]], UNIT["degree", 0.017453292519943295], AXIS["Geodetic longitude", EAST], AXIS["Geodetic latitude", NORTH], AUTHORITY["EPSG","4326"]]
  crs#spatial_ref=GEOGCS["WGS 84", DATUM["World Geodetic System 1984", SPHEROID["WGS 84", 6378137.0, 298.257223563, AUTHORITY["EPSG","7030"]], AUTHORITY["EPSG","6326"]], PRIMEM["Greenwich", 0.0, AUTHORITY["EPSG","8901"]], UNIT["degree", 0.017453292519943295], AXIS["Geodetic longitude", EAST], AXIS["Geodetic latitude", NORTH], AUTHORITY["EPSG","4326"]]
  NC_GLOBAL#Conventions=CF-1.9
  NC_GLOBAL#institution=openEO platform
  NETCDF_DIM_EXTRA={t}
  NETCDF_DIM_t_DEF={1,4}
  NETCDF_DIM_t_VALUES=10592
  occurrence#grid_mapping=crs
  occurrence#long_name=occurrence
  occurrence#units=
  occurrence#_FillValue=0
  occurrence#_Unsigned=true
  t#axis=T
  t#long_name=t
  t#standard_name=t
  t#units=days since 1990-01-01
  x#long_name=longitude
  x#standard_name=longitude
  x#units=degrees_east
  y#long_name=latitude
  y#standard_name=latitude
  y#units=degrees_north
Corner Coordinates:
Upper Left  (  32.3700000, -18.6690000) ( 32d22'12.00"E, 18d40' 8.40"S)
Lower Left  (  32.3700000, -20.9410000) ( 32d22'12.00"E, 20d56'27.60"S)
Upper Right (  34.7580000, -18.6690000) ( 34d45'28.80"E, 18d40' 8.40"S)
Lower Right (  34.7580000, -20.9410000) ( 34d45'28.80"E, 20d56'27.60"S)
Center      (  33.5640000, -19.8050000) ( 33d33'50.40"E, 19d48'18.00"S)
Band 1 Block=256x256 Type=Int16, ColorInterp=Undefined
  NoData Value=0
  Metadata:
    grid_mapping=crs
    long_name=occurrence
    NETCDF_DIM_t=10592
    NETCDF_VARNAME=occurrence
    units=
    _FillValue=0
    _Unsigned=true

In [21]:
data_array = subdataset.ReadAsArray()
In [22]:
print(data_array.min(), data_array.max())
0 100
In [23]:
plt.figure(figsize=(8,10))
plt.imshow(data_array, cmap = 'jet', vmin=0, vmax=10)
plt.title('Water extent')
plt.colorbar(shrink=0.4)
plt.show();
In [24]:
subdatasets = dataset.GetSubDatasets()
print(subdatasets)
[('NETCDF:"/home/eoafrica/surface_flood/data/GSW.nc":change', '[1x9088x9552] change (16-bit integer)'), ('NETCDF:"/home/eoafrica/surface_flood/data/GSW.nc":extent', '[1x9088x9552] extent (16-bit integer)'), ('NETCDF:"/home/eoafrica/surface_flood/data/GSW.nc":occurrence', '[1x9088x9552] occurrence (16-bit integer)'), ('NETCDF:"/home/eoafrica/surface_flood/data/GSW.nc":recurrence', '[1x9088x9552] recurrence (16-bit integer)'), ('NETCDF:"/home/eoafrica/surface_flood/data/GSW.nc":seasonality', '[1x9088x9552] seasonality (16-bit integer)'), ('NETCDF:"/home/eoafrica/surface_flood/data/GSW.nc":transitions', '[1x9088x9552] transitions (16-bit integer)')]
In [25]:
#read second subdataset, retrieve the name as given in the first portion of the tuple, here 'NETCDF:"data/GSW.nc":occurrence'
subdataset = gdal.Open(subdatasets[1][0])

nc_info = gdal.Info(subdataset)
print (nc_info)
Driver: netCDF/Network Common Data Format
Files: /home/eoafrica/surface_flood/data/GSW.nc
Size is 9552, 9088
Coordinate System is:
GEOGCRS["WGS 84",
    DATUM["World Geodetic System 1984",
        ELLIPSOID["WGS 84",6378137,298.257223563,
            LENGTHUNIT["metre",1]],
        ID["EPSG",6326]],
    PRIMEM["Greenwich",0,
        ANGLEUNIT["degree",0.0174532925199433],
        ID["EPSG",8901]],
    CS[ellipsoidal,2],
        AXIS["longitude",east,
            ORDER[1],
            ANGLEUNIT["degree",0.0174532925199433]],
        AXIS["latitude",north,
            ORDER[2],
            ANGLEUNIT["degree",0.0174532925199433]]]
Data axis to CRS axis mapping: 1,2
Origin = (32.369999999999997,-18.669000000000000)
Pixel Size = (0.000250000000000,-0.000250000000000)
Metadata:
  crs#crs_wkt=GEOGCS["WGS 84", DATUM["World Geodetic System 1984", SPHEROID["WGS 84", 6378137.0, 298.257223563, AUTHORITY["EPSG","7030"]], AUTHORITY["EPSG","6326"]], PRIMEM["Greenwich", 0.0, AUTHORITY["EPSG","8901"]], UNIT["degree", 0.017453292519943295], AXIS["Geodetic longitude", EAST], AXIS["Geodetic latitude", NORTH], AUTHORITY["EPSG","4326"]]
  crs#spatial_ref=GEOGCS["WGS 84", DATUM["World Geodetic System 1984", SPHEROID["WGS 84", 6378137.0, 298.257223563, AUTHORITY["EPSG","7030"]], AUTHORITY["EPSG","6326"]], PRIMEM["Greenwich", 0.0, AUTHORITY["EPSG","8901"]], UNIT["degree", 0.017453292519943295], AXIS["Geodetic longitude", EAST], AXIS["Geodetic latitude", NORTH], AUTHORITY["EPSG","4326"]]
  extent#grid_mapping=crs
  extent#long_name=extent
  extent#units=
  extent#_FillValue=0
  extent#_Unsigned=true
  NC_GLOBAL#Conventions=CF-1.9
  NC_GLOBAL#institution=openEO platform
  NETCDF_DIM_EXTRA={t}
  NETCDF_DIM_t_DEF={1,4}
  NETCDF_DIM_t_VALUES=10592
  t#axis=T
  t#long_name=t
  t#standard_name=t
  t#units=days since 1990-01-01
  x#long_name=longitude
  x#standard_name=longitude
  x#units=degrees_east
  y#long_name=latitude
  y#standard_name=latitude
  y#units=degrees_north
Corner Coordinates:
Upper Left  (  32.3700000, -18.6690000) ( 32d22'12.00"E, 18d40' 8.40"S)
Lower Left  (  32.3700000, -20.9410000) ( 32d22'12.00"E, 20d56'27.60"S)
Upper Right (  34.7580000, -18.6690000) ( 34d45'28.80"E, 18d40' 8.40"S)
Lower Right (  34.7580000, -20.9410000) ( 34d45'28.80"E, 20d56'27.60"S)
Center      (  33.5640000, -19.8050000) ( 33d33'50.40"E, 19d48'18.00"S)
Band 1 Block=256x256 Type=Int16, ColorInterp=Undefined
  NoData Value=0
  Metadata:
    grid_mapping=crs
    long_name=extent
    NETCDF_DIM_t=10592
    NETCDF_VARNAME=extent
    units=
    _FillValue=0
    _Unsigned=true

In [26]:
data_array = subdataset.ReadAsArray()
In [27]:
print(data_array.min(), data_array.max())
0 1
In [28]:
plt.figure(figsize=(8,10))
plt.imshow(data_array, interpolation='nearest', cmap = 'winter_r', vmin=0, vmax=1)
plt.title('Water Extent')
plt.show()

Retrieve water bodies¶

In [29]:
con.describe_collection('WATER_BODIES')#web listing
Out[29]:
In [30]:
collection_id = 'WATER_BODIES'
In [31]:
new_temporal_extent = ["2023-07-01", "2023-08-01"]
In [32]:
cube = con.load_collection(
    collection_id, 
    temporal_extent=new_temporal_extent,
    spatial_extent=spatial_extent
)
In [33]:
%%time 
cube.download(work_dir + '/WB_month.nc')
CPU times: user 9.37 ms, sys: 2.24 ms, total: 11.6 ms
Wall time: 12 s
In [34]:
input_file = (work_dir + '/WB_month.nc')
print(input_file)
/home/eoafrica/surface_flood/data/WB_month.nc
In [35]:
dataset = gdal.Open(input_file)
nc_info = gdal.Info(dataset)
print (nc_info)
Driver: netCDF/Network Common Data Format
Files: /home/eoafrica/surface_flood/data/WB_month.nc
Size is 512, 512
Metadata:
  NC_GLOBAL#Conventions=CF-1.9
  NC_GLOBAL#institution=openEO platform
Subdatasets:
  SUBDATASET_1_NAME=NETCDF:"/home/eoafrica/surface_flood/data/WB_month.nc":WB
  SUBDATASET_1_DESC=[1x2454x2580] WB (16-bit integer)
  SUBDATASET_2_NAME=NETCDF:"/home/eoafrica/surface_flood/data/WB_month.nc":QUAL
  SUBDATASET_2_DESC=[1x2454x2580] QUAL (16-bit integer)
  SUBDATASET_3_NAME=NETCDF:"/home/eoafrica/surface_flood/data/WB_month.nc":dataMask
  SUBDATASET_3_DESC=[1x2454x2580] dataMask (16-bit integer)
Corner Coordinates:
Upper Left  (    0.0,    0.0)
Lower Left  (    0.0,  512.0)
Upper Right (  512.0,    0.0)
Lower Right (  512.0,  512.0)
Center      (  256.0,  256.0)

In [36]:
subdatasets = dataset.GetSubDatasets()
print(subdatasets)
[('NETCDF:"/home/eoafrica/surface_flood/data/WB_month.nc":WB', '[1x2454x2580] WB (16-bit integer)'), ('NETCDF:"/home/eoafrica/surface_flood/data/WB_month.nc":QUAL', '[1x2454x2580] QUAL (16-bit integer)'), ('NETCDF:"/home/eoafrica/surface_flood/data/WB_month.nc":dataMask', '[1x2454x2580] dataMask (16-bit integer)')]
In [37]:
#read first subdataset, retrieve the name as given in the first portion of the tuple, here 'NETCDF:"data/WB_month.nc":WB'
subdataset = gdal.Open(subdatasets[0][0])

nc_info = gdal.Info(subdataset)
print (nc_info)
Driver: netCDF/Network Common Data Format
Files: /home/eoafrica/surface_flood/data/WB_month.nc
Size is 2580, 2454
Coordinate System is:
GEOGCRS["WGS 84",
    DATUM["World Geodetic System 1984",
        ELLIPSOID["WGS 84",6378137,298.257223563,
            LENGTHUNIT["metre",1]],
        ID["EPSG",6326]],
    PRIMEM["Greenwich",0,
        ANGLEUNIT["degree",0.0174532925199433],
        ID["EPSG",8901]],
    CS[ellipsoidal,2],
        AXIS["longitude",east,
            ORDER[1],
            ANGLEUNIT["degree",0.0174532925199433]],
        AXIS["latitude",north,
            ORDER[2],
            ANGLEUNIT["degree",0.0174532925199433]]]
Data axis to CRS axis mapping: 1,2
Origin = (32.369999999999997,-18.669000000000000)
Pixel Size = (0.000925925925926,-0.000925925925926)
Metadata:
  crs#crs_wkt=GEOGCS["WGS 84", DATUM["World Geodetic System 1984", SPHEROID["WGS 84", 6378137.0, 298.257223563, AUTHORITY["EPSG","7030"]], AUTHORITY["EPSG","6326"]], PRIMEM["Greenwich", 0.0, AUTHORITY["EPSG","8901"]], UNIT["degree", 0.017453292519943295], AXIS["Geodetic longitude", EAST], AXIS["Geodetic latitude", NORTH], AUTHORITY["EPSG","4326"]]
  crs#spatial_ref=GEOGCS["WGS 84", DATUM["World Geodetic System 1984", SPHEROID["WGS 84", 6378137.0, 298.257223563, AUTHORITY["EPSG","7030"]], AUTHORITY["EPSG","6326"]], PRIMEM["Greenwich", 0.0, AUTHORITY["EPSG","8901"]], UNIT["degree", 0.017453292519943295], AXIS["Geodetic longitude", EAST], AXIS["Geodetic latitude", NORTH], AUTHORITY["EPSG","4326"]]
  NC_GLOBAL#Conventions=CF-1.9
  NC_GLOBAL#institution=openEO platform
  NETCDF_DIM_EXTRA={t}
  NETCDF_DIM_t_DEF={1,4}
  NETCDF_DIM_t_VALUES=12234
  t#axis=T
  t#long_name=t
  t#standard_name=t
  t#units=days since 1990-01-01
  WB#grid_mapping=crs
  WB#long_name=WB
  WB#units=
  WB#_FillValue=0
  WB#_Unsigned=true
  x#long_name=longitude
  x#standard_name=longitude
  x#units=degrees_east
  y#long_name=latitude
  y#standard_name=latitude
  y#units=degrees_north
Corner Coordinates:
Upper Left  (  32.3700000, -18.6690000) ( 32d22'12.00"E, 18d40' 8.40"S)
Lower Left  (  32.3700000, -20.9412222) ( 32d22'12.00"E, 20d56'28.40"S)
Upper Right (  34.7588889, -18.6690000) ( 34d45'32.00"E, 18d40' 8.40"S)
Lower Right (  34.7588889, -20.9412222) ( 34d45'32.00"E, 20d56'28.40"S)
Center      (  33.5644444, -19.8051111) ( 33d33'52.00"E, 19d48'18.40"S)
Band 1 Block=256x256 Type=Int16, ColorInterp=Undefined
  NoData Value=0
  Metadata:
    grid_mapping=crs
    long_name=WB
    NETCDF_DIM_t=12234
    NETCDF_VARNAME=WB
    units=
    _FillValue=0
    _Unsigned=true

In [38]:
data_array = subdataset.ReadAsArray().astype(np.int16)
print(data_array.dtype)
int16
In [39]:
print(data_array.min(), data_array.max())
0 251
In [40]:
plt.figure(figsize=(8,10))
plt.imshow(data_array, cmap = 'winter', interpolation='nearest')#, vmin=data_array.min(), vmax=data_array.max())
plt.title('Water Bodies')
#plt.colorbar(shrink=0.4)
plt.grid(False)
plt.show()
In [ ]: