In [18]:
# 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
In [19]:
import openeo
openeo.client_version()
Out[19]:
'0.25.0'
In [20]:
!rm -rf /home/eoafrica/.local/share/openeo-python-client/refresh-tokens.json
In [21]:
con = openeo.connect(url="openeo.vito.be")
In [23]:
con.authenticate_oidc()
Authenticated using refresh token.
Out[23]:
<Connection to 'https://openeo.vito.be/openeo/1.2/' with OidcBearerAuth>
In [24]:
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']
In [25]:
con.describe_collection('GLOBAL_LAND_COVER')#web listing
Out[25]:
In [26]:
collection_id = 'GLOBAL_LAND_COVER'
In [27]:
spatial_extent = {"west": 32.37, "south": -20.941, "east": 34.758, "north": -18.669}
In [28]:
cube = con.load_collection(
    "GLOBAL_LAND_COVER",
    spatial_extent=spatial_extent,
    temporal_extent=["2015-01-01", "2019-12-31"],
    bands=["Discrete_Classification_map"],
)
In [38]:
%%time 
cube.download('data/Buzi_LCC.nc')
In [36]:
input_file = (r'data/Buzi_LCC.nc')
print(input_file)
data/Buzi_LCC.nc
In [37]:
dataset = gdal.Open(input_file)
nc_info = gdal.Info(dataset)
print (nc_info)
Driver: netCDF/Network Common Data Format
Files: data/Buzi_LCC.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"]]
  Discrete_Classification_map#grid_mapping=crs
  Discrete_Classification_map#long_name=Discrete_Classification_map
  Discrete_Classification_map#units=
  Discrete_Classification_map#_FillValue=0
  Discrete_Classification_map#_Unsigned=true
  NC_GLOBAL#Conventions=CF-1.9
  NC_GLOBAL#institution=openEO platform
  NETCDF_DIM_EXTRA={t}
  NETCDF_DIM_t_DEF={5,4}
  NETCDF_DIM_t_VALUES={9131,9496,9862,10227,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.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=Discrete_Classification_map
    NETCDF_DIM_t=9131
    NETCDF_VARNAME=Discrete_Classification_map
    units=
    _FillValue=0
    _Unsigned=true
Band 2 Block=256x256 Type=Int16, ColorInterp=Undefined
  NoData Value=0
  Metadata:
    grid_mapping=crs
    long_name=Discrete_Classification_map
    NETCDF_DIM_t=9496
    NETCDF_VARNAME=Discrete_Classification_map
    units=
    _FillValue=0
    _Unsigned=true
Band 3 Block=256x256 Type=Int16, ColorInterp=Undefined
  NoData Value=0
  Metadata:
    grid_mapping=crs
    long_name=Discrete_Classification_map
    NETCDF_DIM_t=9862
    NETCDF_VARNAME=Discrete_Classification_map
    units=
    _FillValue=0
    _Unsigned=true
Band 4 Block=256x256 Type=Int16, ColorInterp=Undefined
  NoData Value=0
  Metadata:
    grid_mapping=crs
    long_name=Discrete_Classification_map
    NETCDF_DIM_t=10227
    NETCDF_VARNAME=Discrete_Classification_map
    units=
    _FillValue=0
    _Unsigned=true
Band 5 Block=256x256 Type=Int16, ColorInterp=Undefined
  NoData Value=0
  Metadata:
    grid_mapping=crs
    long_name=Discrete_Classification_map
    NETCDF_DIM_t=10592
    NETCDF_VARNAME=Discrete_Classification_map
    units=
    _FillValue=0
    _Unsigned=true

In [39]:
data_array = dataset.ReadAsArray().astype(np.int16)
print(data_array.dtype)
int16
In [40]:
print(data_array.min(), data_array.max())
0 200
In [41]:
plt.figure(figsize=(8,10))
plt.imshow(data_array[0], cmap = 'jet', interpolation='nearest')
plt.title('Land Cover CLasses')
plt.colorbar(shrink=0.4)
plt.grid(False)
plt.show()
In [42]:
#histogram
plt.hist(data_array[0], bins=10)
plt.show()
In [ ]: