Aerial Data Analysis for Integration Into a Green Cadastre—An Example From Aarhus, Denmark

Fostering urban resilience and adaptation to climate change pose new demands on the knowledge of land use and land cover (LULC) in heterogeneous urban spaces. High-resolution urban mapping is a valuable tool, which serves to map detailed categories. Such semantic data are integrated in national and regional administration as public goods. In the light of many countries around the globe making their data publicly available, we present a method to map urban areas based on multitemporal orthophotos and LiDAR-derived digital surface model, and extract information about vegetation in an automated processing chain. This approach is threshold driven and relies on an automatic generation of spectral thresholds and existing real-world-based classifications. We included cadastral data to add land-use information for specific categories, such as agricultural land use and to assess the product's accuracy. Adding these data creates an LULC product and makes a seamless integration into urban planning routines possible. The results of the study provide a detailed LULC map for the municipality of Aarhus in 2015 with a spatial resolution of 20 cm and ten thematic classes. Depending on the reference data, we achieved thematic overall accuracies of 34% and 64% using a polygon-based approach. Our study has found that utilizing both multitemporal orthophotos and elevation data can enhance the LC mapping of urban landscapes. The methodology could be transferred to other areas in Denmark or to countries providing similar datasets, and lends itself to a repeatable LULC mapping with minimal user interaction.


I. INTRODUCTION
I N THE light of global climate change, more and more countries worldwide have to evaluate and adapt their available ecosystem services (ES) for the human habitat, especially in urban areas where people live at highest densities [1]. The existence of and access to urban green spaces (UGS) are one of the most valued remedies to respond to environmental pressures from heat, noise, and air pollution, may they be driven by climate change, and/ or by urbanization processes [2], [3].
Multiple studies have shown the positive influence of UGS on mental and physical health, be it through opportunities to exercise, better air quality, mitigating urban heating, shielding from noise, meeting places, or aesthetics in greenery [4], [5], [6], [7]. Assessment of biodiversity in cities, with a focus on structural biodiversity and flood mitigation are even more fields, in which remote sensing (RS) can provide additional information and help with the planning of adaption strategies [8], [9].
At the same time, the evaluation of urban pressures was mostly based on datasets, which could not provide a robust basis for future planning on a neighborhood level, but rather at a masterplan level with regards to their spatial extent and resolution as well as their thematic detail [10]. The pressure to improve data quality has increased. Very high-resolution (VHR) imagery has been a basis for detailed land-cover (LC) datasets with meaningful land-use and land-cover (LULC) categories for more than two decades, and a small minimum mapping unit (MMU) is, thus, in order, as this is required for decision making on a city level. These more detailed LULC products can either aim at aspects of climate change resilience, biodiversity, or the wellbeing of citizens [11].
Over the last decades, not only camera techniques and sensor systems have improved greatly but the analytical tools to process enhanced data have made a big leap forward. With the advent of object-based image analysis (OBIA) in RS, many studies have utilized digital orthophotos (DOP) to classify urban LC with a higher spatial resolution. In addition, the linkage with auxiliary light detection and ranging (LiDAR) datasets has had a decisive influence on the potential of 3-D analysis, not least in improving the accuracy of the resulting mapping products [12]. When analyzing VHR data, the OBIA approach has the advantage of using objects rather than single pixels for the analysis. By means of segmentation, homogeneous pixels are grouped into objects, which can then be classified. Aside from the resulting products showing less spatial and thematic noise, the objects themselves can yield information about their area, circumference, texture, and relation to neighboring objects among other traits [13], [14]. With OBIA, many of the inherent problems of urban RS can be overcome. Among these are the heterogeneous distribution of spectral characteristics, occlusion and shading, as well as features being represented by multiple pixels [11].
With the expansion of open data policies, many European countries have opened their archives to allow access to DOPs and even digital surface models (DSMs) and digital terrain models (DTM) [15]. DOPs are captured in regular sequences, often along with LiDAR data, or they are used to create DSMs through This work is licensed under a Creative Commons Attribution 4.0 License. For more information, see https://creativecommons.org/licenses/by/4.0/ photogrammetric algorithms. These datasets, however, are not necessarily usable in the scientific context of urban and regional planning, as the access is not always user friendly. Furthermore, the lack of machine-readable, consistent, and searchable metadata can be another barrier to data acquisition and analysis [16], [17].
Another drawback of aerial RS imageries, when compared with other optical RS data, especially Earth observation satellites, is the lack of radiometric, spectral, and temporal resolutions. While the latter is not of importance in the context of this study, satellites, such as WorldView 3, feature more spectral bands at a higher radiometric resolution. Apart from the three common visible red, green, and blue (RGB) bands, they also feature information about the red edge and two near infrared (NIR) regions of the spectrum at a bit depth of 12 bits. Other satellites make use of a higher spectral resolution as well, such as the Pleiades and SkySat family, in combination with fixed revisit times. In contrast, DOPs are more often supplied in annual or longer time intervals, in the case of Denmark with alternating leaf-off and leaf-on situations. Their spectral bands usually cover the visible RGB range and additional NIR information. Nevertheless, the often limited radiometric and spectral resolutions of DOPs do not necessarily hinder use in research, since the spatial resolution of RS datasets is deemed more important by analysts than spectral information [18]. The added information from extended radiometric resolution does not necessarily improve the information about smaller objects in the heterogeneous urban landscape.
In the Danish context, existing studies exploiting DOPs and/or DSM/DTM data with an OBIA approach include landslide mapping [19], updating cadastre maps [20], and monitoring of shrub and tree cover [12]. Furthermore, two recent examples show the mapping of individual trees [21] and the classification of sealed surfaces [22]. However, to the best of authors' knowledge, the potential of mapping the diverse LULC in urban areas by means of OBIA has not yet been studied.
Considering these findings, this study aims at developing a workflow for an automated LULC product based on DOPs from different time points and a normalized digital surface model (nDSM). The main goals are the following: 1) extract high-resolution LULC information for local decision-making based on available datasets; 2) examine and apply a fitting accuracy assessment; 3) create a basis for the analysis of green space distribution, composition, and accessibility, and filling gaps in existing cadastral data; 4) deliver a refined-scale dataset for ecosystem-service modeling.

II. APPLIED RS IN CADASTRAL PRACTICE
One of the reasons why this study was undertaken is a lack of high-resolution LC data for urban areas covering both public and private areas. While many cities provide cadastral data with extensive information about the built environment, the totality of green infrastructure is often not mapped. One part, which is usually covered in Europe, is trees on public ground, as a database on them is important for both management and liability issues [23]. Studies from the U.K. and New Zealand show that private gardens on the other hand comprise as much as 36% of a given city's area, with a green coverage of up to 50%, and are not part of cadastral elicitation [24], [25]. A study for Padua, Italy, found that as much as 80% of UGS are located on private ground, concluding that this can be a major influence in assessing ES [26]. In addition, green spaces within business areas are also often not captured in such surveys. The datasets from the aforementioned studies on urban green infrastructure, however, are not well distributed and focus on specific countries. These areas of urban green provide a significant amount of ES to the city and need to be taken into account when planning green infrastructure in cities.
In addition, this study aims to implement a per polygon approach for its accuracy assessment. A conventional confusion matrix based on the thematic classification of single pixels may not be sufficient when working with OBIA techniques [27], [28]. As it only takes into account the true or false assignment of the thematic classes per pixel, and pixel-based mapping with VHR data tends to introduce noise to its products, such confusion matrices become less meaningful. To pay tribute to the classified objects, many different studies have developed frameworks, which encompass, but are not limited to, the extent, position, shape of polygon, and the evaluation thereof. These methodological approaches have also been compared and reviewed [28], out of which this article opted to go with the STEP approach [29]. This process evaluates the Shape, Theme, Edge, and Position of a given polygon based on a reference dataset. To assess the correctness of a given shape, the normalized perimeter index is calculated as an equal area circle. Thematic accuracy is measured as a percentage of intersection, as is the reference object's boundary agreement with the classified object's boundary for the edge metric. Finally, the positional accuracy is expressed as the distance between corresponding centroids, normalized by the diameter of a combined area circle [29]. Fig. 1 illustrates the different metrics used to evaluate the accuracy of a given classification.
Overall accuracy values for this evaluation are reported to be lower than their pixel-based counterparts, may it be linked to the lack of high-resolution reference data or oversegmentation, i.e., the splitting of real-world objects into multiple segments [30], [31], [32].
Since the intended uses of the resulting dataset are manifold, we also opted to create a hierarchical classification structure, which is based on the relation of the LC classes and the actual processing of the data. This translates into multiple categories of LC, with different levels of thematic aggregation, and with additional information, such as height or area about an individual object stored in separate parameters.

III. STUDY AREA AND DATASETS
The study area consists of the municipality of Aarhus, located in the eastern part of the Jutland peninsular in Denmark, as seen in Fig. 2. The municipality is located at the bay of Aarhus, part of the Kattegat, and covers an area of 468 km 2 it consists of the city of Aarhus and its rural hinterland. With a total population of 282 910 [33], Aarhus is the second largest city in Denmark after Copenhagen. In 2021, the total population within Aarhus municipality was 352 751, and it is expected to increase to almost 400 000 in 2032 [34], posing challenges to the future provision of ES. The landscape was formed under the last glaciation and is characterized by moraine sediments with elevations between 0 and 128 m a.s.l. As the city center is located in a valley bottom it is consequently prone to flooding, particularly under the expected future increase of extreme precipitation events [33], [34].
The datasets used in this study are listed in Table I. They consist of DOPs for multiple years and phenological stages, DTM and DSM data, as well as cadastral maps and ancillary data containing information about vegetation, building footprints, and agricultural parcels.
The DOPs cover the red, green, and blue range separately, with an additional channel in the NIR infrared spectrum at a spatial resolution of 12 to 20 cm. The scenes for 2014 and 2018 are from leaf-off season in early spring, whereas the scenes for 2015 and 2019 are from leaf-on season in summer [35], [36], [37], [38]. The DSM and the DTM are from leaf-off season in spring 2015 [39].
For the accuracy assessment, various categorical spatial datasets were applied. From the Danish topographical database [40] the polygon, polyline and point data for building footprints, water bodies, trees, groups of trees, hedgerows, and shrub vegetation were taken into account [41]. A tree database containing information on species for around 28 000 trees in public spaces was used as well [42]. Four datasets were applied for postprocessing. From the Danish topographical database [40], the coastline was applied to mask out areas with open sea and water bodies were applied to refine the mapping of the water category. The map of agricultural field parcels for 2018 [43] was applied to identify agricultural land use types. The map contains information on 301 crop types, which, following the method described in Levin [44], were aggregated into three major categories: agriculture, intensive, temporary, which includes annual crops and rotational grassland, agriculture, intensive, permanent, which includes perennial crops, such as energy crops, fruit plantations, and Christmas trees, and agriculture, extensive, including permanent grassland and other extensively managed land. The map of protected habitat types from 2018 was applied for the refinement of the delineation of the grass category [45].

IV. METHODOLOGY
As a first step, the nDSM was generated by subtracting the DTM from the DSM. Using a bilinear resampling method, the nDSM was converted to a spatial resolution of 20 cm, matching the spatial resolution of the DOPs. The DOP from 2019 was resampled to 20 cm as well.
This study used a supervised classification approach, focusing on a minimal input from the user. Fig. 3 describes the schematic workflow of the process designed in Trimble eCognition, version 9.0.0. Limitations in the spectral resolution of the DOPs demanded the use of well-established and simple ratios for feature classification. In summary, the process combines a height-based segmentation with spectral information for further classification, while establishing a robust framework for threshold detection along multiple time intervals with minimum user interaction.
As a first processing step, an initial segmentation of the project was undertaken, classifying the resulting objects purely on height derived from the nDSM. This nDSM is computed by subtracting the DTM from the DSM to derive height information without the measurements including a.s.l. The resulting classification includes the four preliminary classes of the height category, with thresholds being used at 0.3, 1, and 2 m. Those values relate to the mean height information of the objects. This initial classification was used to decrease computing time by looking at different height levels individually. The classes in use are depicted in Table II, along with their hierarchical structure.
Spectral information was then integrated into the workflow by first calculating individual layers for all four years containing information about the NDVI and NDWI, as well as the height based on the nDSM.  By selecting this approach, all LC classes used in this study can be described by their height, a grouping into vegetated and nonvegetated states, and morphological as well as relational parameters.
The four different height levels were first evaluated independently, classifying objects within each domain separately. The classes used in the process were defined by their height, thus LC categories, such as grass, water, and mineral surface, were only assigned to features showing less than 0.3 m of elevation. Shrubs, hedges, flat buildings, and semisubterranean buildings and structures were only classified up to 2 m of mean height above ground, while houses and trees were solely present in the highest objects above 2 m. This mean value reflects findings from previous studies, which set the thresholds for shrubs and trees at 1.8 m [12], [46]. Other values do, however, exist [47], [48]. To account for LC, which spanned multiple height categories, after this classification, objects were merged based on similarity in height and spectral signatures. Buildings were classified based on the nDSM and differences in minimum and maximum pixel values as well as their slope information for improved edge detection. The tree-crown delineation used as local maximum (Max) and local minimum (Min) approach identifies the highest point in tree crowns and the lowest points or boundaries between adjacent trees. These seed points were then grown downhill, with a limit at the local Min points. Afterward, individual objects were adjusted based on their shapes so as to counteract oversegmentation. The identification of evergreen and deciduous trees was achieved by using the bitemporal band ratio with the leaf-on/leaf-off situation in spring 2015 and 2019 and leaf-on in 2014 and 2018; see (4a) and (4b) [49]. All objects of the class tree were evaluated based on the number of neighboring trees and the number of trees in close proximity to classify them as solitary trees or trees in a forest or a stand.
Iterative merging and splitting of objects were finally translated into different vegetation heights, as to classify the structural biodiversity.
Since the DOPs used in this study were nontrue orthophotos, i.e., not all pixels were of a perfect nadir view, different years and different areas in each image presented changing viewing angles, which occluded different parts of the scene. This shift could be corrected in elevated LC categories using the nDSM by classifying the elevated image objects based on areas included in the nDSM footprint. In contrast, shifting viewing angles and differing shadow situations along with shifting water body extents and different vegetation states for low-standing flora made the use of agreement maps necessary.
Based on all four mosaicked DOP datasets, the classification was assigned using the agreement between NDVI and NDWI based on four different years, and objects were classified binary into water/not water and vegetation/not vegetation based on the division by the automatic thresholds. Agreement levels of > 2, that is, the threshold criteria fulfillment in more than two years, were subdivided accordingly. For objects with no agreement, classes were solely assigned on the threshold (NDVI and NDWI) in 2014. In addition, morphological criteria along with seed points for areas with steady spectral response were used to determine the prevailing LC. This greatly reduced the influence of shadow as well as different states in the phenological cycle of the plants.
The refinement of the classified water bodies was carried out using all values of the NDWI range over four years and the respective NGI values, thus enhancing the contrast for water and, simultaneously, eliminating artifacts based on shadow. As the viewing and illumination angles differed between the four DOPs, shaded areas mistakenly classified as water could be assigned to another LC class by taking into account the other three DOPs. Based on this, seed points for water bodies were identified. They were then grown into adjacent objects and subsequently checked for coherence. Water areas in direct neighborhood to elevated structures, such as buildings, were removed as they were most likely misclassified as shadows. Bigger objects were reviewed for the existence of shadows based on their proximity to higher nDSM values and the differences in RGB values.
After the classification of water bodies, the three remaining categories-grass, mineral surface, and bare soil-were introduced. Objects with a stable NDVI and an elevation of less than 0.3 m were assigned as LC class grass, and all other objects were divided into the two categories of mineral surface and bare soil. Mineral surface encompasses pavement and tarmac with an elevation of less than 0.3 m. Bare soil was separated from mineral surfaces using NIR brightness levels and the REI ratio in (5).
The remaining objects with an elevation of more than 0.3 m were subsequently classified into lower buildings and shrubs, using their respective NDVI values. Objects of the class shrub bordering trees were then considered as a candidate class, and using the difference in height and NIR signatures were conditionally merged with trees to achieve a more precise crown delineation.
In a final step, smaller objects with less than 10 m 2 were examined and either merged with neighboring objects based on elevation and ratio similarities or were kept. This step also featured morphological operations, where shrubs were subdivided into linear and round features, representing managed linear hedgerows and other untrimmed shrubs. All trees with less than three tree objects in their immediate neighborhood were additionally assigned a value for single trees. The remaining trees were given the attribute of belonging to a stand or forest. This additional information is supposed to distinguish street trees from others once the cadastral information is linked to the LC product.
The first ratio used in this study was the NDVI, defined as represents the normalized ratio of the NIR and red spectral signature where NIR describes the NIR channel and R the red channel [50]. It has been widely used in RS since its first implementation and has proven to distinguish vegetation from other LC features [51], [52]. To account for the water bodies in Aarhus, the NDWI, defined as was used, where G denotes the green channel [53]. Normalizing the ratio of the NIR and green channel, it provides a well-established methodology to separate water from other land cover features. In addition, the normalized green index (NGI) was adapted to the four bands of the DOPs by removing the red edge channel from the equation resulting in [54] Here, B represents the blue channel. For the identification of deciduous and evergreen trees, the bitemporal band ratio (BTBR) was selected [49]. The original equation is described using with on and off denoting the leaf-on and leaf-off phenological stages when the DOPs were captured. However, to exploit the NIR capabilities of the sensor used in capturing the aerial images, the original ratio was adapted to [49] This new ratio exploits the reflections in the NIR spectrum, making it akin to the NDVI. Finally, to account for impervious surfaces such as tarmac, pavement and concrete, the road extraction index (REI) was used [55]. Its main characteristic is a better performance on shaded LC, which can cause a misclassification and confusion with water bodies.
By combining all input data and the delineated equations, the image processing was set up and a high-resolution LC (HRLC) map was generated.
The results were then exported with their attributes in a vector format for further postprocessing. This additional postprocessing was needed to find inconsistencies in object geometries as they had to be corrected for the accuracy assessment. In addition, the datasets containing protected, agriculture and marine areas were used to reclassify intersecting areas of the map. This reclassification was based on the originally processed data, so that only thematically overlapping classes were changed. The process relied, thus, heavily on the signature of the nDSM data, where classes, such as, for instance, trees and shrubs, could not be reclassified as agriculture. With this, we achieved the retaining of smaller distinct areas, such as woodlots or hedgerows within agricultural fields. These reclassified objects were, however, excluded from the accuracy assessment, so as not to introduce a bias to the results.
The reference datasets discussed in Table I were gradually reclassified and merged into a single vector dataset, which served as a reference dataset for the accuracy assessment. To pay tribute to the presence or absence of certain LC classes in the thematic classification, the reference dataset had to be structured carefully. Thus, vector datasets had, according to their reported accuracies, higher priority in areas where they overlapped with other, comparable classes, and elevated LC features were gradually introduced, starting from the lowest to the highest.
For the sampling design, publicly not accessible areas of the Aarhus municipality were excluded, as they lacked sufficient high-resolution LC reference data. The remaining parts of the municipality were covered by cadastral data showing buildings, roads, water bodies, and single street trees as well as shrubs. However, since some of the data were only available as line or point features, it had to be buffered in order to be spatially explicit. These data constraints, however, limit the accuracy of the generated reference dataset. Simple extents were assumed for shrubs and trees, ranging from 1 to 3 m in radius and width. Furthermore, areas with artificial grass were provided through visual interpretation of different sports fields around Aarhus.
A second reference dataset based on manual digitization within the study area was also used for comparison. It is based on 562 sample sites, ranging from 30 to 135 samples per class. The sampling sites were chosen based on a random stratified sampling to provide sufficient insight into the accuracy of each individual LC class.
The accuracy assessment used in this study differs necessarily from medium-or low-resolution pixel-based approaches as the geometries of the objects and their positions had to be evaluated along with the thematic accuracy. As discussed in the outline of this study, the work by Lizarazo [29] was implemented for the evaluation of the results. The STEP procedure is readily implemented under current QGIS versions (version 3.18.2) and was used in this study.
For the purpose of the accuracy assessment using STEP, the different reference datasets were reprojected and merged. A stratified sampling was then applied to the dataset. As the geometries of the objects play a vital role, they were checked for inconsistencies and repaired using the toolsets in QGIS.

V. RESULTS
The resulting LC map includes the municipality of Aarhus at a resolution of 20 cm, with ten thematic classes and three additional classes derived from cadastral data, namely, intensive permanent, intensive temporary, and extensive agriculture. Fig. 4 gives an overview of this dataset. The distribution of the LC classes is given in Table III. As described in the introduction, the authors opted not to use a pixel-based confusion matrix for the accuracy assessment, but rather a polygon-based approach as it is more in line with the OBIA approach.
In comparison with the dataset compiled from existing sources mentioned in the methodology section, the reported accuracies suggested a low agreement between the reference data and the classification. However, visual comparison between these two suggests a higher agreement. This discrepancy was traced to four different reasons, namely, inconsistencies in spatial resolution, i.e., medium-resolution raster datasets being compared with high-resolution objects, temporal inconsistencies, oversegmentation, and a partial lack of reference data. Tables IV and V tabulate the two accuracy assessments, where  Table IV uses available reference data, and Table V manually classified objects as a reference. Since the LU class agriculture was not generated in this classification procedure, and only added in a post processing, it was excluded from the accuracy assessment.
In the accuracy assessment Tables IV and V, the values can range from 0 to 1, where 1 indicates the highest agreement for each metric. The average values of the metrics shape, theme, A STEP similarity matrix assigns a value of 1 to a correctly classified object and 0 to a misclassified object, indicating poor shape, thematic, edge, or positional similarity of the object. The analysis reveals that the approach with manual reference data performed better, especially for natural features or features without clear boundaries.
The LC class bare soil shows some of the lowest values throughout the metrics. However, all other LC classes show higher values, with the exception of both shrub classes. Overall, the accuracy assessment using the manual reference dataset shows higher agreement.
Tables VI and VII give more insight into the thematic accuracy of the classes. Again, Table VI shows the confusion matrix based on the available reference dataset, whereas Table VII shows the values based on the manual reference. Due to the setup of the STEP metrics system, the values represent the overlap between the reference object and the classified object as a percentage.
The errors of commission and omission are clearly visible. As seen from the comparison between Tables IV and V, Table VII shows higher thematic agreement. The errors of commission and omission are mostly limited to classes within hierarchy level, as, for example, the confusion between evergreen trees and round shrubs.
As mentioned in the introduction, similar studies have also reported lower accuracies when using polygon-based accuracy measurements, ranging from 42% to 96% overall accuracy [28]. With the exception of bare soils, this study is in line with the previous ones. Upon visual comparison, the final LC product and the underlying DOPs in this study suggest a high accuracy for which we accept the results.

VI. DISCUSSION
This study aims at a dual purpose, 1) to create a highresolution LULC dataset based on publicly accessible data, with a methodology that is potentially transferable to other areas, and 2) to assess its accuracy in a way that is more adapted to the OBIA approach than a conventional confusion matrix based on properly assigned or misclassified pixels.
The advantage lies in a semiautomated classification procedure, where users have to provide reference polygons, as available datasets did not prove to be a reliable source for reference objects. Using open access datasets makes this approach transferrable to other Danish cities and shows opportunities for other countries.
In terms of the shape-based accuracy assessment, the concept is more suitable for the OBIA approach. The limitation is set in the reference data base although as such reference data demands  III  DISTRIBUTION OF LAND-COVER CLASSES IN THE STUDY AREA, RELATIVE,  ABSOLUTE, AND TOTAL   TABLE IV  ACCURACY ASSESSMENT OF THE LAND-COVER PRODUCT, BASED ON  AVAILABLE REFERENCE a high accuracy at the same temporal and spatial extents. Most likely, either temporal and/or spatial processing status differs from the created product. Moreover, there is a high degree of dependence on such reference data, on its documented quality and on it representing the same categories as the ones defined for the new LC product. As an OBIA classification for urban LULC might not have this corresponding counterpart in reference data, it leads to a need for manually digitized reference polygons. This step requires expert knowledge, ideally in both DOP interpretation and knowledge of the local landscape. However,  this additional task translates into a more meaningful accuracy assessment. Studies with a similar setting have undertaken their accuracy assessment with a conventional confusion matrix, albeit with other classification approaches. Overall accuracies are reportedly higher, usually within 70% to 90%, with an inherently missing focus on polygon metrics [56], [57]. The derived LC classification is very promising as it yields good accuracy from an approach that requires minimal user interaction and is based on openly accessible data that are available for all of Denmark and increasingly also in other countries. A potential deficit is the misclassification of bare soils. Here, the spectral limitations of the DOPs come into play, as other studies often rely on short-wave infrared to identify bare soils [58], [59]. The spectral response of bare soils is more distinctively featured in this spectral range, and separating bare soils from artificial LC classes can be done more accurately. For this reason, the bare soil category is inconsistent and often wrongly assigned to mineral surface. Another difference between the accuracy assessments lies within the agreement of water bodies. While having distinct spectral characteristics and are, thus, easily classified correctly, their extent can change over time and differ between map producers. While the reference dataset for this study focuses on the actual riverbanks and shorelines as seen from the ground, our mapping approach inherently uses a bird's-eye perspective, which can lead to the occlusion of water by, e.g., overhanging trees. This case is in fact true for every other LC class with a low elevation, being potentially occluded by higher structures, such as trees covering houses or streets and roofs overhanging grass or mineral surface. This manifests as different values for the edge shape metric, differing between the manual and the readymade reference dataset. Other differences in the two accuracy assessments can be traced to inconsistencies in the available reference data. Especially elevated vegetation was sometimes misclassified, e.g., lines of trees as linear shrubs. With the use of local knowledge and field trips, as well as Google Maps and Google Street View, we were able to minimize the errors in our manual reference dataset.
With regards to the aspect of interoperability with municipal data, we conclude that the MMU is applicable as it reflects the size of natural and artificial LC within the area of interest. We have, thus, continued the existing mapping of public ground into privately owned parcels, providing additional insight into their LULC composition.
In addition, this study highlights the potential constraints when working with openly accessible data, especially with regard to the provision of sufficient metadata. For this study the seamlines and dates for the individual DOPs, which are critical information for data analysis, were not readily provided but had to be acquired by approaching the Danish Agency for Data Supply and Efficiency. Such lack of readily available metadata may lead to an error-loaded documentation of the resulting product.
For application by cadastral authorities, the results show a high level of categorical differentiation. For this reason, the derived product is a value-added dataset, e.g., the produced LULC map has by Aarhus municipality already been applied for assessing the potential for the rewilding of lanes within privately owned business areas. Furthermore, as part of the EU H2020 project REGREEN, the map is currently also being applied as input for the assessment of a variety of ES, such as flood mitigation, air pollution removal by trees, heat mitigation, recreation, and biodiversity [60]. However, it becomes obvious that the applied methodology is quite complex and demands sophisticated software techniques.
The latter may be a first hindrance when initially working with it. Yet, after having established the sequence of processing steps, possible updates could be carried out faster and more precisely. In addition to technical experiences, the presented process and resulting product greatly add to the detection of green space composition. More generally speaking, the hierarchical approach and the resulting classification are inherently structured in such a way that they can be easily aggregated to different super classes, depending on further aspects of investigation.

VII. CONCLUSION
The method of generating an LULC product with a submeter spatial resolution and a high thematic resolution presented in this study is a novelty, as the approach used is highly transferable and could potentially be used to classify a whole country, e.g., Denmark. With the emphasis on structural vegetation in urban areas, this LULC product provides an addition to existing municipal cadastral data and serves as a basis for urban planning. The product is potentially useful for urban planners as, in combination with information on parcel accessibility, details on the provision of urban green spaces for citizens can be explored and areas with a lack thereof can be identified.
The benefit of the approach proposed in this study is its transferability. Further work may, therefore, focus on time series analysis as the datasets are updated in a fixed cycle, and potential transfer to other countries with national DOP and DSM campaigns is feasible.