How to work with the Copernicus DEM data on Earth Data Hub¶

This notebook will provide you guidance on how to access and use the https://data.earthdatahub.destine.eu/copernicus-dem/GLO-30-v0.zarr dataset on Earth Data Hub (EDH). The Copernicus DEM is a Digital Surface Model (DSM) representing the surface of the Earth including buildings, infrastructure and vegetation.

Goal of this tutorial¶

The goal of this tutorial is to visualize the copernicus DEM for the European area.

What you will learn:¶

  • how to access the data
  • select and reduce the data
  • plot the results

To access datasets on Earth Data Hub you need to instruct your tools (xarray, Zarr, etc.) to use EDH personal access token when downloading the data.

To obtain a personal access token you first need to register to the Destination Earth platform. Then you can go to Earth Data Hub account settings where you can find your default personal access token or create others. After retrieving your personal access token, please cut and paste it below: ⤵

In [1]:
PAT = "your_personal_access_token"

# e.g. PAT="edh_pat_44bbb7e9192a4c6bb47ddf07d07564eee5d17de8dfc48f7118f88e3bc4a4157f8fe2403f5aa0a2d53441b6922ea9a33a"

Working with EDH data¶

Datasets on EDH are often very large and remotely hosted. Typical use imply a selection of the data followed by one or more reduction steps to be performed in a local or distributed Dask environment.

The structure of a workflow that uses EDH data tipically looks like this:

  • data access
  • data selection
  • (optional) data reduction
  • data download
  • further operations and visualization

Xarray and Dask work together following a lazy principle. This means that when you access and manipulate a Zarr store the data is in not immediately downloaded and loaded in memory. Instead, Dask constructs a task graph that represents the operations to be performed. A smart user will first reduce the amount of data that needs to be downloaded and explicitly call compute() on it. Once the compute() operation is complete the data is loaded into memory and available for subsequent fast processing.

1. Data access¶

To access the data, only the dataset metadata must be downloaded. Xarray does this automatically when you access a Zarr dataset:

In [2]:
import xarray as xr

ds = xr.open_dataset(
    f"https://edh:{PAT}@data.earthdatahub.destine.eu/copernicus-dem/GLO-30-v0.zarr",
    chunks={},
    engine="zarr",
    decode_coords="all",
    mask_and_scale=False,
)
ds
Out[2]:
<xarray.Dataset> Size: 3TB
Dimensions:      (lat: 648000, lon: 1296001)
Coordinates:
  * lat          (lat) float64 5MB -90.0 -90.0 -90.0 -90.0 ... 90.0 90.0 90.0
  * lon          (lon) float64 10MB -180.0 -180.0 -180.0 ... 180.0 180.0 180.0
    spatial_ref  int64 8B ...
Data variables:
    dsm          (lat, lon) float32 3TB dask.array<chunksize=(3600, 3600), meta=np.ndarray>
xarray.Dataset
    • lat: 648000
    • lon: 1296001
    • lat
      (lat)
      float64
      -90.0 -90.0 -90.0 ... 90.0 90.0
      long_name :
      latitude
      standard_name :
      latitude
      units :
      degrees_north
      _FillValue :
      nan
      array([-89.999722, -89.999444, -89.999167, ...,  89.999444,  89.999722,
              90.      ])
    • lon
      (lon)
      float64
      -180.0 -180.0 ... 180.0 180.0
      long_name :
      longitude
      standard_name :
      longitude
      units :
      degrees_east
      _FillValue :
      nan
      array([-180.      , -179.999722, -179.999444, ...,  179.999444,  179.999722,
              180.      ])
    • spatial_ref
      ()
      int64
      ...
      crs_wkt :
      COMPD_CS["WGS 84 + EGM2008 height",GEOGCS["WGS 84",DATUM["WGS_1984",SPHEROID["WGS 84",6378137,298.257223563,AUTHORITY["EPSG","7030"]],AUTHORITY["EPSG","6326"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]],AXIS["Latitude",NORTH],AXIS["Longitude",EAST],AUTHORITY["EPSG","4326"]],VERT_CS["EGM2008 height",VERT_DATUM["EGM2008 geoid",2005,AUTHORITY["EPSG","1027"]],UNIT["metre",1,AUTHORITY["EPSG","9001"]],AXIS["Gravity-related height",UP],AUTHORITY["EPSG","3855"]],AUTHORITY["EPSG","9518"]]
      geographic_crs_name :
      WGS 84
      geopotential_datum_name :
      EGM2008 geoid
      grid_mapping_name :
      latitude_longitude
      horizontal_datum_name :
      World Geodetic System 1984
      inverse_flattening :
      298.257223563
      longitude_of_prime_meridian :
      0.0
      prime_meridian_name :
      Greenwich
      reference_ellipsoid_name :
      WGS 84
      semi_major_axis :
      6378137.0
      semi_minor_axis :
      6356752.314245179
      spatial_ref :
      COMPD_CS["WGS 84 + EGM2008 height",GEOGCS["WGS 84",DATUM["WGS_1984",SPHEROID["WGS 84",6378137,298.257223563,AUTHORITY["EPSG","7030"]],AUTHORITY["EPSG","6326"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]],AXIS["Latitude",NORTH],AXIS["Longitude",EAST],AUTHORITY["EPSG","4326"]],VERT_CS["EGM2008 height",VERT_DATUM["EGM2008 geoid",2005,AUTHORITY["EPSG","1027"]],UNIT["metre",1,AUTHORITY["EPSG","9001"]],AXIS["Gravity-related height",UP],AUTHORITY["EPSG","3855"]],AUTHORITY["EPSG","9518"]]
      [1 values with dtype=int64]
    • dsm
      (lat, lon)
      float32
      dask.array<chunksize=(3600, 3600), meta=np.ndarray>
      long_name :
      height above geoid
      standard_name :
      altitude
      units :
      m
      _FillValue :
      0.0
      Array Chunk
      Bytes 3.06 TiB 49.44 MiB
      Shape (648000, 1296001) (3600, 3600)
      Dask graph 64980 chunks in 2 graph layers
      Data type float32 numpy.ndarray
      1296001 648000
    • lat
      PandasIndex
      PandasIndex(Index([-89.99972222222222, -89.99944444444445, -89.99916666666667,
             -89.99888888888889, -89.99861111111112, -89.99833333333333,
             -89.99805555555555, -89.99777777777778,           -89.9975,
             -89.99722222222222,
             ...
                        89.9975,  89.99777777777778,  89.99805555555555,
              89.99833333333333,  89.99861111111112,  89.99888888888889,
              89.99916666666667,  89.99944444444445,  89.99972222222222,
                           90.0],
            dtype='float64', name='lat', length=648000))
    • lon
      PandasIndex
      PandasIndex(Index([             -180.0, -179.99972222222223, -179.99944444444444,
             -179.99916666666667,  -179.9988888888889,  -179.9986111111111,
             -179.99833333333333, -179.99805555555557, -179.99777777777777,
                       -179.9975,
             ...
                        179.9975,  179.99777777777777,  179.99805555555557,
              179.99833333333333,   179.9986111111111,   179.9988888888889,
              179.99916666666667,  179.99944444444444,  179.99972222222223,
                           180.0],
            dtype='float64', name='lon', length=1296001))

âš  At this point, no data has been downloaded yet, nor loaded in memory.

2. Data selection¶

From the Copernicus DEM dataset we select the dsm (digital surface model) variable, and perform a geographical selection corresponding to the European area:

In [3]:
dsm = ds.dsm
europe = dsm.sel(lat=slice(36, 71), lon=slice(-10, 41))
europe
Out[3]:
<xarray.DataArray 'dsm' (lat: 126001, lon: 183601)> Size: 93GB
dask.array<getitem, shape=(126001, 183601), dtype=float32, chunksize=(3600, 3600), chunktype=numpy.ndarray>
Coordinates:
  * lat          (lat) float64 1MB 36.0 36.0 36.0 36.0 ... 71.0 71.0 71.0 71.0
  * lon          (lon) float64 1MB -10.0 -10.0 -9.999 -9.999 ... 41.0 41.0 41.0
    spatial_ref  int64 8B ...
Attributes:
    long_name:      height above geoid
    standard_name:  altitude
    units:          m
    _FillValue:     0.0
xarray.DataArray
'dsm'
  • lat: 126001
  • lon: 183601
  • dask.array<chunksize=(1, 3600), meta=np.ndarray>
    Array Chunk
    Bytes 86.18 GiB 49.44 MiB
    Shape (126001, 183601) (3600, 3600)
    Dask graph 1872 chunks in 3 graph layers
    Data type float32 numpy.ndarray
    183601 126001
    • lat
      (lat)
      float64
      36.0 36.0 36.0 ... 71.0 71.0 71.0
      long_name :
      latitude
      standard_name :
      latitude
      units :
      degrees_north
      _FillValue :
      nan
      array([36.      , 36.000278, 36.000556, ..., 70.999444, 70.999722, 71.      ])
    • lon
      (lon)
      float64
      -10.0 -10.0 -9.999 ... 41.0 41.0
      long_name :
      longitude
      standard_name :
      longitude
      units :
      degrees_east
      _FillValue :
      nan
      array([-10.      ,  -9.999722,  -9.999444, ...,  40.999444,  40.999722,
              41.      ])
    • spatial_ref
      ()
      int64
      ...
      crs_wkt :
      COMPD_CS["WGS 84 + EGM2008 height",GEOGCS["WGS 84",DATUM["WGS_1984",SPHEROID["WGS 84",6378137,298.257223563,AUTHORITY["EPSG","7030"]],AUTHORITY["EPSG","6326"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]],AXIS["Latitude",NORTH],AXIS["Longitude",EAST],AUTHORITY["EPSG","4326"]],VERT_CS["EGM2008 height",VERT_DATUM["EGM2008 geoid",2005,AUTHORITY["EPSG","1027"]],UNIT["metre",1,AUTHORITY["EPSG","9001"]],AXIS["Gravity-related height",UP],AUTHORITY["EPSG","3855"]],AUTHORITY["EPSG","9518"]]
      geographic_crs_name :
      WGS 84
      geopotential_datum_name :
      EGM2008 geoid
      grid_mapping_name :
      latitude_longitude
      horizontal_datum_name :
      World Geodetic System 1984
      inverse_flattening :
      298.257223563
      longitude_of_prime_meridian :
      0.0
      prime_meridian_name :
      Greenwich
      reference_ellipsoid_name :
      WGS 84
      semi_major_axis :
      6378137.0
      semi_minor_axis :
      6356752.314245179
      spatial_ref :
      COMPD_CS["WGS 84 + EGM2008 height",GEOGCS["WGS 84",DATUM["WGS_1984",SPHEROID["WGS 84",6378137,298.257223563,AUTHORITY["EPSG","7030"]],AUTHORITY["EPSG","6326"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]],AXIS["Latitude",NORTH],AXIS["Longitude",EAST],AUTHORITY["EPSG","4326"]],VERT_CS["EGM2008 height",VERT_DATUM["EGM2008 geoid",2005,AUTHORITY["EPSG","1027"]],UNIT["metre",1,AUTHORITY["EPSG","9001"]],AXIS["Gravity-related height",UP],AUTHORITY["EPSG","3855"]],AUTHORITY["EPSG","9518"]]
      [1 values with dtype=int64]
    • lat
      PandasIndex
      PandasIndex(Index([              36.0, 36.000277777777775,  36.00055555555556,
              36.00083333333333,  36.00111111111111,  36.00138888888889,
             36.001666666666665,  36.00194444444445,  36.00222222222222,
                        36.0025,
             ...
                        70.9975,  70.99777777777778,  70.99805555555555,
              70.99833333333333,  70.99861111111112,  70.99888888888889,
              70.99916666666667,  70.99944444444445,  70.99972222222222,
                           71.0],
            dtype='float64', name='lat', length=126001))
    • lon
      PandasIndex
      PandasIndex(Index([             -10.0, -9.999722222222223, -9.999444444444444,
             -9.999166666666667, -9.998888888888889, -9.998611111111112,
             -9.998333333333333, -9.998055555555556, -9.997777777777777,
                        -9.9975,
             ...
                        40.9975,  40.99777777777778,  40.99805555555555,
             40.998333333333335,  40.99861111111111,  40.99888888888889,
              40.99916666666667,  40.99944444444444, 40.999722222222225,
                           41.0],
            dtype='float64', name='lon', length=183601))
  • long_name :
    height above geoid
    standard_name :
    altitude
    units :
    m
    _FillValue :
    0.0

This greatly reduces the amount of data that will be downloaded from Earth Data Hub. However, due to the chunked structure of the DataArray, xarray must still download every chunk that includes a portion of the selected data.

To estimate the size of the download, we can use the costing.py module. This estimate must be done before we apply any reduction operation.

In [5]:
import costing

costing.estimate_download_size(dsm, reduced_europe)
estimated_needed_chunks: 1836
estimated_memory_size: 95.178 GB
estimated_download_size: 9.518 GB

3. Data reduction¶

Copernicus DEM on Earth Data Hub is a high resoluton regular latitide-longitude grid dataset. For the scope or this tutorial, we do not need its full precision. We can downscale the dataset resolution by selecting only one grid point out of 1000:

In [4]:
reduced_europe = europe[::1000, ::1000]
reduced_europe
Out[4]:
<xarray.DataArray 'dsm' (lat: 127, lon: 184)> Size: 93kB
dask.array<getitem, shape=(127, 184), dtype=float32, chunksize=(4, 4), chunktype=numpy.ndarray>
Coordinates:
  * lat          (lat) float64 1kB 36.0 36.28 36.56 36.83 ... 70.44 70.72 71.0
  * lon          (lon) float64 1kB -10.0 -9.722 -9.444 ... 40.28 40.56 40.83
    spatial_ref  int64 8B ...
Attributes:
    long_name:      height above geoid
    standard_name:  altitude
    units:          m
    _FillValue:     0.0
xarray.DataArray
'dsm'
  • lat: 127
  • lon: 184
  • dask.array<chunksize=(1, 4), meta=np.ndarray>
    Array Chunk
    Bytes 91.28 kiB 64 B
    Shape (127, 184) (4, 4)
    Dask graph 1836 chunks in 4 graph layers
    Data type float32 numpy.ndarray
    184 127
    • lat
      (lat)
      float64
      36.0 36.28 36.56 ... 70.72 71.0
      long_name :
      latitude
      standard_name :
      latitude
      units :
      degrees_north
      _FillValue :
      nan
      array([36.      , 36.277778, 36.555556, 36.833333, 37.111111, 37.388889,
             37.666667, 37.944444, 38.222222, 38.5     , 38.777778, 39.055556,
             39.333333, 39.611111, 39.888889, 40.166667, 40.444444, 40.722222,
             41.      , 41.277778, 41.555556, 41.833333, 42.111111, 42.388889,
             42.666667, 42.944444, 43.222222, 43.5     , 43.777778, 44.055556,
             44.333333, 44.611111, 44.888889, 45.166667, 45.444444, 45.722222,
             46.      , 46.277778, 46.555556, 46.833333, 47.111111, 47.388889,
             47.666667, 47.944444, 48.222222, 48.5     , 48.777778, 49.055556,
             49.333333, 49.611111, 49.888889, 50.166667, 50.444444, 50.722222,
             51.      , 51.277778, 51.555556, 51.833333, 52.111111, 52.388889,
             52.666667, 52.944444, 53.222222, 53.5     , 53.777778, 54.055556,
             54.333333, 54.611111, 54.888889, 55.166667, 55.444444, 55.722222,
             56.      , 56.277778, 56.555556, 56.833333, 57.111111, 57.388889,
             57.666667, 57.944444, 58.222222, 58.5     , 58.777778, 59.055556,
             59.333333, 59.611111, 59.888889, 60.166667, 60.444444, 60.722222,
             61.      , 61.277778, 61.555556, 61.833333, 62.111111, 62.388889,
             62.666667, 62.944444, 63.222222, 63.5     , 63.777778, 64.055556,
             64.333333, 64.611111, 64.888889, 65.166667, 65.444444, 65.722222,
             66.      , 66.277778, 66.555556, 66.833333, 67.111111, 67.388889,
             67.666667, 67.944444, 68.222222, 68.5     , 68.777778, 69.055556,
             69.333333, 69.611111, 69.888889, 70.166667, 70.444444, 70.722222,
             71.      ])
    • lon
      (lon)
      float64
      -10.0 -9.722 -9.444 ... 40.56 40.83
      long_name :
      longitude
      standard_name :
      longitude
      units :
      degrees_east
      _FillValue :
      nan
      array([-10.      ,  -9.722222,  -9.444444,  -9.166667,  -8.888889,  -8.611111,
              -8.333333,  -8.055556,  -7.777778,  -7.5     ,  -7.222222,  -6.944444,
              -6.666667,  -6.388889,  -6.111111,  -5.833333,  -5.555556,  -5.277778,
              -5.      ,  -4.722222,  -4.444444,  -4.166667,  -3.888889,  -3.611111,
              -3.333333,  -3.055556,  -2.777778,  -2.5     ,  -2.222222,  -1.944444,
              -1.666667,  -1.388889,  -1.111111,  -0.833333,  -0.555556,  -0.277778,
               0.      ,   0.277778,   0.555556,   0.833333,   1.111111,   1.388889,
               1.666667,   1.944444,   2.222222,   2.5     ,   2.777778,   3.055556,
               3.333333,   3.611111,   3.888889,   4.166667,   4.444444,   4.722222,
               5.      ,   5.277778,   5.555556,   5.833333,   6.111111,   6.388889,
               6.666667,   6.944444,   7.222222,   7.5     ,   7.777778,   8.055556,
               8.333333,   8.611111,   8.888889,   9.166667,   9.444444,   9.722222,
              10.      ,  10.277778,  10.555556,  10.833333,  11.111111,  11.388889,
              11.666667,  11.944444,  12.222222,  12.5     ,  12.777778,  13.055556,
              13.333333,  13.611111,  13.888889,  14.166667,  14.444444,  14.722222,
              15.      ,  15.277778,  15.555556,  15.833333,  16.111111,  16.388889,
              16.666667,  16.944444,  17.222222,  17.5     ,  17.777778,  18.055556,
              18.333333,  18.611111,  18.888889,  19.166667,  19.444444,  19.722222,
              20.      ,  20.277778,  20.555556,  20.833333,  21.111111,  21.388889,
              21.666667,  21.944444,  22.222222,  22.5     ,  22.777778,  23.055556,
              23.333333,  23.611111,  23.888889,  24.166667,  24.444444,  24.722222,
              25.      ,  25.277778,  25.555556,  25.833333,  26.111111,  26.388889,
              26.666667,  26.944444,  27.222222,  27.5     ,  27.777778,  28.055556,
              28.333333,  28.611111,  28.888889,  29.166667,  29.444444,  29.722222,
              30.      ,  30.277778,  30.555556,  30.833333,  31.111111,  31.388889,
              31.666667,  31.944444,  32.222222,  32.5     ,  32.777778,  33.055556,
              33.333333,  33.611111,  33.888889,  34.166667,  34.444444,  34.722222,
              35.      ,  35.277778,  35.555556,  35.833333,  36.111111,  36.388889,
              36.666667,  36.944444,  37.222222,  37.5     ,  37.777778,  38.055556,
              38.333333,  38.611111,  38.888889,  39.166667,  39.444444,  39.722222,
              40.      ,  40.277778,  40.555556,  40.833333])
    • spatial_ref
      ()
      int64
      ...
      crs_wkt :
      COMPD_CS["WGS 84 + EGM2008 height",GEOGCS["WGS 84",DATUM["WGS_1984",SPHEROID["WGS 84",6378137,298.257223563,AUTHORITY["EPSG","7030"]],AUTHORITY["EPSG","6326"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]],AXIS["Latitude",NORTH],AXIS["Longitude",EAST],AUTHORITY["EPSG","4326"]],VERT_CS["EGM2008 height",VERT_DATUM["EGM2008 geoid",2005,AUTHORITY["EPSG","1027"]],UNIT["metre",1,AUTHORITY["EPSG","9001"]],AXIS["Gravity-related height",UP],AUTHORITY["EPSG","3855"]],AUTHORITY["EPSG","9518"]]
      geographic_crs_name :
      WGS 84
      geopotential_datum_name :
      EGM2008 geoid
      grid_mapping_name :
      latitude_longitude
      horizontal_datum_name :
      World Geodetic System 1984
      inverse_flattening :
      298.257223563
      longitude_of_prime_meridian :
      0.0
      prime_meridian_name :
      Greenwich
      reference_ellipsoid_name :
      WGS 84
      semi_major_axis :
      6378137.0
      semi_minor_axis :
      6356752.314245179
      spatial_ref :
      COMPD_CS["WGS 84 + EGM2008 height",GEOGCS["WGS 84",DATUM["WGS_1984",SPHEROID["WGS 84",6378137,298.257223563,AUTHORITY["EPSG","7030"]],AUTHORITY["EPSG","6326"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]],AXIS["Latitude",NORTH],AXIS["Longitude",EAST],AUTHORITY["EPSG","4326"]],VERT_CS["EGM2008 height",VERT_DATUM["EGM2008 geoid",2005,AUTHORITY["EPSG","1027"]],UNIT["metre",1,AUTHORITY["EPSG","9001"]],AXIS["Gravity-related height",UP],AUTHORITY["EPSG","3855"]],AUTHORITY["EPSG","9518"]]
      [1 values with dtype=int64]
    • lat
      PandasIndex
      PandasIndex(Index([              36.0,  36.27777777777778,  36.55555555555556,
             36.833333333333336, 37.111111111111114, 37.388888888888886,
             37.666666666666664,  37.94444444444444,  38.22222222222222,
                           38.5,
             ...
                           68.5,  68.77777777777777,  69.05555555555556,
              69.33333333333333,  69.61111111111111,  69.88888888888889,
              70.16666666666667,  70.44444444444444,  70.72222222222223,
                           71.0],
            dtype='float64', name='lat', length=127))
    • lon
      PandasIndex
      PandasIndex(Index([             -10.0, -9.722222222222221, -9.444444444444445,
             -9.166666666666666,  -8.88888888888889,  -8.61111111111111,
             -8.333333333333334, -8.055555555555555, -7.777777777777778,
                           -7.5,
             ...
             38.333333333333336, 38.611111111111114, 38.888888888888886,
             39.166666666666664,  39.44444444444444,  39.72222222222222,
                           40.0,  40.27777777777778,  40.55555555555556,
             40.833333333333336],
            dtype='float64', name='lon', length=184))
  • long_name :
    height above geoid
    standard_name :
    altitude
    units :
    m
    _FillValue :
    0.0

4. Data download¶

This is the phase where we explicitly trigger the download of the data. Remember to assign the return of the compute() function to a new variable, so that the data is kept in memory.

In [6]:
%%time

reduced_europe_computed = reduced_europe.compute()
CPU times: user 3min 17s, sys: 1min 41s, total: 4min 58s
Wall time: 3min 46s

3. Visualization¶

Now that the data has been loaded in memory, we can visualize the height above geoid for continental Europe. We also add borders, coastlines and ocean. The ocean feature will help us to to distinguish zero quote land (or lower) from the sea.

In [7]:
import numpy as np
import matplotlib.pyplot as plt
import cartopy.crs as ccrs
import cartopy.feature as cfeature

lon = reduced_europe_computed.lon
lat = reduced_europe_computed.lat

# Create the figure and axis with Mercator projection
fig, ax = plt.subplots(figsize=(10, 10), subplot_kw={"projection": ccrs.AlbersEqualArea()})

# Add map features (coastlines, borders, etc.)
ax.add_feature(cfeature.COASTLINE, linewidth=0.5)
ax.add_feature(cfeature.BORDERS, linestyle=':', linewidth=0.5)
ax.add_feature(cfeature.OCEAN, facecolor='lightblue', zorder=2)

# Set the extent of the map
ax.set_extent([min(lon), max(lon)-13, min(lat), max(lat)], crs=ccrs.PlateCarree())

contourf = ax.contourf(lon, lat, reduced_europe_computed, levels=100, cmap='terrain', transform=ccrs.PlateCarree(), vmin=-2000, vmax=4000)

# Add colorbar
cbar = plt.colorbar(contourf, ax=ax, orientation='vertical', pad=0.07)
cbar.set_label('Altitude (m)')  # Label for the colorbar

# Add gridlines for lat/lon
ax.gridlines(draw_labels=True, zorder=3)

# Show the plot
plt.show()
No description has been provided for this image