Spatiotemporal Changes in 3D Building Density with LiDAR and GEOBIA: A City-Level Analysis

Understanding how, where, and when a city is expanding can inform better ways to make our cities more resilient, sustainable, and equitable. This paper explores urban volumetry using the Building 3D Density Index (B3DI) in 2001, 2010, 2019, and quantifies changes in the volume of buildings and urban expansion in Luxembourg City over the last two decades. For this purpose, we use airborne laser scanning (ALS) point cloud (2019) and geographic object-based image analysis (GEOBIA) of aerial orthophotos (2001, 2010) to extract 3D models, footprints of buildings and calculate the volume of individual buildings and B3DI in the frame of a 100 × 100 m grid, at the level of parcels, districts, and city scale. Findings indicate that the B3DI has notably increased in the past 20 years from 0.77 m3/m2 (2001) to 0.9 m3/m2 (2010) to 1.09 m3/m2 (2019). Further, the increase in the volume of buildings between 2001–2019 was +16 million m3. The general trend of changes in the cubic capacity of buildings per resident shows a decrease from 522 m3/resident in 2001, to 460 m3/resident in 2019, which, with the simultaneous appearance of new buildings and fast population growth, represents the dynamic development of the city.


Introduction
As the urban is extending, monitoring change over space and time is important for supporting decisions about appropriate development practices and land resource use. Increasing demands in urban management sectors need the coordinated use of remote sensing and a geographic information system (GIS) for monitoring urban intensification [1]. Building density is one of the most important indices for city monitoring [2]. Urban structure analyzes are carried out in two dimensions (2D), where building footprints are solely considered, or in three dimensions (3D), where building heights and footprints are considered [3]. However, in growing countries, obtaining accurate, complete, and current building information from cadastral data is hard to come by [4,5]. Two metrics often used in the regulations of city planning are the floor area ratio (FAR) and the building coverage ratio (BCR). FAR is the relationship of the gross floor area of a building to the total buildable area of the lot/parcel. BCR is calculated by dividing the total buildable area of a lot/parcel by the total area of the lot/parcel [6].
Remotely sensed imagery and GIS analysis are often used to extract building footprints. Various classification methods based on high-resolution satellite imagery have been used to derive in this case. The question of buildings density is closely connected to urbanization and how cities may grow in the future. Thus, it is important to track the city over time and across space: how does the three-dimensional density of buildings evolve in Luxembourg City over time?
In this paper, we propose a 3D index to detect changes in building volume density in the last 20 years in Luxembourg City employing ALS LiDAR point cloud (2019) and GEOBIA approach of archival orthophoto maps (2001,2010). Our analysis was supported by datasets collected from orthophotos during the last 20 years concerning demolish and newly constructed buildings. We process ALS classified point clouds to extract the points reflecting a rooftop surface and generate 3D models only for buildings class. A building's volume density has been calculated for all buildings in the city including residential and non-residential categories. The detection of archival footprints of buildings was carried out using the object-based method. The complex integration of final datasets is unique in research and gives us very detailed information about the urban spatiotemporal (4D) changes, considering volume, vertical, and horizontal buildings' changes.
Our methodology is general and can be applied to several cities when LiDAR data is available.
In this study, our key objectives are: (1) To automate 3D building modelling and object-based footprints of buildings extraction; (2) To define Building 3D Density Index (B3DI); (3) To quantify changes in volume density in Luxembourg City over the last two decades (between 2001, 2010, and 2019).

The Geographical Localisation of the Study Area
Luxembourg City is the capital of the Grand Duchy of Luxembourg and the country's most populous city (Figure 1). Luxembourg City at municipality level occupies an area of 51 km 2 . The city has a complex topographic heterogeneous configuration, lying at the different levels with parks and forested areas. The city hosts several EU institutions, and it is one of the major international financial hubs. Changes in land use in Luxembourg in recent decades reveal a rapid urban development and land artificialization [41,42]. According to Corine Land Cover (CLC), made available by the European Environment Agency, urban land use in Luxembourg City covered about 35% and about 39% in 1990 and 2018, respectively. At the same time, the industrial, commercial, and transport uses increased by 3.3% whereas agricultural areas decreased by 6.6%, and forest and semi-natural areas decreased by 2.3%.

Datasets and Data Sources
In this study, we used ALS LiDAR point clouds from airborne mission obtained in February 2019 for Luxembourg City (mean point density was 15 points/m 2 with horizontal and vertical accuracies was to ±3 cm and ±6 cm, respectively; in the area of the Luxembourg City, the average density was 25 points/m 2 ) and archival RGB and CIR (color-infrared) aerial orthophoto maps from 2001, 2010, and 2019 (ground sampling distance (GSD) 50/25/20 cm), provided the Luxembourg Government's open data (data.public.lu). The acquired point clouds were previously classified (ASPRS standard). The dataset was supplemented with cadastral information from 2019. However, the footprints of buildings were not updated and contained some errors in the geometries, thus, we generated footprints based on ALS LiDAR point cloud and added information about the use of buildings (residential or non-residential), according to cadastral data. We also used an archival database from the Housing Observatory of an on-going project [43,44] concerning the demolition/reconstruction (2001-2019) of buildings for Luxembourg, to detect buildings where there was a significant height change during the last two decades. This data was collected through the visual comparison of orthophotos every 3 years starting in 2001. To estimate the volume of buildings existing only in 2010, we acquired a layer of "Copernicus Urban Atlas-Building Height" from the European Environment Agency, Remote Sens. 2020, 12, 3668 4 of 23 within the framework of the Copernicus program. The product is a 10 m resolution and provides the height of buildings information that is generated for urban cores of capitals as part of the Urban Atlas project. Height information is based on IRS-P5 stereo images acquired as close as possible to the defined reference year and derived datasets like the digital surface model (DSM), the digital terrain model (DTM), and the normalized digital surface model (nDSM). Table 1 lists the data used in the current study.

Datasets and Data Sources
In this study, we used ALS LiDAR point clouds from airborne mission obtained in February 2019 for Luxembourg City (mean point density was 15 points/m 2 with horizontal and vertical accuracies was to ±3 cm and ±6 cm, respectively; in the area of the Luxembourg City, the average density was 25 points/m 2 ) and archival RGB and CIR (color-infrared) aerial orthophoto maps from 2001, 2010, and 2019 (ground sampling distance (GSD) 50/25/20 cm), provided the Luxembourg Government's open data (data.public.lu). The acquired point clouds were previously classified (ASPRS standard). The dataset was supplemented with cadastral information from 2019. However, the footprints of buildings were not updated and contained some errors in the geometries, thus, we generated footprints based on ALS LiDAR point cloud and added information about the use of buildings (residential or non-residential), according to cadastral data. We also used an archival database from the Housing Observatory of an on-going project [43,44] concerning the demolition/reconstruction (2001-2019) of buildings for Luxembourg, to detect buildings where there was a significant height change during the last two decades. This data was collected through the visual comparison of orthophotos every 3 years starting in 2001. To estimate the volume of buildings existing only in 2010, we acquired a layer of "Copernicus Urban Atlas-Building Height" from the European Environment Agency, within the framework of the Copernicus program. The product is a 10 m resolution and provides the height of buildings information that is generated for urban cores of capitals as part of the Urban Atlas project. Height information is based on IRS-P5 stereo images acquired as close as possible to the defined reference year and derived datasets like the digital surface model (DSM), the digital terrain model (DTM), and the normalized digital surface model (nDSM). Table 1 lists the data used in the current study.

Methods
Monitoring of spatiotemporal changes in the volume of buildings in Luxembourg City in 2001, 2010 and 2019 required integration of multi-source datasets and different methods of data analysis. This section presents a detailed description of all steps of analysis ( Figure 2).

Methods
Monitoring of spatiotemporal changes in the volume of buildings in Luxembourg City in 2001, 2010 and 2019 required integration of multi-source datasets and different methods of data analysis. This section presents a detailed description of all steps of analysis ( Figure 2).

Generation of 3D Models of Buildings
ALS LiDAR point clouds were used to obtain precise information about the 3D configuration of the terrain and buildings in 2019. The processing was started by generating the DTM and the DSM of buildings from the classified point cloud (ASPR standard): class ground and building with 0.5 m resolution using Area Processor in FUSION software, ver. 3.70 (USDA Forest Service, Pacific Northwest Research Station). The nDSM, which represents the relative heights of buildings, was generated based on the difference of DSM and DTM rasters. The nDSM, generated from all classes of the point cloud, contains objects that could disturb the height or volume calculations of buildings (e.g., tree branches next to buildings or cranes on newly constructed buildings-example in Figure 3). For this reason, filtered nDSM, created only from buildings class, was used in further analysis. resolution using Area Processor in FUSION software, ver. 3.70 (USDA Forest Service, Pacific Northwest Research Station). The nDSM, which represents the relative heights of buildings, was generated based on the difference of DSM and DTM rasters. The nDSM, generated from all classes of the point cloud, contains objects that could disturb the height or volume calculations of buildings (e.g., tree branches next to buildings or cranes on newly constructed buildings-example in Figure  3). For this reason, filtered nDSM, created only from buildings class, was used in further analysis.  Buildings vectorization was performed fully-automatically to generate non-planar roof shapes of buildings and 3D models of buildings in LoD2 (Level of Details; CityGML standard) based on classified ALS point clouds with two classes: ground and building ( Figure 4) in TerraScan (Terrasolid) and MicroStation V8i (Bentley) software. These tools check and correct buildings' topology. The quality of the automatic buildings' vectorization depends not only on the quality of the LiDAR data processing that is done in preparation of the vectorization but also on the point density of the data. Based on 3D models, building footprints were generated and exported in the shapefile format. The minimum mapping area units (MMU) were 20 m 2 , and smaller objects were not included. The nDSM of buildings (GSD 0.5 m) was shrunk to a roof polygon generated from ALS point cloud to eliminate errors on the height model. The model of buildings was used to determine the total volume of buildings [m 3 ] and including the classification into residential and non-residential uses in Luxembourg City. The information about the category of buildings from cadastral data has been added to the shapefile of buildings in 2019. The volume of buildings was calculated for single buildings, at the districts-level and for standardization in cells (100 × 100 m) using ArcMap software (Esri) with the Zonal Statistics and Surface Volume tool, which calculates the volume of the region between a surface and a reference plane. Buildings vectorization was performed fully-automatically to generate non-planar roof shapes of buildings and 3D models of buildings in LoD2 (Level of Details; CityGML standard) based on classified ALS point clouds with two classes: ground and building ( Figure 4) in TerraScan (Terrasolid) and MicroStation V8i (Bentley) software. These tools check and correct buildings' topology. The quality of the automatic buildings' vectorization depends not only on the quality of the LiDAR data processing that is done in preparation of the vectorization but also on the point density of the data. Based on 3D models, building footprints were generated and exported in the shapefile format. The minimum mapping area units (MMU) were 20 m 2 , and smaller objects were not included. resolution using Area Processor in FUSION software, ver. 3.70 (USDA Forest Service, Pacific Northwest Research Station). The nDSM, which represents the relative heights of buildings, was generated based on the difference of DSM and DTM rasters. The nDSM, generated from all classes of the point cloud, contains objects that could disturb the height or volume calculations of buildings (e.g., tree branches next to buildings or cranes on newly constructed buildings-example in Figure  3). For this reason, filtered nDSM, created only from buildings class, was used in further analysis.  Buildings vectorization was performed fully-automatically to generate non-planar roof shapes of buildings and 3D models of buildings in LoD2 (Level of Details; CityGML standard) based on classified ALS point clouds with two classes: ground and building ( Figure 4) in TerraScan (Terrasolid) and MicroStation V8i (Bentley) software. These tools check and correct buildings' topology. The quality of the automatic buildings' vectorization depends not only on the quality of the LiDAR data processing that is done in preparation of the vectorization but also on the point density of the data. Based on 3D models, building footprints were generated and exported in the shapefile format. The minimum mapping area units (MMU) were 20 m 2 , and smaller objects were not included. The nDSM of buildings (GSD 0.5 m) was shrunk to a roof polygon generated from ALS point cloud to eliminate errors on the height model. The model of buildings was used to determine the total volume of buildings [m 3 ] and including the classification into residential and non-residential uses in Luxembourg City. The information about the category of buildings from cadastral data has been added to the shapefile of buildings in 2019. The volume of buildings was calculated for single buildings, at the districts-level and for standardization in cells (100 × 100 m) using ArcMap software (Esri) with the Zonal Statistics and Surface Volume tool, which calculates the volume of the region between a surface and a reference plane. The nDSM of buildings (GSD 0.5 m) was shrunk to a roof polygon generated from ALS point cloud to eliminate errors on the height model. The model of buildings was used to determine the total volume of buildings [m 3 ] and including the classification into residential and non-residential uses in Luxembourg City. The information about the category of buildings from cadastral data has been added to the shapefile of buildings in 2019. The volume of buildings was calculated for single buildings, at the districts-level and for standardization in cells (100 × 100 m) using ArcMap software (Esri) with the Zonal Statistics and Surface Volume tool, which calculates the volume of the region between a surface and a reference plane.

Extraction of Building Footprints Using GEOBIA and Volume Calculation
The object-based approach was applied to extract all buildings in Luxembourg City with RGB (2001, 2010) and CIR (2010) orthophoto maps using ruleset prepared in eCognition Developer 9.3 (Trimble Geospatial) software. The ruleset creation was started with generating the derivative layers Remote Sens. 2020, 12, 3668 7 of 23 based on spectral channels of CIR orthophotos. Using bands NIR (near-infrared) and Red, a normalized difference vegetation index (NDVI) layer was generated for separate areas with different spectral characteristics. In order to improve the contrast between the objects and the background and reduce the noise, which occurs in very high-resolution images, median filter images of NIR, Red, and Green bands were created as an additional input layer for the segmentation algorithm. This produces image segments that are smoother and better representing the object. principal component analysis (PCA) was performed based on orthophotos bands to generate three derivatives components for recognizing objects in a higher-level concept of computer vision. We used PCA to reduce the dimensionality of the data and to increase the quality of processing, transforming linearly correlated variables into uncorrelated. On the analyzed orthophotos in the borders of the city of Luxembourg, the greatest correlation was found between bands 2 (Red) and 3 (Green). Hence, there was a redundancy in information between these bands, which means that the reflection coefficients are somewhat correlated between the bands. For this reason, we used the second PCA component to improve buildings' classification process, because in this layer, data compression was achieved, reducing data correlation. After initial image processing, we applied a multi-resolution segmentation algorithm [45] with scale parameter: 10, shape: 0.5, and compactness value: 0.7 in eCognition (Trimble Geospatial) software. The input layers (NIR, Red, Green, PCA_2) were smoothed with a median filter and have been assigned the following image layer weights respectively: 2, 2, 1, 2 according to experience from previous work on image classification [46,47]. Using an iterative algorithm, the segments were generated, whereby objects (starting with individual pixels) are grouped until a threshold generated representing the upper object variance is reached. Thereafter, we executed the spectral difference segmentation, where neighboring objects with a spectral mean threshold: 10 (maximum spectral difference) were merged to produce the final objects and then to improve classification accuracy.
The process of classification was modular, in the first phase all segments were divided into two temporary classes (vegetated and non-vegetated) using threshold-based approach of value in Red band and NDVI layer. Next, parameters of brightness and mean values of NIR, Red, and Green bands were used to capture buildings from non-vegetated areas. In addition to spectral values, geometry, texture features, and context information was also used. Finding the optimal value of features such as asymmetry, boundary index, object size, and relative boundary to class object features were helpful for capturing buildings. In the final phase of the GEOBIA process, we removed single objects from buildings class with an area smaller than 20 m 2 that, during visual interpretation, did not correspond to a building, but were rather parts of walls, terraces, or small garden gazebos. The shape of the final objects has been smoothed using "Morphology" and "Pixel-based object resizing" algorithms ( Figure 5).
Remote Sens. 2020, 12, x FOR PEER REVIEW 8 of 24    The number of reference points required to generate the error matrix and estimate the accuracy assessment was calculated using the binomial distribution [48]. The 298 test points were generated using a random sampling pattern approach. For each class, buildings, and unclassified, the 149 segments were placed and matched with the final classification. The quality control of the GEOBIA classification was performed by visual analysis using reference data from high-resolution orthophotos. The overall accuracy (OA), producer (PA), user accuracies (UA), and Kappa coefficient were generated by using an error confusion matrix (i.e., cross-tabulation matrix between observed and predicted values) [49].
To calculate the volume of buildings in 2010, we used demolition-reconstruction shapefile [44] to detect plots of land where buildings that have changed between 2010 and 2019. For buildings that have not changed at all, we used nDSM (2019) to calculate the volume. Buildings that changed (for The number of reference points required to generate the error matrix and estimate the accuracy assessment was calculated using the binomial distribution [48]. The 298 test points were generated using a random sampling pattern approach. For each class, buildings, and unclassified, the 149 segments were placed and matched with the final classification. The quality control of the GEOBIA classification was performed by visual analysis using reference data from high-resolution orthophotos. The overall accuracy (OA), producer (PA), user accuracies (UA), and Kappa coefficient were generated by using an error confusion matrix (i.e., cross-tabulation matrix between observed and predicted values) [49].
To calculate the volume of buildings in 2010, we used demolition-reconstruction shapefile [44] to detect plots of land where buildings that have changed between 2010 and 2019. For buildings that have not changed at all, we used nDSM (2019) to calculate the volume. Buildings that changed (for example, buildings that have had a significant roof change or a different shape) during that time or that only existed in 2010, the volume was calculated from the Urban Atlas layer (2011).
To estimate the volume of buildings in 2001, several GIS operations were performed on available data. From the shapefile obtained in the GEOBIA analysis in 2001, all buildings that changed after 2001 were separated using the intersection with Housing Observatory demolition-reconstruction data. For unchanged buildings between 2001-2019, the height was calculated using LiDAR 2019 data (layer 1). The buildings that have changed between 2001 and 2019 were divided into two layers: the height of the buildings changed before 2010 (layer 2) and the height of the buildings changed after 2010 (layer 3). For layer 2, as it was not available in our datasets, the volume of buildings was calculated using the estimated height values. For this purpose, an additional layer was made with the difference in height of modified buildings between 2010 and 2019 and the mean value of change in building height during this period was calculated (based on 100 randomly selected buildings on nDSM in 2019 and Urban Atlas-Building Height in 2011), which was 4.6 m. This value was used to estimate the height value for 2001 buildings that have changed by subtracting the value 4.6 m from the Urban Atlas-Building Height raster. For the second group of buildings, the height value was calculated directly from the Urban Atlas data. Layer 4 was created for these buildings that existed only in 2001 and were demolished after without any reconstruction structures. To do so, it was necessary to remove all buildings that existed or have changed after 2001 with a 3-step process to keep only buildings in 2001 (layer 4). Based on vector buildings layer obtained in GEOBIA process (2001), we used historical data from demolition-reconstruction layer and clip tool (Arc Map function) to separate all buildings unmodified (step 1), modified buildings before 2010 (step 2), and modified ones after 2010 (step 3). The result was layer 4 containing the unique buildings demolished after 2001. The mean of buildings heights was computed using the Urban Atlas layer.

Data Fusion to Calculate Building 3D Density Index
To quantify the built-up area, we proposed a spatial index for assessment of urban form in three dimensions. The 3D index describes the density of buildings in 2001, 2010, 2019. To calculate the quantitative spatial volume density of urban buildings, to make it closer to realistic status [50], a volumetric descriptor of Building 3D Density Index (B3DI) was proposed, expressed as (Equation (1)): where V B demonstrates the total volume of buildings [m 3 ], V i volume of a single building in the area of interest (AOI), n number of buildings, and S represents AOI land area [m 2 ]. The B3DI was calculated at entire city and district scale (residential and non-residential buildings separately in 21 districts according to cadaster data (Luxembourg City cadastral districts: Basse Pétrusse, Beggen, Bonnevoie, Cessange, Clausen, Dommeldange, Eich, Gasperich, Grund, Hamm, Hollerich, Kockelscheuer, Limpertsberg, Merl-Nord, Merl-Sud, Neudorf, Pfaffenthal, Pulvermuehl, Rollingergrund, Ville Haute, Weimerskirch.), per plot of land in the city and in a grid with cell dimension of 100 × 100 m, which has been adapted to Luxembourg municipal buildings. A higher B3DI means higher 3D density of the building, and therefore a higher intensity of regional land use.
Due to the lack of archival cadastral data of building footprints in 2001 and 2010 and incomplete data from 2019, an object algorithm and automatic 3D vectorization of LiDAR point cloud were developed to extract buildings in the city of Luxembourg at different years.  (Figures 7 and 8).     Based on ALS point cloud acquired in 2019 with a high density about 25 pts/m 2 , we generated models (LoD2) which include detailed modeling of the shape of roofs (without architectural elements such as skylights, chimneys, antennas, solar panels, etc.) with vertical projection to the ground ( Figure 9). The high density of ALS point cloud >10 pts/m 2 allows to generate accurate models with details and roof constructions [51]. At a point cloud density of 25 pts/m 2 , the XY accuracy of constructed 3D building models was approximately 20 cm. The accuracy of the roof planar surface for the Z coordinate was approximately 10 cm and for the edges 20 cm. The roof shape of the building has been exported to shapefile format and divided into residential and non-residential buildings ( Figure 10). A visual assessment of the cadastral shapefile was performed based on orthophoto (2019) to check the correctness of the administrative data.  Based on ALS point cloud acquired in 2019 with a high density about 25 pts/m 2 , we generated models (LoD2) which include detailed modeling of the shape of roofs (without architectural elements such as skylights, chimneys, antennas, solar panels, etc.) with vertical projection to the ground ( Figure 9). The high density of ALS point cloud > 10 pts/m 2 allows to generate accurate models with details and roof constructions [51]. At a point cloud density of 25 pts/m 2 , the XY accuracy of constructed 3D building models was approximately 20 cm. The accuracy of the roof planar surface for the Z coordinate was approximately 10 cm and for the edges 20 cm. The roof shape of the building has been exported to shapefile format and divided into residential and non-residential buildings ( Figure 10). A visual assessment of the cadastral shapefile was performed based on orthophoto (2019) to check the correctness of the administrative data. Based on ALS point cloud acquired in 2019 with a high density about 25 pts/m 2 , we generated models (LoD2) which include detailed modeling of the shape of roofs (without architectural elements such as skylights, chimneys, antennas, solar panels, etc.) with vertical projection to the ground ( Figure 9). The high density of ALS point cloud >10 pts/m 2 allows to generate accurate models with details and roof constructions [51]. At a point cloud density of 25 pts/m 2 , the XY accuracy of constructed 3D building models was approximately 20 cm. The accuracy of the roof planar surface for the Z coordinate was approximately 10 cm and for the edges 20 cm. The roof shape of the building has been exported to shapefile format and divided into residential and non-residential buildings ( Figure 10). A visual assessment of the cadastral shapefile was performed based on orthophoto (2019) to check the correctness of the administrative data.

Urban Volumetry and Building 3D Density Index (B3DI)
The volume calculation and 3D spatial analyses performed on the model of buildings from 2019 show that total building volume in Luxembourg City is about 56 million m 3 , which is 460 m 3 per resident. The mean volume of all types of buildings is 2349 m 3 per building. The tallest building (in the business center of Kirchberg) reaches 114 m in height (84,683 m 3 ), and the largest volume of a single building is 580 thousand m 3 in Gasperich district (Cloche d'Or Shopping Center). The other largest buildings in the city are hypermarkets in Kirchberg, European Union edifices, some banks, ArcelorMittal company, European School, and hospitals.
The volume in 2019 was calculated separately for the residential and non-residential buildings. The volume of all residential buildings in Luxembourg City is 30 million m 3 . The mean volume of a residential building is 1460 m 3 and the mean height is 12.5 m. Buildings marked as non-residential were assigned to a sub-category: commercial, industrial, or agricultural buildings (86%) and public buildings (14%). The total volume of non-residential buildings is 25 million m 3 . The mean volume of a non-residential building is about 13 thousand m 3 and the mean height is 13.5 m. The volume of individual buildings in the city and in conversion to a 100 × 100 m grid was presented in Figures 11  and 12.

Urban Volumetry and Building 3D Density Index (B3DI)
The volume calculation and 3D spatial analyses performed on the model of buildings from 2019 show that total building volume in Luxembourg City is about 56 million m 3 , which is 460 m 3 per resident. The mean volume of all types of buildings is 2349 m 3 per building. The tallest building (in the business center of Kirchberg) reaches 114 m in height (84,683 m 3 ), and the largest volume of a single building is 580 thousand m 3 in Gasperich district (Cloche d'Or Shopping Center). The other largest buildings in the city are hypermarkets in Kirchberg, European Union edifices, some banks, ArcelorMittal company, European School, and hospitals.
The volume in 2019 was calculated separately for the residential and non-residential buildings. The volume of all residential buildings in Luxembourg City is 30 million m 3 . The mean volume of a residential building is 1460 m 3 and the mean height is 12.5 m. Buildings marked as non-residential were assigned to a sub-category: commercial, industrial, or agricultural buildings (86%) and public buildings (14%). The total volume of non-residential buildings is 25 million m 3 . The mean volume of a non-residential building is about 13 thousand m 3 and the mean height is 13.5 m. The volume of individual buildings in the city and in conversion to a 100 × 100 m grid was presented in Figures 11  and 12.      (Table A2).  (Table A2).

Changes in City Volumetry over the Past 20 Years
The dynamics of change in building urban volumetry over the last 20 years in Luxembourg City was analyzed by comparing the Building 3D Density Index in grid of 100 m and at the district-level in 2001, 2010, and 2019 (Figures 14 and 15). Overall change in the past 20 years in Buildings 3D Density Index was +0.

Changes in City Volumetry over the Past 20 Years
The dynamics of change in building urban volumetry over the last 20 years in Luxembourg City was analyzed by comparing the Building 3D Density Index in grid of 100 m and at the district-level in 2001, 2010, and 2019 (Figures 14 and 15). Overall change in the past 20 years in Buildings 3D Density Index was +0.

Discussion
The spatiotemporal monitoring of changes in the volume of built-up areas and 3D modeling of buildings is a very actual topic that can reshape the perception of the urban landscape and the environment in which we live. Our research shows that the combination of up-to-date remote sensing data, spatial data, and archival aerial orthophotos offers an accurate estimation of changes that have occurred in urban volumetry over time and across space. Hence, our methodology allows supplementing incomplete time series data for exploring changes in 3D urban structure using ALS LiDAR point cloud and GEOBIA approach of the archival orthophoto maps. The concepts of the spatial index based on 3D information contained in the ALS point clouds allow for a synthetic representation of the spatial features such as buildings volume and 3D density. The automatic 3D modeling was used to calculate the volume, and by combining the results with data of demolitionreconstruction in plots of land (Housing Observatory project, LISER-Luxembourg), it was possible to use the models in previous years on unmodified buildings, where 3D information was not available.

Discussion
The spatiotemporal monitoring of changes in the volume of built-up areas and 3D modeling of buildings is a very actual topic that can reshape the perception of the urban landscape and the environment in which we live. Our research shows that the combination of up-to-date remote sensing data, spatial data, and archival aerial orthophotos offers an accurate estimation of changes that have occurred in urban volumetry over time and across space. Hence, our methodology allows supplementing incomplete time series data for exploring changes in 3D urban structure using ALS LiDAR point cloud and GEOBIA approach of the archival orthophoto maps. The concepts of the spatial index based on 3D information contained in the ALS point clouds allow for a synthetic representation of the spatial features such as buildings volume and 3D density. The automatic 3D modeling was used to calculate the volume, and by combining the results with data of demolitionreconstruction in plots of land (Housing Observatory project, LISER-Luxembourg), it was possible to use the models in previous years on unmodified buildings, where 3D information was not available.

Discussion
The spatiotemporal monitoring of changes in the volume of built-up areas and 3D modeling of buildings is a very actual topic that can reshape the perception of the urban landscape and the environment in which we live. Our research shows that the combination of up-to-date remote sensing data, spatial data, and archival aerial orthophotos offers an accurate estimation of changes that have occurred in urban volumetry over time and across space. Hence, our methodology allows supplementing incomplete time series data for exploring changes in 3D urban structure using ALS LiDAR point cloud and GEOBIA approach of the archival orthophoto maps. The concepts of the spatial index based on 3D information contained in the ALS point clouds allow for a synthetic representation of the spatial features such as buildings volume and 3D density. The automatic 3D modeling was used to calculate the volume, and by combining the results with data of demolition-reconstruction in plots of land (Housing Observatory project, LISER-Luxembourg), it was possible to use the models  [52]. Built volume was over 204 million m 3 , which indicates a much higher index (B3DI = 2.4) than in other European cities. The results show that in Luxembourg City the B3DI index and volume of buildings per resident are higher than most of the Cities in the UK and lower than German Cities (Figures 17 and 18). The high value of B3DI for Luxembourg City compared to the UK cities (city borders at municipality level) can be explained by the fact that Luxembourg is a relatively small-sized city (51 km 2 ). However, the separation of residential and non-residential buildings in other European cities is needed for better comparability.
Previous work [53,54], performed the detection of buildings using a similar GEOBIA approach, showed comparable accuracy to our results (KIA: 0.95 and 0.92). Tiwari et al. (2020) developed an object-based algorithm on orthophoto (GSD 0.25 m) to extract all buildings in a small city in Israel with a Kappa index of agreement of 0.95 (UA and PA of the building classes are 99.38% and 98.03%) [53]. Our methodology is not limited to aerial orthophotos, very high resolution (VHR) satellite imagery is often used for building detection and offers high flexibility when it comes to data available.  [55]. The growing number of the upcoming new Earth Observation missions and accompanying increase in VHRS imagery (e.g., Pléiades Neo [56] and WorldView Legion [57]) ensure an average daily revisit time and thus flexibility in terms of available data. Therefore, our approach does not have to depend on the available open data, if we consider the satellite image fee, we can get data from a specific day. Additionally, 3D data is also becoming more common and open access. The wall-to wall LiDAR ALS point clouds coverage is available for many European countries (e.g., in Denmark, Estonia, Finland, Latvia, Luxembourg, the Netherlands, Poland, Slovenia, Spain, Sweden, Switzerland, UK), but time series analysis is still not always possible due to some gaps in the data from several years.  Our study has some limitations. First, lack of ALS LiDAR data for the three years 2001, 2010, and 2019. For this reason, only in 2019, we were able to carry out 3D modeling of buildings. For detecting changes in the last 20 years, we used available data with various resolutions, in 2010 CIR and RGB orthophotos with 0.25 m, and in 2001 only RGB with 0.5 m spatial resolution, but this did not significantly affect the overall accuracy (the difference was 1.3%). To obtain the height of the buildings for the archival data, the Urban Atlas Buildings Height layer (2011) was used to compute the volume of the buildings that had been modified and could cause a significant demolition-  Our study has some limitations. First, lack of ALS LiDAR data for the three years 2001, 2010, and 2019. For this reason, only in 2019, we were able to carry out 3D modeling of buildings. For detecting changes in the last 20 years, we used available data with various resolutions, in 2010 CIR and RGB orthophotos with 0.25 m, and in 2001 only RGB with 0.5 m spatial resolution, but this did not significantly affect the overall accuracy (the difference was 1.3%). To obtain the height of the buildings for the archival data, the Urban Atlas Buildings Height layer (2011) was used to compute the volume of the buildings that had been modified and could cause a significant demolition- Our study has some limitations. First, lack of ALS LiDAR data for the three years 2001, 2010, and 2019. For this reason, only in 2019, we were able to carry out 3D modeling of buildings. For detecting changes in the last 20 years, we used available data with various resolutions, in 2010 CIR and RGB orthophotos with 0.25 m, and in 2001 only RGB with 0.5 m spatial resolution, but this did not significantly affect the overall accuracy (the difference was 1.3%). To obtain the height of the buildings for the archival data, the Urban Atlas Buildings Height layer (2011) was used to compute the volume of the buildings that had been modified and could cause a significant demolition-restoration error. For buildings that did not change between 2001-2019 and 2010-2019, LiDAR data was used to eliminate errors caused by lower data resolution. The comparison of the results (volume of buildings in 2001, 2010, 2019) obtained by data fusion was possible by using the historical data from the demolition-reconstruction layers, which contained information about buildings that had changed (with details of the date) or buildings which never changed. This allowed the use of LiDAR data wherever no changes were recorded. In 2001, 90.9% of the building volume could be calculated with LiDAR, as no changes were recorded for these buildings between 2001 and 2019. In 2010, as much as 94.5% of the volume of buildings could be calculated with LiDAR, as there was no change between 2010 and 2019 in buildings. Therefore, due to the availability of the demolition-reconstruction layer, the error in comparing the building volumes in 2001, 2010, and 2019 does not affect the results to a large extent. Second, we were unable to assess how index changes for residential buildings because in 2001 and 2010 no cadastral information of this type was available.
Future work will simulate urban densification and assess its socioeconomic, population, environmental impacts, and will monitor urban greenery and its relation to built-up areas. The debate about the urban growth of Luxembourg for the coming future should be open as soon as possible. As the different strategies elaborated 20 years ago concerning polycentric concentration or deconcentrated concentration [58] are not efficient enough and should be revised. Even if the third industrial revolution is coming quickly, it does not resolve all problems of a growing population, urban development, and densification, or big pressure on the land. Luxembourg City is still maintaining its principal role in the country and attracting people to live in the capital. The discussion surrounding the distribution of densities across the urban area can certainly be one of the most controversial issues in urban planning. The measure of the current situation and geovisualization of 3D changes over years assists the analysis of cities and can help urban planners for smart solutions in urban form for future developments. From this perspective, defining volumetric ratios that explain the dynamics of urban areas can bring more clarity to the city's evolution debate and management of urban space.

Conclusions
The methodology presented in is paper that employed ALS LiDAR point cloud and archival aerial orthophotos allows quantifying changes in urban volumetry over time and space. Further, our methodology can provide monitoring based on buildings 3D density index as an effective tool that can be used to characterize the built-up configurations in our cities. Due to the increasingly widespread use of LiDAR technologies, point clouds present valuable spatial information about buildings in a clear and synthetic way, thus being a tool for extracting 3D features difficult to determine by traditional methods. Our findings demonstrate the value of ALS LiDAR data for monitoring effectively the city volume density at some specific dates. Despite airborne laser scanning data is widely recognized as a highly precise data source, time-series data are not available in many cities to detect temporal changes. Therefore, we proposed the use of GEOBIA classification based on archival orthophotos. The ruleset produced for this study to extract buildings can be used in other cities based on CIR aerial photos or VHR satellite images supported by other technology like SLS-satellite laser scanning (e.g., GEDI, ICESat-2 missions NASA) or radar technology (e.g., TerraSAR-X; IceEye-1).
We believe that this paper provides an important methodological approach for urban planners and city authorities. Tracking changes in city volume and 3D index are critical for urban policies that concern equality in the distribution of benefits, e.g., equal access to green spaces and amenities, for city residents. Due to the increasing availability of VHR satellite images and the potential of using unmanned aerial vehicles (UAVs), or the rapidly evolving high-altitude platform station (HAPS) systems for cost-effective mapping of urban areas in 2D/3D, it could ensure constant access to quickly obtain data in an area of interest. Acknowledgments: We thank the Housing Observatory, a collaboration between the Luxembourg Housing Ministry and the Luxembourg Institute for Socio-Economic Research, for providing access to their data.

Conflicts of Interest:
The authors declare no conflict of interest.
Appendix A Table A1. Volume of buildings and Building 3D Density Index (B3DI) for residential (R) and non-residential (NR) buildings in Luxembourg City in 2019.