High-Resolution Land Use Land Cover Dataset for Meteorological Modelling – Part 1: ECOCLIMAP-SG+ an Agreement-Based Dataset

1. Introduction

To estimate parameters required for calculating turbulent, radiative, heat, and moisture fluxes from the surface of the Earth, NWP models need information describing the Earth’s surface. Physiographic databases are used to provide such information. These databases comprise a LULC map, and a complementary set of geophysical datasets such as leaf area index, albedo, tree heights, and lake depths. A LULC map gathers a set of identifiable features from observations (generally remote-sensing) into classes that the map producer desires. The LULC map is a pivotal element of the physiographic database used in NWP models as it triggers the use of different physical parametrisations and biophysical arrays associated with each LULC class. Thus, the classes of the LULC map used in each NWP model need to address the surface physics requirements of the model.

Meteorological organisations use different LULC maps in their NWP models [1]. These maps vary in the number of land-cover classes and grid size (or resolution), but all aim to represent the sub-grid surface heterogeneity of the NWP model correctly. To proceed, the resolution of the LULC map must be significantly smaller than the NWP model grid spacing. Met Éireann, the Irish Meteorological Service, uses the HARMONIE-AROME canonical model configuration of the shared ALADIN-HIRLAM NWP system for short-range operational weather forecasting with a grid spacing of 2.5 km [2]. The LULC map used with HARMONIE-AROME is the latest version of ECOCLIMAP [3,4] called ECOCLIMAP -Second Generation (ECOSG)1 at 300 m resolution [5], less than eight times the HARMONIE-AROME grid spacing. Over the past few years, sub-kilometer resolution NWP experiments have emerged with grid spacings of up to 100 m [6,7,8]. Thus, the current resolution of ECOSG is insufficient, and using a higher-resolution LULC map with ECOSG labelling is necessary. It is essential to retain the ECOSG labelling to avoid the need to rewrite the surface physics code.

Many high-resolution maps exist [9,10,11,12]. As described in [13], these maps can be global, continental [9,10], national [12] or even local [14], and all have limitations due to spatial coverage, spatial accuracy, semantic accuracy and how up-to-date they are. The ECOSG physiographic database has a land cover with 33 labels, which is much less than the 215 labels of ECOCLIMAP-II [4] but much more than many existing LULC. These limitations mean that no existing high-resolution (less than 100 m) datasets with ECOSG labels are currently available globally, or even for Europe.

Machine learning (ML) is well suited to help tackle these limitations. Previous work on improving ECOSG with ML showed potential for producing a LULC map with simpler labels [1] and for extracting building heights to describe urban areas [15] better. However, until now, these ML techniques could not be applied to identify the 33 ECOSG labels to produce a higher resolution map, as no existing reference dataset was available to be used to train and validate the ML model.

The paper, as part one of a two-part publication, introduces an agreement-based method for producing a 60 m resolution reference LULC map with ECOSG labels, leveraging on existing suitable datasets, a map we will call ECOCLIMAP-SG+ or ECOSG+ hereafter. Moreover, the agreement between maps is used to quantify the uncertainty on the land cover labels, which is rarely provided information, although very useful in some applications. In particular, it can be used to build a reference dataset for machine learning easily. This will be done in the second part of this publication which focuses on improving ECOSG+ where its uncertainty is high. Section 2 presents the datasets and methodology used to produce ECOSG+. The results follow in Section 3 while the conclusions are provided in Section 4.

2. Materials and Methods

ECOSG+ is built by mixing a large number of existing land cover maps based on the agreement between them. This section details which land cover maps were used and how the agreement-based mix is performed.

2.1. Material

2.1.1. Primary and Secondary Labels

In the construction and evaluation of ECOSG+, we use two sets of labels: primary labels, ${\mathcal{L}}_1$, and secondary labels, ${\mathcal{L}}_2$. The sets ${\mathcal{L}}_1$ and ${\mathcal{L}}_2$ are such that for each $l_2\in {\mathcal{L}}_2$ there is a unique $l_1\in {\mathcal{L}}_1$. The function for the hierarchical link is denoted h as in Equation (1).

$$h:l_2\in {\mathcal{L}}_2\mapsto h\left(l_2\right)=l_1\in {\mathcal{L}}_1$$

(1)

Figure 1 gives the sets of labels ${\mathcal{L}}_1$ and ${\mathcal{L}}_2$ and the hierarchy between these.

While we aim to create ECOSG+ with the set of secondary labels ${\mathcal{L}}_2$, the method starts with ${\mathcal{L}}_1$ to aid the comparison of maps with heterogeneous semantic coverage. More specifically, some maps focus on primary labels (we refer to these as backbone maps), and some maps focus on some secondary labels (referred to as specialist maps), although the focus is not always exclusive (some maps are both backbone and specialist maps).

2.1.2. Land Cover Maps

The data used in the construction of ECOSG+ is a set of LULC maps. We define a map as a function from a geographical location to a label as in Equation (2).

$$M:x\in \mathcal{D}\mapsto l\in \mathcal{L}$$

(2)

with the following definitions: x is a geographical location, $\mathcal{D}$ is the geographical domain of definition of the map, l is a land cover label, $\mathcal{L}$ is the set of labels for this map (i.e., the semantic domain of definition).

For ECOSG+, the geographical domain of definition is the entire globe, denoted $\mathcal{D}$. However, most of the maps used in the building of ECOSG+ are focused on Europe. Therefore, the evaluation is only done over Europe, and visualisation is restricted to the EURAT domain (longitudes: -32 to 42 degrees, latitudes: 20 to 72 degrees). The semantic domain of definition of ECOSG+ is the set of secondary labels ${\mathcal{L}}_2$.

A total of 43 land cover maps were used in the construction of ECOSG+. They are listed in Table A1–Table A3. For details on the pre-processing and exceptions, we refer to Appendix D.1.

Backbone Maps

We define a backbone map as any map that provides all2 the primary labels ${\mathcal{L}}_1$ (preferably at a high resolution). We denote the total number of backbone maps available as $N_{bb}$. For any $i\in \{1,...,N_{bb}\}$, the i-th backbone map $M_{bb}^i$ is, therefore, a function satisfying Equation (3).

$$M_{bb}^i:x\in {\mathcal{D}}_{bb}^i\mapsto l_1\in {\mathcal{L}}_1$$

(3)

where ${\mathcal{D}}_{bb}^i$ is the geographical domain of the map $M_{bb}^i$. Note that the geographical domain varies from one backbone map to another. In this work, $N_{bb}=16$ backbone maps were used. These are identified as $M_{bb}$ in the Table A1–Table A3.

Specialist Maps

We define a specialist map as any map that provides secondary labels. We denote the total number of specialist maps available as $N_{sp}$. For any $j\in \{1,...,N_{sp}\}$, the j-th specialist map $M_{sp}^j$ satisfies:

$$M_{sp}^j:x\in {\mathcal{D}}_{sp}^j\mapsto l_2\in L_2^j$$

(4)

where ${\mathcal{D}}_{sp}^j$ is the geographical domain of definition and $L_2^j\subset {\mathcal{L}}_2$ is the semantic domain of definition. For backbone maps, the geographical coverage varies from one specialist map to another, as does the semantic domain of definition. In this work, $N_{sp}=33$ specialist maps were used. They are identified as $M_{sp}$ in Table A1–Table A3 which also contain the secondary labels indices that are in their semantic domain of definition.

2.2. Methods

2.2.1. Construction of ECOSG+

This section only describes the main steps of the construction, which are original to this work. For technical details and exceptions, please see Appendix D. The input data used for this method are the backbone maps $\{M_{bb}^1,...,M_{bb}^{N_{bb}}\}$ and the specialist maps $\{M_{sp}^1,...,M_{sp}^{N_{sp}}\}$ with their definition domains, as introduced in Section 2.1.2.

Definition of a Specialist Agreement Score

For any position $x\in \mathcal{D}$ and secondary label $l_2\in {\mathcal{L}}_2$, we define the specialist agreement score, $S_{sp}(x,l_2)$, as in Equation (5).

$$S_{sp}(x,l_2)={\displaystyle \frac{{\sum }_{j=1}^{N_{sp}}⊮\left(x\in {\mathcal{D}}_{sp}^j\wedge l_2\in L_2^j\wedge M_{sp}^j\left(x\right)=l_2\right)}{{\sum }_{j=1}^{N_{sp}}⊮\left(x\in {\mathcal{D}}_{sp}^j\wedge l_2\in L_2^j\right)}}$$

(5)

with $⊮$ the indicator function returning 1 if its argument is true and 0 otherwise, and ∧ the logical "and" operator. Therefore, the score $S_{sp}(x,l_2)$ is the ratio between the number of specialist maps that agree with the label $l_2$ at x versus the number of maps that could provide this information. For example, if we have a position x for which 4 maps can give the secondary label $l_2$="19. Winter C3 crops" (i.e., ${\sum }_{j=1}^{N_{sp}}⊮\left(x\in {\mathcal{D}}_{sp}^j\wedge l_2\in L_2^j\right)=4$) but only 3 maps actually give this label (i.e., ${\sum }_{j=1}^{N_{sp}}⊮\left(x\in {\mathcal{D}}_{sp}^j\wedge l_2\in L_2^j\wedge M_{sp}^j\left(x\right)=l_2\right)=3$), then we have a specialist agreement score of 0.75 (i.e., $S_{sp}(x,l_2)=3/4$). The specialist agreement score ranges between 0 and 1 and reflects the confidence level on the information provided by the specialist maps: the higher the value of $S_{sp}(x,l_2)$, the more confident we are that the label $l_2$ is correct at x. Exceptions to Equation (5) are listed in Appendix D.2.

Refinement of Backbone Maps

For each backbone map, we create a map giving secondary labels instead of primary labels. This process results in a so-called refined map for each backbone map. For any $i\in \{1,...,N_{bb}\}$ we write $M_{rf}^i$, the refined map of the i-th backbone map, as:

$$M_{rf}^i:x\in {\mathcal{D}}_{rf}^i\mapsto \underset{l_2\in {\mathcal{L}}_2,h\left(l_2\right)=M_{bb}^i\left(x\right)}{\mathrm{argmax}}{S}_{sp}(x,l_2)\in {\mathcal{L}}_2$$

(6)

with $S_{sp}$ the specialist agreement score (see Equation (5)), h the hierarchical link between secondary and primary labels (see Equation (1)) and

$${\mathcal{D}}_{rf}^i=\left\{x\in {\mathcal{D}}_{bb}^i:\underset{l_2\in {\mathcal{L}}_2,h\left(l_2\right)=M_{bb}^i\left(x\right)}{max}{S}_{sp}(x,l_2)\ne 0\right\}$$

(7)

The refined map $M_{rf}^i$ returns the secondary label with the highest specialist agreement score, while satisfying the hierarchical link with the primary label given by the backbone map $M_{bb}^i$. For example, let us consider the ESA WorldCover v200 backbone map and a position x where the map gives the primary label "Water bodies" (i.e $M_{bb}\left(x\right)=$"Water bodies"). The associated refined map will return whichever of the following secondary labels $l_2\in$ {"1. Sea and oceans", "2. Lakes", "3. Rivers"} (i.e $h\left(l_2\right)=M_{bb}\left(x\right)$) has the highest $S_{sp}(x,l_2)$.

The domain of definition of a refined map, ${\mathcal{D}}_{rf}^i$, is smaller than that of the corresponding backbone map, ${\mathcal{D}}_{bb}^i$, because refinement is not possible everywhere. When the highest specialist agreement score is zero, there are no specialist maps that provide the relevant secondary labels. Therefore, refinement is impossible.

We define a refined agreement score, $S_{rf}(x,l_2)$, for any position x and secondary label $l_2$ as in Equation (8).

$$S_{rf}(x,l_2)={\displaystyle \frac{{\sum }_{i=1}^{N_{bb}}⊮(x\in {\mathcal{D}}_{rf}^i\wedge M_{rf}^i\left(x\right)=l_2)}{{max}_x\left\{{\sum }_{i=1}^{N_{bb}}⊮(x\in {\mathcal{D}}_{rf}^i)\right\}}}$$

(8)

Therefore, $S_{rf}(x,l_2)$ is the ratio between the number of refined maps that agree with the label $l_2$ at x versus the maximum number of overlapping maps. For example, in our case $N_{bb}=16$, but because the definition domains of some refined maps are mathematically disjoint, no more than 9 refined maps overlap (i.e., ${max}_x{\sum }_{i=1}^{N_{bb}}⊮(x\in {\mathcal{D}}_{rf}^i)=9$). If we have a position x for which 4 refined maps give the secondary label $l_2$="19. Winter C3 crops" (i.e., ${\sum }_{i=1}^{N_{bb}}⊮(x\in {\mathcal{D}}_{rf}^i\wedge M_{rf}^i\left(x\right)=l_2)=4$), then we have a refined agreement score of 0.44 (i.e., $S_{rf}(x,l_2)=4/9$). Note that the denominator is constant. This choice ensures that an area with more available refined maps gets a higher score. Exceptions in the refinement process are listed in Appendix D.3.

Best-Guess Map

The ensemble of refined maps is then used to create a single best-guess map. This is done by taking the label $l_2$ with the best-refined agreement score $S_{rf}(x,l_2)$. The resulting map, $M^{*}$, is defined as follows3:

$$M^{*}:x\in {\mathcal{D}}^{*}\mapsto \underset{l_2\in {\mathcal{L}}_2}{\mathrm{argmax}}{S}_{rf}(x,l_2)\in {\mathcal{L}}_2$$

(9)

with $S_{rf}$ the refined agreement score (see Equation (8)) and

$${\mathcal{D}}^{*}=\left\{x\in \mathcal{D}:\underset{l_2\in {\mathcal{L}}_2}{max}{S}_{rf}(x,l_2)\ne 0\right\}$$

(10)

For the latter, we extend the definition domain of $M^{*}$ to the whole globe $\mathcal{D}$ by inserting the label "0. No data" for all $x\notin {\mathcal{D}}^{*}$. The extended map is also denoted $M^{*}:\mathcal{D}\to {\mathcal{L}}_2$.

Quality Assessment

The quality of $M^{*}\left(x\right)$ depends on the refinement process and the construction of the best-guess map. The uncertainties of these steps are represented by the specialist agreement score $S_{sp}$ and the refined agreement score $S_{rf}$. For any position x, we define the quality score, $S\left(x\right)$, as in Equation (11).

$$S\left(x\right)=\sqrt{S_{rf}(x,M^{*}\left(x\right))S_{sp}(x,M^{*}\left(x\right))}$$

(11)

with the following definitions:

• $M^{*}$ is the best-guess map, defined in Equation (9),
• $S_{rf}$ is the refined agreement score, defined in Equation (8),
• $S_{sp}$ is the specialist agreement score, defined in Equation (5),

The quality score $S\left(x\right)$ is the geometric mean of the uncertainty caused by disagreement in the backbone maps (represented by their refined counterpart) and the uncertainty due to disagreement in the specialist maps. This is the score we use to estimate the uncertainty of the label given by ECOSG+ at any position x. Note that when $M^{*}\left(x\right)$="0. No data", the score $S\left(x\right)=0$ because $S_{sp}(x,\text{"}0.\mathrm{No}\mathrm{data}\text{"})=0$ as defined in Equation (5).

Assembling

The last step in producing ECOSG+ is to assemble the available information, namely: the best-guess map $M^{*}$, the quality score S and the ECOSG map. The assembling of ECOSG+ takes the best-guess map $M^{*}$ where the agreement score S is higher than a threshold $S_{min}$, and uses ECOSG elsewhere. If we denote the map returning the labels of ECOSG (resp. ECOSG+) by $M_{sg}$ (resp. $M_{sg+}$), we have:

$$M_{sg+}\left(x\right)=\left\{\begin{array}{cc}{M}^{*}\left(x\right)\hfill & \mathrm{if}S\left(x\right)>S_{min}\hfill \\ M_{sg}\left(x\right)\hfill & \mathrm{else}\hfill \end{array}\right.$$

(12)

The $S_{min}$ threshold was determined using a histogram of $S\left(x\right)$ for x covering the EURAT domain at 0.1º resolution. Looking for a gap in the histogram using the Otsu method, the value attained for $S_{min}$ was 0.525.

2.2.2. Evaluation of ECOSG+

As for all land cover maps, the evaluation of ECOSG+ is made complicated by the heterogeneity of resolutions, semantic, and geographical coverage. Existing methods for evaluating land cover maps include:

• Comparison to derived measurable quantities. For example, [16] trained a machine learning model to derive the skin temperature from the land cover, and compared the derived skin temperature to the measured skin temperature. This method allows quantitative evaluation of the land cover maps of large domains but requires measured quantities at a comparable resolution over a comparable domain.
• Human validation. For example, LUCAS [17] and CLC+ [18] are validated by human experts. In the case of LUCAS, experts went to designated sites to verify the land cover. In the case of CLC+, experts validated the land cover classes by photo-interpretation. In both cases, human validation requires a sufficient number of trained staff and a carefully designed validation procedure.
• Comparison to trusted land cover maps. For example, [1,19] trained and validated a machine learning model using the CORINE land cover map. In this method, the quality of the evaluation is dependent on the quality of the trusted map. It, therefore, requires the existence of a trusted map of proven quality with an appropriate set of labels, and preferably a higher spatial resolution and greater detail than the map being assessed. Although this is considered less accurate than human validation [20], this method validates every pixel on the trusted map domain.

In our case, human validation was not possible because of limitations in time and expertise. The comparison to derived measurable quantities was not investigated because of the lack of measured quantities at 60 m resolution covering significant parts of Europe. Therefore, we chose to perform a comparison with trusted land cover maps used as references.

Reference Maps

The trusted land cover maps used are the following:

• LUCAS 2022: In-situ data at validated sites over all of Europe translated to primary labels (see Table A7). A translation to secondary labels is not possible.
• NLC 2018: A raster map providing primary labels at 10 m resolution, created by the National Mapping Division of Tailte Éireann (formerly Ordnance Survey of Ireland) in partnership with the Irish Environmental Protection Agency (EPA) and translated to primary labels (see Table A8). A translation to secondary labels is not possible, and the map only covers Ireland.
• ECOSGIMO: A raster map at 25 m (provided at 60 m) resolution providing secondary labels, using national datasets and expert rules. Most of the covers over nature come from a habitat classification map [21] from the Icelandic Institute of National History (IINH) based on the EUNIS classification system. The habitat types were translated to secondary labels for Snow, Water bodies, Bare land, Grassland, Crops and Flooded vegetation. The secondary labels for Forests and Shrubs are based on data from the Icelandic Forest Service - Icelandic Forest Research, Mógilsá. Two maps were used, i.e., a map of native birch forests and shrubs, and a map of afforestation with different coniferous and broadleaf species. The urban local climate zone labels come from the CORINE Land Cover 2018 [22] with a few updates in Reykjavík city. Recent lava fields were added as rocks with data from Icelandic Meteorological Institute and other national institutes in Iceland.

None of the reference maps were included in the construction process, which avoids bias in the evaluation. Three components of ECOSG+ are tested with the three reference maps: accuracy across Europe (LUCAS), accuracy of small-scale features (NLC 2018), accuracy of secondary labels (ECOSGIMO). However, these three components are not evaluated together, and each has limitations (only discrete sites in LUCAS 2022, only the geographical region of Ireland in NLC 2018, and only Icelandic secondary labels in ECOSGIMO).

Comparison Scores

To quantify the similarity between two maps, we compute the confusion matrix on the labels4 for all pixels. Then, we compute the overall accuracy (OA) for an overall comparison, and the F1-score for a per-label comparison.

The confusion matrices shown in Section 3 are normalized row-wise by the number of pixels with the row’s label in the reference map. Therefore, the diagonal of these matrices shows the recall (or producer accuracy) value for each label.

Baseline Maps

Once the comparison scores were computed for ECOSG+ against the reference maps, the values of the scores were compared to the ones obtained using baseline maps. The baseline maps considered are ECOSG, ECOSG+300 and ESA WorldCover v200 ([23], hereafter ESA WorldCover).

• ECOSG is used as a baseline to quantify the improvement on the land cover map currently used in the HARMONIE-AROME NWP model.
• ECOSG+300 refers to ECOSG+ resampled at ECOSG´s native resolution of 300 m. This baseline aims to show whether the improvement is due to the increase in resolution or the correction of some labels.
• ESA WorldCover v200 is one of the most commonly used land cover maps, and therefore makes a good standard.

In addition to the quantitative evaluation described here, a qualitative evaluation is also carried out. It consists of visualizing ECOSG+ and its quality score at different scales and calculating basic statistics about the labels.

3. Results

The results are presented as introduced in Section 2.2.2: first, a qualitative evaluation of ECOSG+ and its quality score (Section 3.1), then a quantitative evaluation. The quantitative evaluation involves testing three components of ECOSG+ separately, marked by their reference datasets: a continental scale evaluation is performed with LUCAS 2022 as a reference (Section 3.2.1), a small scale evaluation is performed with NLC 2018 as a reference (Section 3.2.2), and a secondary label evaluation is performed with ECOSGIMO as a reference (Section 3.2.3).

3.1. Qualitative Evaluation

3.1.1. Overview of the ECOSG+ Map and Its Quality Score

Figure 2 shows an overview of the ECOSG+ map (upper panel) and its associated quality score map (lower panel). It is obtained by taking the nearest neighbour value on a regular longitude–latitude grid (EPSG:4326) covering the EURAT domain at 0.1º. The colormap of the land cover labels is the same as for ECOSG. Only labels actually present on the map are listed in the color bar. The colormap of the scores transitions from red (low values, which means low confidence on the label) to green (high values, which means high confidence). The cut-off value in the colormap is the same as the threshold $S_{min}$ used in assembling the map (see Equation (12)): 0.525. Therefore, the labels of ECOSG+ come from ECOSG where the scores are in red, and they come from the best-guess map $M^{*}$ (see Equation (9)) where the scores are in green.

On the land cover map (upper panel of Figure 2), little can be said at this scale, except that no obviously wrong labels were found. We retrieve the expected main features: credible coastlines, a large area of desert in Northern Africa, and a complex tiling of covers everywhere else.

On the score map (lower panel of Figure 2), the main red areas are over the Atlantic Ocean, Eastern Europe and Mediterranean countries. The low scores over the Atlantic Ocean are due to the common practice of removing data over large sea areas. Consequently, despite the low score values, we are confident that the "1. Sea and oceans" label taken from ECOSG over this area is correct. In Eastern Europe, low scores are expected as fewer datasets representing these areas were included in the method. This is also the case for some Mediterranean countries, such as Turkey, Morocco and Algeria. However, some regions are in red, despite having a good coverage of datasets. For example, southern France and Portugal are depicted in red despite using the national land cover maps there, which suggests disagreement between all the datasets used in these areas.

The main green areas are the coastline, the deserts and most European countries. The coastline is usually well represented in all land cover maps, and some of the maps used in this work have a resolution as high as 10 m. Therefore, this explains the good agreement on the coastline. A large number of maps cover most European countries, explaining the high confidence level in these areas. Surprisingly, the deserts are shown in light green, which suggests a satisfactory confidence level, despite the limited number of datasets available. This result can be explained by the limited number of labels covering the deserts ("4. Bare land" and "5. Bare rock"), which reduces the risk of disagreement compared to areas covered by more labels. Therefore, the numerator in the specialist agreement score (see Equation (5)) is most likely high. This high numerator is divided by a small number of maps available, leading to a large score and, consequently, a high confidence level.

Other noticeable features on the score map are the numerous artefacts due to the manipulation of many maps with heterogeneous projections and boundaries. We can see stitching artefacts (regularly spaced red horizontal lines) and reprojection artefacts (patterns in the North and Arctic Seas). Noticeably, these artefacts are not visible on the land cover map, at least at this scale, which is an encouraging result for the ECOSG+ construction method.

3.1.2. Distribution of Labels

Figure 3 represents the distribution of labels over the EURAT domain for ECOSG+ (outer ring) and ECOSG (inner ring). In addition, the proportion of pixels with scores exceeding the $S_{min}=0.525$ threshold has been estimated as 33.79%, which means that 33.79% of the pixels of ECOSG+ take their label from a source other than ECOSG. The pie chart names only the most common labels (those with 2% coverage or more). The first visible feature is the similar distribution of labels between ECOSG+ and ECOSG, which is expected at this scale. Therefore, although 33.79% of pixels have changed, the distribution of labels is barely modified, which is an element of ECOSG+ validation. The second visible feature is the strong imbalance among the land cover labels. On the one hand, four labels cover more than 81% of the pixels: "1. Sea and oceans" (52%), "4. Bare land" (17%), "19. Winter C3 crops" (7%) and "12. Boreal needleleaf evergreen" (5%). On the other hand, the 15 least dominant labels cover less than 1% of the pixels. In particular, urban areas (LCZs 1 to 10) represent 0.9% of the pixels.

3.1.3. Zoom on a Few Patches

Figure 4 shows examples of land cover patches for several land cover maps (one per column) at several locations (one per row). Each row represents a geographical area and is identified on the left-hand side by a toponym, the country it is in, and the longitude-latitude coordinates of the central point. The first row is the Snaefell Glacier in Iceland. The second row is centred on Nanterre, France, in the north-western part of the Paris urban area. The third row shows the small islands called Kihdinluoto in the south-west of Finland. The fourth row is a rural part of Portugal, around the small town of Pinhel. The fifth and last row is the oasis town of El Menia in the Sahara desert (Algeria). The geographical areas have been chosen to display various landscapes and latitudes within the EURAT domain. All patches are 0.0833º in size, representing approximately 8 km at low latitudes.

Each column represents a different land cover map. The first column is ESA World Cover, one of the backbone maps with global coverage and 10 m resolution, shown here to verify the primary labels. The second column is ECOSG, currently used in NWP models and is the baseline to improve upon. The colormap for the ECOSG labels is the same as in Figure 1, Figure 2 and Figure 3. The third column is ECOSG+, $M_{sg+}$ (see Equation (12)), the final map of this work. The fourth column is the best-guess map, $M^{*}$ (see Equation (9)). The fifth and last column is the ECOSG+ quality score map, S (see Equation (11)), used in the assembling of ECOSG+. The colormap for the score values is the same as in Figure 2, with green indicating where ECOSG+ takes values from $M^{*}$ instead of ECOSG.

Figure 4 illustrates the gain in resolution between ECOSG and ECOSG+ and the label correction (e.g., over Paris –line 2). The best-guess map shows even more pronounced differences with ECOSG very often in agreement with ESA WorldCover, which is good. However, some missing values remain where the score is 0 (deep red, e.g., the crops over Portugal –line 4– or the Algerian lake and town –line 5). This justifies the need for a cut-off value for the score. As ECOSG is taken where the score is below the cut-off value, the areas with the lowest scores also have a lower resolution.

3.2. Quantitative Evaluations

3.2.1. Europe-Wide Evaluation against LUCAS

The ESA WorldCover, ECOCLIMAP-SG map and the new ECOSG+ map were evaluated over the European Union using the LUCAS 2022 dataset [17]. In this comparison, the LUCAS points were expanded to a 60 m radius, and the LUCAS labels were translated according to the C3 classification and Table A7. To rule out resolution changes as a factor for the differences between ECOSG and ECOSG+, ECOSG+ was downsampled to 300 m (ECOSG+300) using the most frequent land cover type in each 300 m square grid.

Figure 5 shows the confusion matrices of ECOSG, ESA WorldCover, ECOSG+ and ECOSG+300. Primary labels of the reference dataset (LUCAS) are on the y-axis and those of the evaluated dataset are on the x-axis. To compensate for label imbalance, the matrices have been normalised row-wise by the number of pixels per label in the reference dataset. The figure also shows each map’s associated overall accuracies (OA). For every dataset, "Forest", "Grassland", and "Crops" display the most off-diagonal spread in values, with a clear confusion between "Bare land" and "Crops". ECOSG+, ECOSG+300 and ESA WorldCover exhibit similar behaviour with higher values on the diagonal for "Grassland", and more "Grassland" versus "Shrubs" confusion than in ECOSG. The OA reflects these two groups with ECOSG+, ECOSG+300 and ESA WorldCover having a value of over 0.5, while ECOSG has a value of 0.42. These observations suggest that the agreement-based method improved the overall representation of primary labels, and these observations are not only due to the resolution increase, as ECOSG+300 exhibits similar behaviour to ECOSG+.

Table 1 displays the F1-score for ECOSG, ESA WorldCover, ECOSG+ and ECOSG+300 with LUCAS 2022 as a reference. ECOSG+ has the highest F1-score for "Bare land", "Snow", and "Flooded vegetation", and outperforms ECOSG for every label. While the impact of the downsampling leads to lower F1-scores, ECOSG+300 outperforms ECOSG for every label except "Shrubs". This shows that ECOSG+ improves the overall representation of ECOSG primary labels over Europe.

3.2.2. Small Scale Feature Evaluation against NLC 2018

In Section 3.2.1, we assessed the overall representation of the primary labels over Europe. In this section, we compare ECOSG+ against a dataset for Ireland to assess whether ECOSG+ represents the high-resolution detail well. Produced by the EPA and Taillte Ireland, NLC 2018 covers the Republic of Ireland. It has two classification levels, with thematic accuracies of 78.5% and 88.7% at Level 2 and Level 1, respectively, while the geometric accuracy (i.e., the area outline) is 87.2%.

To perform the analysis, the NLC 2018 data were rasterised on a 60 m grid and converted to primary labels following Table A8. We note that the NCL 2018 labels are not a perfect fit for primary labels, especially for the Flooded vegetation types. Thus, an assessment of the NLC 2018 primary label conversion is needed. This assessment was performed over Ireland using LUCAS as a reference. The assessment, presented in Appendix E, has determined that NLC 2018 is a suitable reference map for the Republic of Ireland after being transformed into primary labels. However, it has been observed that NLC 2018 classifies "Shrubs" as "Flooded vegetation", meaning that the outcomes for these two primary labels should be viewed with caution.

Figure 6 shows the row-wise normalised confusion matrices for ECOSG+, ECOSG and ESA WorldCover with NLC 2018 as a reference. The dark blue column over grassland indicates the large spread in the grassland classification, which is consistent with the LUCAS confusion matrix (Appendix E). ECOSG+ and ESA WorldCover classify most of the NLC 2018 "Flooded vegetation" as "Grassland", while there is good agreement between ECOSG and NLC 2018 for flooded grassland. As "Flooded vegetation" is overestimated in NLC 2018, the relative underestimation in ECOSG+ and ESA WorldCover suggests a better representation of "Flooded vegetation" than in ECOSG. The forest land cover in ECOSG+ and ESA WorldCover corresponds well to NLC 2018 forest, while ECOSG overestimates "Grassland".

Table 2 shows the F1-scores over the Republic of Ireland for ECOSG+, ECOSG and ESA WorldCover, using NLC 2018 as a reference. ECOSG+ has the highest F1-score, ahead of ESA WorldCover, for "Forests", "Grassland", "Crops" and "Urban" areas. Meanwhile, ECOSG is a better fit for "Bare land" and "Flooded vegetation", which were the least accurate in NLC 2018 compared to LUCAS (Appendix E). This suggests that the agreement-based method correctly captures most primary labels over Ireland, even when compared to a higher resolution reference. The quality of the "Bare land", "Shrubs" and "Flooded vegetation" labels remains questionable in all tested maps.

3.2.3. Secondary Label Evaluation against ECOSGIMO

In the previous sections, we evaluated ECOSG+ primary labels over Europe and its level of detail over the Republic of Ireland. This comparison is incomplete, as this work aims to create a reference dataset with the 33 secondary labels. An evaluation at a secondary label level is necessary. In this section, we will compare ECOSG+ secondary labels with a high-resolution version of ECOSG called ECOSGIMO. ECOSGIMO is the only high-resolution version of ECOSG available. However, ECOSGIMO does not have an accuracy assessment and its coverage over Iceland is outside the LUCAS domain, or that of any known in-situ datasets. Hence, we cannot verify the accuracy of ECOSGIMO. Rather than providing an absolute assessment of the ECOSG+ secondary labels, this section aims to demonstrate how closely the ECOSG+ agreement-based method can approximate a more traditional technique.

Figure 7 shows the row-wise normalised confusion matrices for ECOSG+ and ECOSG with ECOSGIMO as a reference. For visibility purposes, labels not represented over Iceland were removed from Figure 7. Therefore, we can only evaluate 20 labels. The overall accuracies are similar and high compared to the previous section due to the imbalance of label distribution in the Icelandic domain. Indeed, 48% of the domain is covered by "1. Sea and oceans" and "6. Permanent snow", and the dark blue pixel on the diagonal for these labels reveals a strong agreement of ECOSG+ or ECOSG with ECOSGIMO, which increases the proportion of correctly represented labels.

The dashed squares in Figure 7 indicate the primary labels, and the diagonal in the "Water bodies" square indicates that Water bodies in both ECOSG+ and ECOSG agree with ECOSGIMO at a primary level. Still, at a secondary level, rivers in ECOSG+ are better represented compared to in ECOSG. In the "Bare land" square, ECOSG+ agrees more at a primary level but is more confused at a secondary level. For "Forests", "Shrubs" and "Grassland", both maps overestimate "Grassland", but ECOSG+ preferentially returns "16. Boreal grassland" while ECOSG preferentially returns "17. Temperate grassland" and the latter is the most frequent in ECOSGIMO, which explains a significant part of the overall accuracy gap between ECOSG and ECOSG+. For "Flooded vegetation", ECOSG is in good agreement with ECOSGIMO, while ECOSG+ overestimates "Grassland". For "Urban", ECOSG confuses sparsely built labels with "Grassland" compared to ECOSG+. The observed underestimation of sparsely built areas in ECOSG is consistent with previous studies, and ECOSG+ seems to resolve this issue partially.

Table 3 and Table 4 show the F1-score of ECOSG and ECOSG+ over Iceland for primary and secondary labels with ECOSGIMO as a reference. ESA WorldCover was added at the primary level for information. ECOSG+ has the highest score at the primary level, except for "Forest", "Shrubs", and "Flooded vegetation". Interestingly, ECOSG+ has the worst score for "Shrubs", which is also a secondary label, consistent with other analyses showing that "Shrubs" are still problematic in ECOSG+. At the secondary level, ECOSG+ scores better than ECOSG for 14 labels, while ECOSG has a better score for 4. Neither map identifies "13. Boreal needleleaf deciduous" and "28. LCZ5: open midrise" explaining the NaN values.

Among the labels where ECOSG+ is better than ECOSG in terms of the F1-score, for 8 of the labels the gap between the two scores is more than 0.1. Of the 4 labels, ECOSG has a better F1-score than ECOSG+. The gap between the two scores is more than 0.1 for two labels: "4. Bare land", "17. Temperate grassland". For "4. Bare land ", the low F1-score is due to a distribution issue between the "Bare land" primary label and the two secondary labels "4. bare land" and "5. bare rock". This issue is due to the limited number of specialist maps distinguishing "4. Bare land" and "5. Bare rock". The bioclimatic distribution between the grassland primary label and the "16. Boreal grassland" and "17. Temperate grassland" secondary labels seems to be off in Figure 7 explaining the relatively low boreal grassland F1-score. This shows the limits of solely relying on [24] for the bioclimatic classification (see Appendix D.3 for more detail).

4. Conclusions

In this work, we presented an agreement-based method for producing ECOSG+, an LULC map at 60 m with ECOSG labels. ECOSG+ results from the assembly of ECOSG and a so-called best-guess map $M^{*}$ representing the most probable ECOSG labels combining information from many suitable datasets. The quality score S represents the confidence level of $M^{*}$ and defines the ECOSG+ assembly rules. A threshold value of 0.525 was defined for $S_{min}$ using the Otsu method, to define pixels where $M^{*}$ replaces ECOSG in ECOSG+.

Although 33.79% of the resulting ECOSG+ map pixels are sourced from $M^{*}$, ECOSG+ keeps a similar label distribution to the original ECOSG. A qualitative evaluation of ECOSG+ shows the gain in resolution between ECOSG+ and ECOSG and the label correction. The qualitative evaluation also exhibits the limits of $M^{*}$, with some missing values, and therefore the need for the cut-off value $S_{min}$ despite losing higher resolution information.

We performed the first quantitative evaluation of ECOSG, in parallel with ECOSG+, against LUCAS 2018, NLC 2018 and ECOSGIMO. The evaluation against LUCAS revealed the superiority of ECOSG+ across Europe for every primary label at the LUCAS sites. The downsampling of ECOSG+ to 300 m (ECOSG+300) showed that the improvements between ECOSG+ and ECOSG are not only due to resolution. The most improved labels are "Forest" and "Grassland" while "Bare land" and "Shrubs" have small improvements and are the least well-represented labels.

The evaluation of small-scale features against NLC 2018 revealed that ECOSG+ improves upon ECOSG for every label except "Bare land" and "Flooded vegetation", for which the NLC 2018 translation to primary labels has been identified as being less reliable. "Shrubs" were confirmed as poorly represented in every map.

The evaluation of secondary labels against ECOSGIMO showed that ECOSG+ is superior for 14 of the 20 labels represented in Iceland. However, two secondary labels were inferior: "4. Bare land", "17. Temperate grassland". While their respective primary labels, "Bare land" and "Grassland", were superior in ECOSG+, the lack of a specialist dataset for "4. Bare land" and the sole reliance on the [24] distribution for the bioclimatic distribution showed the limits of the method.

The evaluation showed limitations in representing "Bare land", "Shrubs", and "Flooded vegetation". Therefore, implementing new specialist maps covering these labels would greatly benefit these aspects of ECOSG+. The classification of "Urban" areas into local climate zones needs further evaluation. Due to the limited number of datasets explicitly representing local climate zones, the method could be improved by looking at the agreement based on building height and density, similar to what was done for trees (see Appendix D.3). The current study focused on the EURAT domain. However, the method to create ECOSG+ is flexible enough to implement other datasets and to expand ECOSG+ and its verification domain beyond EURAT.

This study showed that ECOSG+ outperforms ECOSG at the primary level regarding overall accuracy and F1-score per label. Still, it would have been desirable to use a larger domain with more secondary labels to ensure the quality of all 33 labels. One way to assess this product could be to compare it to derived measurable quantities.

However, the aim of ECOSG+ is not only to provide an accurate land cover map but also a quality score on the given land cover. Such information is rarely provided along with LULC maps although it enables further use of the LULC information. For example, the quality score can be used to build a trustworthy reference dataset in machine learning applications. This work is done in [25], where ECOSG+ trains the AI model to replace the areas with the lowest S scores.

As land cover is only one aspect of the ECOCLIMAP-SG physiography database, additional work will be needed to ensure the complementarity between land cover and geophysical parameters such as leaf area index, albedo, tree heights, and lake depths.

Appendix D.1. Exceptions in the Land Cover Maps

Pre-processing is applied and consists of uploading all required maps onto GEE and reprojecting these onto the GEE default grid (see [61] for details). Most of the data are hosted in the GEE data catalogue [62] and the so-called "awesome" GEE community catalogue [63]. Additional datasets not in these catalogues, such as CLC+, were manually uploaded to GEE. Another reprojection is made when exporting the ECOSG+ map and its quality estimation in EPSG:4326 at 0.000539° resolution (approximately 60 m). These operations are dealt with by GEE commands that are beyond the scope of this paper. We refer the reader to Wu [61], Wu et al. [64] for comprehensive documentation on the Python API and the geemap python package.

Some particular cases in Table A1–Table A4 are noteworthy:

• The Geoclimate dataset consists of a map of LCZ obtained by running the Geoclimate tool [45] on the main European urban areas.
• The Copernicus Imperviousness HRL does not provide secondary labels but distinguishes (secondary labels "4. Bare land", "5. Bare rocks" ) and concrete runways ("31. LCZ8: large low-rise") (see Appendix D.2). We extract the artificial imperviousness density, denoted $\varphi \left(x\right)$, from this dataset.
• ESAGHSurban is a combination of ESA WorldCover v200 and the GHS built-up surface (GHS-BUILT-S) converted to primary labels, where GHS-BUILT-S was resampled to the target grid, and missing urban areas have been added from ESA WorldCover following Table A5
• MACATECOSG is a merge of the Portuguese DGterritorio MACAT and COSC 2020 [51] maps. Using the rules in Table A6 this merge locates the labels "19. Winter C3 Crops", "20. Summer C3 Crops", and "21. C4 crops" over Portugal.

Appendix D.2. Exceptions in the Specialist Agreement Score

Exceptions to Equation (5) are made in the following cases:

• $l_2=$"0. No data". The specialist agreement score is set to 0 everywhere for this label.
• Null denominator: ${\sum }_{j=1}^{N_{sp}}⊮\left(x\in {\mathcal{D}}_{sp}^j\wedge l_2\in L_2^j\right)=0$. When no specialist map provides the label $l_2$ at x, the specialist agreement score is set to 0.
• $l_2\in${"4. Bare land", "5. Bare rocks"}. Preliminary experiments showed that confusion often occurs between sand, rocks (secondary labels "4. Bare land", "5. Bare rocks" ) and concrete runways ("31. LCZ8: large low-rise"), which can be disambiguated thanks to the artificial imperviousness density. Therefore, maps providing the labels "4. Bare land " and "5. Bare rocks" see their score penalised by the imperviousness density (sand and rocks with high imperviousness is likely to be wrong).

Therefore, a more general (but more complicated) formula for the specialist agreement score is:

$$S_{sp}(x,l_2)=\alpha (x,l_2)\left(\sum _{j=1}^{N_{sp}}⊮\left(x\in {\mathcal{D}}_{sp}^j\wedge l_2\in L_2^j\wedge M_{sp}^j\left(x\right)=l_2\right)\right)+\beta (x,l_2)$$

(A1)

with the following coefficients:

$$\alpha (x,l_2)=\left\{\begin{array}{cc}1/N_{sp}^{\prime }(x,l_2)\hfill & \mathrm{if}{N}_{sp}^{\prime }(x,l_2)>0\mathrm{and}{l}_2\ne \text{"}0.\mathrm{no}\mathrm{data}\text{"}\hfill \\ 0\hfill & \mathrm{else}\hfill \end{array}\right.$$

(A2)

$$\beta (x,l_2)=\left\{\begin{array}{cc}-\varphi \left(x\right)\hfill & \mathrm{if}{l}_2=\text{"}4.\mathrm{bare}\mathrm{land}\text{"},\text{"}5.\mathrm{bare}\mathrm{rocks}\text{"}\hfill \\ 0\hfill & \mathrm{else}\hfill \end{array}\right.$$

(A3)

where $\varphi \left(x\right)$ is the artificial imperviousness density (ranging between 0 and 1, with 1 being totally human-made impervious ground).

Appendix D.3. Exceptions in the Refinement Process

The refinement process as described in Section 2.2.1 suffers three exceptions: a joint maxima in score, the labels with bioclimatic classification and the forest primary label.

Joint Maximum in Score

In Equation (6), in the case of a joint maximum in the specialist agreement score, the argmax returns the one with the lowest label number. For example, if we have a position x, a backbone map $M_{bb}^i$ giving $M_{bb}^i\left(x\right)=$"Crops", and a joint maximum in the specialist agreement score (e.g $S_{sp}(x,$"21. C4 crops"$)=S_{sp}(x,$"20. Summer C3 crops"$)>S_{sp}(x,$"19. Winter C3 crops")), then the refined map $M_{rf}^i\left(x\right)$ will return "20. Summer C3 crops" because it has the lowest label number in the joint maximum. This default behavior is arbitrary, as no better solution was found for now.

Bioclimatic Classification

Bioclimatic classification (i.e "boreal", "temperate", or "tropical") is used in the "Forest" and "Grassland" primary labels. However, few datasets provide the bioclimatic class. Therefore, the secondary labels that differ only by their bioclimatic class are distinguished following the Beck et al. [24] bioclimatic map and the score is calculated independently of its bioclimatic class. For example, "8. Temperate broadleaf deciduous" and "9. Tropical broadleaf deciduous" are two types of "broadleaf deciduous". They will have the same score values: $S_{sp}(x,$"8. Temperate broadleaf deciduous"$)=S_{sp}(x,$"9. Tropical broadleaf deciduous") for all x.

The Forest Primary Label

The "Forest" primary label includes secondary labels having a more regular structure than the others primary labels. Indeed, the seven classes of forest are a combination of

• The bioclimatic classification (i.e "boreal", "temperate", or "tropical")
• The dominant leaf type (i.e "broadleaf", "needleleaf")
• The leaf cycle (i.e "deciduous", "evergreen")

Therefore, the secondary labels that fall under the "Forest" primary label are obtained using a deeper hierarchy of labels than other primary labels. The bioclimatic class is given by Beck et al. [24], as explained in the previous paragraph. The dominant leaf type and the leaf cycle are jointly taken from datasets specific to these information. The refinement score $S_{rf}$ (see Equation (8)) is then calculated ignoring the bioclimatic class.