Mapping past land cover on Poitiers in 1993 at very high resolution using GEOBIA approach and open data

The land cover data presented here is a reconstruction of the past landscape (1993) at Very High Resolution (VHR) for the city of Poitiers, France. This reconstruction is based on multiple sources of images and data. We combined the strengths of both mono-temporal and multi-temporal classifications. Orthophotos were created at a spatial resolution of 0.5 m using aerial raw images from the French National Geographic Institute (IGN), taken during two aerial missions in July and August 1993. These orthophotos were merged at a spatial resolution of 5 m to conduct a first object-based classification using Landsat-5 TM images. The goal was to identify croplands, grasslands, coniferous and deciduous forests, urban areas, water bodies, and shadows. This learning-based classification employed a dataset consisting of 1371 polygons and demonstrated strong classification performances, achieving an overall accuracy of 86.31% and a kappa index of 0.832. On the other hand, mono-temporal classifications at a 0.5 m spatial resolution were carried out on each orthophoto to extract trees and herbaceous vegetation, especially in urban contexts. As mono-temporal classifications contained less information, we used a larger number of polygons for the learning step: 3849 and 5173 polygons for the northern and southern classifications, respectively. The segmentation step performed better in urban areas compared to rural areas. Consequently, the performance of classifications was evaluated separately for both contexts. Urban areas exhibited excellent performances, achieving kappa indices of 0.897 and 0.881 for the northern and southern classifications, respectively, whereas only tree vegetation was accurately detected in rural areas. To compensate for the lack of information such as buildings, railways, or roads, we modified the BD TopoⓇ dataset from IGN. This land cover map provides highly detailed information, facilitating the understanding of urban sprawl and changes in urban and rural vegetation surrounding the city of Poitiers. Due to these reasons, this freely accessible map can be utilized by researchers, land managers, and private companies for addressing urban and ecological challenges.

extract trees and herbaceous vegetation, especially in urban contexts.As mono-temporal classifications contained less information, we used a larger number of polygons for the learning step: 3849 and 5173 polygons for the northern and southern classifications, respectively.The segmentation step performed better in urban areas compared to rural areas.Consequently, the performance of classifications was evaluated separately for both contexts.Urban areas exhibited excellent performances, achieving kappa indices of 0.897 and 0.881 for the northern and southern classifications, respectively, whereas only tree vegetation was accurately detected in rural areas.To compensate for the lack of information such as buildings, railways, or roads, we modified the BD Topo R dataset from IGN.This land cover map provides highly detailed information, facilitating the understanding of urban sprawl and changes in urban and rural vegetation surrounding the city of Poitiers.Due to these reasons, this freely accessible map can be utilized by researchers, land managers, and private companies for addressing urban and ecological challenges. ©

Value of the Data
• The Geographic Object-Based Image Analysis (GEOBIA) method used here demonstrates the powerful combination of i) multi-temporal classification using satellite images to differentiate croplands from grasslands for example and ii) mono-temporal classification to recover small objects like hedgerows or isolated tree in urban and rural areas.
• This land cover map is valuable for Poitiers that have been experienced a strong expansion of its urban areas since 1993 showing now a High-Speed-Railway, more housing areas, larger industrial areas, and larger surrounding towns.Such a map allows to monitor the landscape through time, representing thus a useful data in land management.• As land cover changes is crucial to land management, this map will help to understand changes from 1993 to now for urban, agricultural issues but also their impact on ecological processes.Data will be easily used in GIS applications for any users.

Data Description
The land cover analysis of Poitiers in 1993 used aerial images at a Very High Resolution (VHR) and Landsat-5 TM images at High Resolution (HR).In addition, training and validating polygons were created to conduct learning-based classifications (i.e., Random Forest classification) et evaluate their performance ( Figs. 1-6 ).
Digital aerial images were provided by the IGN, contained a NIR, red and green channel and were acquired in summer 1993 over the area of interest.They were acquired during two different missions: 20 tiles for the north of Poitiers were taken on 16 August 1993 and 30 tiles for the south of Poitiers were taken on 28 July 1993 ( Table 1 ).It is important to note that different date acquisition can influence spectral information due to phenological variation of the vegetation.Orthorectification was conducted on each tile using ERDAS IMAGINE version 2015 (Appendix A).The images were projected into the Lambert-93 (EPSG: 2154) projection based on the RGF93 geoid (IAG GRS 1980 ellipsoid) and were resampled to a spatial resolution of 0.5 m.Then, mosaics were performed from these orthorectified tiles for each mission (Northern and Southern   images).The topographic correction, radiometric calibration, atmospheric and anisotropy correction were not performed because the needed parameters (e.g., atmospheric composition, aerosol types) were not available.Even though the pre-processing step is not optimal as in many other conditions or studies, the main objective was to get reliable information to conduct a classification on these images.
We used multi-temporal Landsat scenes to catch variations of vegetation.Landsat images have been widely used to consider the multi-seasonal variation of spectral information of vegetation due to the different phenology of plant species [18 , 23] .We used seven Landsat Thematic Mapper 5 scenes acquired between January 1992 and August 1993 ( Table 1 ) to differentiate, for example, cropland from grassland or deciduous from coniferous.The satellite images were   The OTB open-source project [9] was used to conduct the remote sensing analyses and QGIS project [17] was used for the remote sensing and GIS analyses.

Steps of multi-temporal classification
Two classifications were conducted successively using multi-temporal data to get the main land cover categories and using mono-temporal data to extract small vegetation elements such as isolated trees or hedgerows.To identify the main objects such as impervious areas, forests, grasslands, or croplands, the two orthophotos were merged and resampled to a spatial resolution of 5 m using bilinear method to facilitate the segmentation.The following algorithms were used from the OTB project [9] .

Segmentation
The segmentation is an important step of the object-based classification as its results influence the classification accuracy ( [8] ; Jian [11 , 21] ).The objective is to create groups of spatially close pixels with similar radiometric characteristics [1 , 4] .We used the MeanShiftSegmentation algorithm, this procedure contains four steps: i) a smoothing step facilitating the segmentation, ii) the segmentation step creating segments, iii) the merging step that fusion small segments to similar neighbour objects, iv) the vectorization step to convert pixel segments into vector segments and calculate mean and standard deviation of each raster band within segments [9] .The spatial radius and the range radius manage the smoothing effect and the capacity to preserve edges between objects.These parameters should not be too high to avoid under-segmentation and not too low to avoid over-segmentation.The minimum segment size avoids very small objects.In addition, the spatial resolution of the input image impacts the segmentation step.For the LargeScaleMeanShift algorithm, we selected a spatial radius of 5 pixels, a range radius of 15 pixels and a minimum segment size of 12 pixels.

Feature extraction
The feature extraction allows to distinguish categories during the classification step.Thus, segments must contain information related to spectral, texture, 3D or geometric characteristics [15 , 20] .For each segment, we computed the mean and standard deviation value of the NIR, Red and Green channel of the orthophoto, the mean value of raster bands of Landsat images (i.e., Blue, Green, Red, NIR, MIR channel) (Appendix B).To improve the detection of water, the modified normalized difference water index (MNDWI) [22] was calculated for the 7 Landsat scene and their mean value were computed within segments.Moreover, the number of pixels of each segment was extracted for the classification step.

Random forest classification
Random forest has been introduced by Breiman [2] as an ensemble of classification methods.This approach uses classification and regression trees (CART; [3] ) where each tree gives a classification to identify the most prevalent class.Then, an out-of-bag (OOB) error is calculated to estimate the rate of misclassified samples.This learning approach is widely used in remote sensing classification for its capacity to manage many variables, its stability and robustness.Moreover, there is few parameters to tune which are i) the number of trees within the forest and ii) the depth of each tree.The number of trees in the forest was set to 250 and the depth of tree to 7, corresponding to the square root of the number of variables (i.e., 49; Appendix B) [18] .

Training and validating polygons
The selection of training and validating polygons is crucial for the classification accuracy.The best classification accuracy is obtained when statistics of categories show low intra-class variation and high inter-class variation.Thus, they must represent statistical variation within each class.Training and validating polygons were chose across all the studied area, to avoid spatial autocorrelation, which can have an influence on the evaluation of classification accuracy ( Fig. 1 ; [12] ).From the segmentation, we selected segments by photointerpretation (based on the orthophotos and Landsat images) to construct the training dataset.Seven target classes were identified for the classification: Impervious surface contained buildings, roads, bare soils; cropland contained annual crops; grassland contained meadows and herbaceous vegetation; coniferous forest; deciduous forest; water contained the main water bodies and courses; and shadow ( Table 2 ).The number of polygons was roughly based on the relative importance of each category considering the classification objectives.
The evaluation of classification accuracy is based on the number of pixels well assigned.Thus, we designed validating polygons with similar area corresponding to approximately 12 pixels ( Fig. 2 ; [18] ).To evaluate the global performance of the classification, the overall accuracy and the kappa index were computed.Three metrics were also calculated at the class-level: the precision (also called Producer's accuracy), the recall (also called User's accuracy) and the F-Score (which is the harmonic mean of the precision and the recall).

Steps of mono-temporal classifications
We conducted classifications at VHR from orthophotos to a 0.5 m of spatial resolution using a similar workflow as previously.However, each of the original orthophoto were treated separately because of their spectral differences induced by their acquisition condition, in order to increase the classification performances.Conversely, the main information for the multi-temporal classification was contained in the Landsat images, so the differences between the two orthophotos were negligible.Consequently, two classifications were realized (i.e., Northern orthophoto & Southern orthophoto).Due to the nature of the input data, this procedure was highly similar to the method proposed by Morin et al. [14] .This approach was mainly developed for highly heterogeneous and complex landscapes such as the urban environment.

Segmentation
This segmentation step aims to create appropriate segments for small objects such as the canopy of a tree.For the LargeScaleMeanShift algorithm, we selected a spatial radius of 5 pixels, a range radius of 28 pixels and a minimum segment size of 60 pixels.A large value of the range radius induces a strong smooth effect facilitating the segmentation that is particularly interesting for very high resolution images.Moreover, the minimum size of 60 pixels is a good compromise for detecting isolated trees and avoiding oversegmentation.

Feature extraction
There are some limits to conduct a classification only from the spectral information of a mono-temporal image.For example, the absence of the height of each object can impedes the identification of vegetation strata.However, textures have been used with success to distinguish vegetation strata or urban objects [5 , 7 , 13] .Textures have been described as the smoothness, the regularity, or the coarseness of an image (Gonzalez & Woods 2002).This aspect of an image can be statistically captured from a grey level co-occurrence matrix [10 , 16] that represent the spatial homogeneity around the pixels using a specific size of moving window.The matrix quantifying the variation in grey levels is then used to compute the textures.Thus, spectral and texture information were calculated to improve the classification accuracy.The Normalized Difference Vegetation Index (NDVI) [19] was computed, a widely used index to detect vegetation in remote sensing applications.The brightness was calculated as the mean of NIR, Red and Green band.In addition, four Haralick's textures were computed from the brightness using the Haral-ickTextureExtraction algorithm with a window size of 7 × 7 pixels: Energy, Entropy, Correlation, Contrast .The average of the NIR, Red, Green, the NDVI and the four texture was included within segments as well as the standard deviation of the Green and the four textures (Appendix B).Segments were filtered before the classification to remove the shadow segments which showed a brightness value lower than 95 for the Northern orthophoto and 80 for the Southern orthophoto.

Random forest classification
We used the Random Forest algorithm for the two classifications.For the parameter values of the TrainVectorClassifier algorithm, we chose a value of 250 for the number of trees and a value of 4 for the depth of tree (i.e., the square root of the number of the input variables) [18] .A stratified subsampling followed by a cross-validation process was also used by dividing the training dataset into five roughly equal subsets to conduct 5 different classifications using the VectorClassifier algorithm [6 , 14] .Then, these classifications were merged with a majority voting (i.e., FusionOfClassifications algorithm).The overall accuracy and the Kappa index were computed to estimate the global performance of the two classifications.The precision, recall and F-Score were calculated for each land cover category.

Training and validating polygons
The herbaceous vegetation, the wooded vegetation and the artificial area were the three target categories.Because less information was available compared to the first classification that used a multi-temporal approach, a cross validation process with a larger number of training segments were chose ( Table 3 ; Figs. 3 , 5 ).To estimate the accuracy of each classification, we randomly designed validating polygons with similar area corresponding to approximately 60 pixels ( Figs. 4 , 6 ; [18] ).As this method was optimized for urban context but was conducted on the whole area of interest in this paper, we evaluated the accuracy of the two classifications in urban context and rural context separately and together.

Segmentation results for the two different approaches
The segmentation is an important process to create semantic objects.The two segmentations were conducted from the same mono-temporal images but at two different spatial resolutions (i.e., 0.5 m and 5 m).Results showed that the resolution and the chosen parameters influenced the created segments.Differences rely on the complexity and heterogeneity of the environment.In rural areas where heterogeneity and complexity are low, segmentation at HR was more relevant to delimitate crops ( Fig. 7 ).Conversely in urban areas which are very complex and hetero- geneous landscapes, VHR segmentation well performed to identify artificial from tree and grass objects ( Fig. 7 ).The segmentation at HR produced 359,750 segments for the whole area of interest and the segmentations at VHR produced 710,688 segments for the Northern images and 1514,964 for the Southern images.

Analysis of accuracy assessments
The multi-temporal classification showed good classification results with an overall accuracy and a kappa index of 86,31% and 0.832 respectively.Regarding the accuracy at the class-level, impervious surface, cropland and shadow were well identified with high F-Score (i.e., around 90%) ( Table 4 ).Grassland showed relatively high F-Score but had confusions with impervious surface (in urban context) and deciduous forest (mainly due to meadows close to forest and for which Landsat pixel overlapped) ( Table 4 ).Deciduous forest was well detected but had some true deciduous pixels classified in other categories (recall value of 77%) such as coniferous forest or grassland.Coniferous forest was overestimated by having a low precision value and showing a confusion with deciduous forest and water and well classified true coniferous pixels (i.e., recall value of 100%).The water category was underestimated with a good precision value but a low recall value identifying a confusion mainly with the impervious surface category ( Table 4 ).
It was more contrasted for the VHR classifications.The accuracy of the northern classification was quite low, considering only three land cover categories, with an overall accuracy of 80.74% and a kappa index of 0.713 for the northern classification but quite better for the southern classification with an overall of 88.36% and a kappa index of 0.814.However, the use of independent validating polygons in rural and urban areas allows to estimate the accuracy of both classifications in these different contexts.
In urban context, the classification accuracy was high with an overall accuracy of 93.20% and kappa index of 0.897 for the northern classification and an overall accuracy of 92.20% and a kappa index of 0.881 for the southern classification.Conversely, in rural context results were lower with an overall accuracy of 6 8.4 9% and a kappa index of 0.527 for the northern classification and an overall accuracy of 85.25% and a kappa index of 0.746 for the southern classification.Wooded vegetation was well detected in both context by showing high F-Score values ( Table 5 , Appendix C).

Compiling data
We compiled the created data and available institutional databases to build a map at very high spatial resolution ( Table 6 ; Fig. 8 ).Buildings, roads, railways and water courses and bodies were manually modified from existing data (BD Topo of IGN) by photointerpretation.Deciduous and coniferous forests, croplands, grasslands and artificial areas were extracted from the multi-temporal classification.Water from the HR classification was not used due to its coarse resolution and because this category was mainly integrated into the classification to not overestimate the impervious surfaces which could have close radiometric information.Urban vegetation strata, detected within the artificial areas of the multi-temporal classification, were extracted from the mono-temporal classifications.In addition, the wooded vegetation of mono-temporal classifications was also integrated in rural areas du to its good detection.Shrublands were added by creating polygons using photointerpretation from orthophotos.

Limitations
Despite the high accuracy of the land use map, shrub areas are difficult to classify and can be classified as either tree or herbaceous vegetation.

Appendix A
We used the method of geometric modelling based on the collinearity equation and implemented in the Leica Photogrammetry Suite (integrated in ERDAS Imagine).The ERDAS Imagine software is a simple and useful tool to visualize and manipulate geographic imagery data.The ERDAS Imagine software can perform a wide variety of tasks such as orthorectification (Long & Srihar, 2004) for which it is considered the most reliable.Moreover, during orthorectification, the software incorporates the following parameters: (i) The reference system in which the image will be georeferenced, (ii) The Digital Terrain Model (DTM) used to derive the elevation of the Ground Points Control (GCP), (iii) The camera and raw image parameters (see Fig. A1 ): -the focal length in millimeters, -the coordinates in millimeters in the system of the photo of the main point of symmetry, which symmetry point, which are ideally [0,0] but not necessarily, -the number of fiducial points used for the interior orientation, -the coordinates of the fiducials in millimeters in the coordinate system of the photo.According to Rocchini et al., 2012, The ortho-rectification steps ( Fig. A2 ) include: • Interior orientation: The interior orientation of the photograph is needed to relate the image coordinate system (in pixels) to the photo coordinate system (camera sensor).This is done by graphically marking the fiducial marks into the scanned image.• Exterior orientation: Exterior orientation of the aerial photograph is based on the selection of n GCPs, taking into account their x and y (coordinates of each GCP of the incorrect input image in the image coordinate system), their x0 and y0 (corresponding GCP coordinates of the rectified image in the target coordinate system defined by projection parameters), and z coordinates.
• Rectification and resampling: The parameters computed during the interior and exterior orientation steps are used to rectify the image, ultimately passing from the input x, y image coordinates to the x0, y0 standard map coordinates in the target projection.Finally, after having retrieved the parameters from the IGN and the DTM, we have orthorectified the 50 tiles taken in 1993 in the surroundings of Poitiers.We have, compiled these tiles into two mosaics (corresponding to the two periods in which the photos were taken: July in the South and August in the North).We trimmed these two mosaics in order to obtain a global study area of 15 km by 15

Fig. 1 .
Fig. 1.Distribution of the training polygons for the multi-temporal classification.

Fig. 2 .
Fig. 2. Distribution of the validating polygons for the multi-temporal classification.

Fig. 3 .
Fig. 3. Distribution of the training polygons for the mono-temporal classification on northern images.

Fig. 4 .
Fig. 4. Distribution of the validating polygons for the mono-temporal classification on northern images.

Fig. 5 .Fig. 6 .
Fig. 5. Distribution of the training polygons for the mono-temporal classification on southern images.

Fig. 7 .
Fig. 7. Results of segmentation and classification at both scales.
2023 The Author(s).Published by Elsevier Inc.This is an open access article under the CC BY license ( http://creativecommons.org/licenses/by/4.0/ ) Data source location Poitiers is an important historical metropolis of the west of France (46.34 °N, 0.20 °E) with 131,499 inhabitants (census in 2017).The area of the land cover map is about 225km ².

Table 1
Characteristics of used images.

Table 2
Training dataset used for the multitemporal classification.

Table 3
Training and validating polygons used for each mono-temporal classification.

Table 5
Performances of both classifications evaluated in urban context, rural context and urban-rural context.