Integrate landcover info to aggregate accuracy statistics

[loreabad6]

Landcover can be integrated within the final accuracy statistics of the generated DEMs. This would serve to have an overview of how the error is influenced by landcover type. Ideally we would have a default landcover layer, for example: https://esa-worldcover.org/, and also give the user the option to add their own layers. This workflow should also work to aggregate quality measures for other zonal characteristics like aspect.

[loreabad6]

Regarding the WOLRDCOVER layer, check for the AWS S3 data bucket that would allow access to the data by an AOI. The example script extracts per country, but that is easily modifiable.

import geopandas as gpd

s3_url_prefix = "https://esa-worldcover.s3.eu-central-1.amazonaws.com"

## ---- LA: This can be skipped
# load natural earth low res shapefile
ne = gpd.read_file(gpd.datasets.get_path("naturalearth_lowres"))

# get AOI geometry (Italy in this case)
country = 'Italy'

## ---- LA: This can be replaced with own AOI geom
geom = ne[ne.name == country].iloc[0].geometry

# load worldcover grid
url = f'{s3_url_prefix}/v100/2020/esa_worldcover_2020_grid.geojson'
grid = gpd.read_file(url)

# get grid tiles intersecting AOI
tiles = grid[grid.intersects(geom)]

# use requests library to download them
import requests
from tqdm.auto import tqdm  # provides a progressbar

for tile in tqdm(tiles.ll_tile):
    url = f"{s3_url_prefix}/v100/2020/map/ESA_WorldCover_10m_2020_v100_{tile}_Map.tif\n"
    r = requests.get(url, allow_redirects=True)
    out_fn = f"ESA_WorldCover_10m_2020_v100_{tile}_Map.tif"
    with open(out_fn, 'wb') as f:
        f.write(r.content)