Monitoring of Urban Sprawl and Densification Processes in Western Germany in the Light of SDG Indicator 11.3.1 Based on an Automated Retrospective Classification Approach

By 2050, two-third of the world’s population will live in cities. In this study, we develop a framework for analyzing urban growth-related imperviousness in North Rhine-Westphalia (NRW) from the 1980s to date using Landsat data. For the baseline 2017-time step, official geodata was extracted to generate labelled data for ten classes, including three classes representing low, middle, and high level of imperviousness. We used the output of the 2017 classification and information based on radiometric bi-temporal change detection for retrospective classification. Besides spectral bands, we calculated several indices and various temporal composites, which were used as an input for Random Forest classification. The results provide information on three imperviousness classes with accuracies exceeding 75%. According to our results, the imperviousness areas grew continuously from 1985 to 2017, with a high imperviousness area growth of more than 167,000 ha, comprising around 30% increase. The information on the expansion of urban areas was integrated with population dynamics data to estimate the progress towards SDG 11. With the intensity analysis and the integration of population data, the spatial heterogeneity of urban expansion and population growth was analysed, showing that the urban expansion rates considerably excelled population growth rates in some regions in NRW. The study highlights the applicability of earth observation data for accurately quantifying spatio-temporal urban dynamics for sustainable urbanization and


Introduction
An urban sprawl that often results in the conversion of natural vegetation cover and the expansion of the spatial footprint of impervious surfaces is one of the most striking landuse changes. Urbanization, which is caused by economic development and continuous population growth, as well as changing lifestyles, and residential and retail uses within the urban fringes is one of the most important drivers of terrestrial change [1,2]. As the global trend of urbanization is expected to further increase the urban share of the world population from about 55 percent in 2018 to 68 percent by 2050 [3], understanding transformations of urban areas, and navigating those transformations towards more sustainable path-ways is and have global coverage. Furthermore, at a European level, Copernicus high-resolution layers provide several thematic layers, including the imperviousness level [33]. Although the data have high spatial resolution and relevant thematic context, it is only available for five-time steps (2006)(2007)(2008)(2009)(2010)(2011)(2012)(2013)(2014)(2015)(2016)(2017)(2018), thus, limiting the analysis of densification processes to a smaller, recent time interval. Besides the limitation of different continental to global land cover products, with the limited coverage over time, available long-term products often lack thematic context and present data that support the analysis of single class (built-up areas/urban/settlement) changes [9]. Therefore, within-class densities and their dynamics are often not assessed.
Urbanization has been one of the most prominent change processes in Central Europe coupled with simultaneous sociodemographic changes [34,35]. Around 73% of the European population lives in cities, which is projected to reach 82% by 2050. The expansion of built-up areas was observed in most regions of Europe, even in areas where the population has declined [36]. In Germany, the land consumption and reallocation of natural and agricultural land to built-up and transportation-related areas exceeds 60 ha per day [24,37,38]. Although in recent years, artificial development in Germany has shown decrease, especially for diffuse residential sprawl, the area around industrial and commercial sites has increased significantly during 2006-2012, which was followed by decline in the period 2000-2006 [39]. Therefore, the need to monitor and accurately quantify these changes increases. In Germany, in general, but also in the German federal state of North Rhine-Westphalia (NRW), parallel development has been observed, with a distinct growth of settlements and transportation infrastructure. In the national sustainability strategy, presented in 2002, the 30-ha goal was introduced by the German government aimed at reducing Germany's land consumption to a maximum of 30 ha per day until 2030 [4]. Although the daily land consumption has declined during the past years, it was still at 10 ha per day in NRW [40]. Furthermore there were changes applied to the cadastral data recording, which made the year-to-year comparison challenging. Alongside the introduction of ALKIS (Official Real Estate Cadaster Information System) as the official real estate cadaster information system, all German federal states have replaced their previous records. In NRW, the last cadastral authorities switched over in 2015. With the changes in some of the categories of cadastral system [41], the areas belonging to the category "urban" have not changed much in terms of quantity, but a lot in terms of location.
Based on this, the objectives of this study are; first, to implement an automated retrospective classification approach for observing land use and imperviousness dynamics during the last 30 years; second, to analyze the spatial trends of urban sprawl in NRW, and; third to assess the intensity of densification and population dynamics to support the monitoring of SDG indicator 11.3.1.

Study Area
The study area covers the German federal state of North-Rhine Westphalia (NRW) (Figure 1). It is the most densely populated state in Germany, counting 526 persons per km 2 [42]. NRW has 396 municipalities, including 22 independent cities. Different LC/LU characterizes this area. Around half of the state area is covered by seminatural arable lands, and forests cover 20%. Built-up areas comprise around 20% of the area. The area includes the Rhine-Ruhr Metropolitan Area (RRMA) (Figure 1), which is the largest metropolitan region in Germany with over 10 million inhabitants [42], and which has the largest area of continuously high urban sprawl [36]. Although the population has decreased during the last decade in NRW, the population has grown from 1985 to 2017 by more than 7% [42]. Unlike the trends within the RRMA mentioned above, some areas, such as Cologne and Düsseldorf, which are the largest cities in NRW, had a continuous population increase. During the same period, the land area within NRW dedicated to settlement and transportation infrastructure increased [40]. Accordingly, the study area is heterogeneous and challenging when it comes to monitoring urban change dynamics, Remote Sens. 2021, 13, 1694 4 of 17 along with the highly urbanized areas, it comprises rather rural areas that are currently experiencing high urbanization rates. est area of continuously high urban sprawl [36]. Although the population has decreased during the last decade in NRW, the population has grown from 1985 to 2017 by more than 7% [42]. Unlike the trends within the RRMA mentioned above, some areas, such as Cologne and Düsseldorf, which are the largest cities in NRW, had a continuous population increase. During the same period, the land area within NRW dedicated to settlement and transportation infrastructure increased [40]. Accordingly, the study area is heterogeneous and challenging when it comes to monitoring urban change dynamics, along with the highly urbanized areas, it comprises rather rural areas that are currently experiencing high urbanization rates.
As the primary source of reference information, official geodata from the German Authoritative Topographic-Cartographic Information System (ATKIS) was extracted, in order to generate labeled data for 10 classes. For the classes water, coniferous, mixed, and
As the primary source of reference information, official geodata from the German Authoritative Topographic-Cartographic Information System (ATKIS) was extracted, in order to generate labeled data for 10 classes. For the classes water, coniferous, mixed, and deciduous forests, arable land, grassland, and open-pit mine samples were extracted from the Digital-base Landscape Model (BASIS-DLM) for 2017. In ATKIS DLM, the topographic features of the landscape are described by point, line, and shape in vector format with information describing residential area, traffic, vegetation, water, administrative area, and relief forms, with each category containing additional levels (e.g., different types of roads) [16]. ATKIS data are available since 1990, and the data are continuously updated. Unfortunately, it is not possible to reconstruct past time steps from the official release. Therefore, reference data exists only for 2017 in this study. In order to further interpret the change patterns in the study area, population data for administrative regions in NRW were extracted and analyzed [42].

Methods
The proposed method can be summarized in three major steps ( Figure 2). The first step is data preparation and pre-processing of Landsat images. The second step is the baseline classification, implemented for the year 2017, which included the generation of labeled data for three imperviousness classes as well as the classification with a land-use map as an output. In the third step, we perform a per-pixel change detection analysis to detect pixels that are stable between two-time steps, and thus, can be used as an input for classification. This procedure was iterated for the years that preceded 2015, resulting in five change detection analysis. Detailed procedures of the approach are described below.
information describing residential area, traffic, vegetation, water, administrative area, and relief forms, with each category containing additional levels (e.g., different types of roads) [16]. ATKIS data are available since 1990, and the data are continuously updated. Unfortunately, it is not possible to reconstruct past time steps from the official release. Therefore, reference data exists only for 2017 in this study. In order to further interpret the change patterns in the study area, population data for administrative regions in NRW were extracted and analyzed [42].

Methods
The proposed method can be summarized in three major steps ( Figure 2). The first step is data preparation and pre-processing of Landsat images. The second step is the baseline classification, implemented for the year 2017, which included the generation of labeled data for three imperviousness classes as well as the classification with a land-use map as an output. In the third step, we perform a per-pixel change detection analysis to detect pixels that are stable between two-time steps, and thus, can be used as an input for classification. This procedure was iterated for the years that preceded 2015, resulting in five change detection analysis. Detailed procedures of the approach are described below.

Data Processing
To overcome the impact of cloud and cloud shadow, we used the FMASK algorithm for Landsat cloud masking [43]. Images with more than 80% clouds were excluded from further analyses. Following the cloud mask, we calculated spectral indices, namely the Normalized Difference Vegetation Index (NDVI) [44,45] and the Normalized Difference

Data Processing
To overcome the impact of cloud and cloud shadow, we used the FMASK algorithm for Landsat cloud masking [43]. Images with more than 80% clouds were excluded from further analyses. Following the cloud mask, we calculated spectral indices, namely the Normalized Difference Vegetation Index (NDVI) [44,45] and the Normalized Difference Water Index (NDWI) [46], as well as Tasseled cap wetness, greenness and brightness (Table 1) [47]. We implemented a pixel-based compositing approach before the classification procedure and calculated temporal metrics [48][49][50][51] (e.g., 10th, 25th, 50th, 75th, 80th and 90th percentiles of annual time series, as well as a minima and maxima for both annual and growing season time series), which resulted in composites for each target year. We extracted eleven types of image features from each compositing result, including six surface reflectance (green, blue, red, NIR, SWIR1 and SWIR2) bands and the spectral indices (NDVI, NDWI, TC greenness, wetness, brightness). As there were differences in the coverage of cloud free pixels during different years, the imagery was reduced to the growing season (e.g., 1985, 2010). Therefore, the total number of input features that was used in a single year classification was 99. Since it was not possible to generate a gapless dataset for 1985 the imagery from 1986 was used Remote Sens. 2021, 13, 1694 6 of 17 as well. Due to the low quality of the composite in 1995, this timestep was excluded from further analysis.

Land Cover and Land Use Classification
To create the land-use maps for seven timesteps, the Random Forest (RF) [52] classifier was selected, as it is known to be robust for large-area land cover classification [53][54][55], and specifically useful with a large number of features [2]. For all classifications, the RF classifier was implemented with 500 trees and the square root of the total number of input variables as the number of variables to split at each node. The analysis was mainly carried out with R statistical software (randomForest package [56]) [57].
As one of the aims of the study was to estimate the changes with the differentiation of several imperviousness classes, an additional step was carried out for the generation of samples for the baseline year. First, regression was run based on ATKIS data. The data in ATKIS DLM, which represents residential areas, as well as roads were used. On a generated 30 m grid, the percentage of the cover of each grid cell with different built-up surfaces was estimated. Then RF regression was used to create a continuous layer at 30 m resolution showing the degree of imperviousness. Based on this layer, three classes of high, middle and low imperviousness samples were separated based on the percentage of the impervious area, corresponding to higher than 80%, 40-80% and less than 40% impervious surfaces per pixel. [22,58,59]. Similar classification schemes can be observed in other products, such as Corine Land cover, where continuous urban fabric class corresponds to areas with more than 80% imperviousness [60]. This approach has been shown to be suitable for deriving consistent classification results in urban areas with the required high precision. Furthermore, it identifies which regions experience scattered, diffusive, and dispersed urban growth. Combined with the thematic information, representing the other classes, RF classifier was run for 2017 with the use of 20,000 samples. The classes represent both environmentally and socio-economically important classes that can also capture the changes during the time span of the study. These classes include vegetated areas, such as forests, croplands and grasslands, as well as impervious areas. We separated water class to capture water bodies in the study area. As NRW includes highly industrial areas with formerly intensive coal mining, we included a separate class for open pit mine, in order to reduce the misclassification of these areas with urban classes.
Following the generation of accurate output for the baseline year, the Iteratively Reweighted Multivariate Alteration Detection (IR-MAD) was used [61,62], in order to derive the areas with no change, which are the basis for the extraction of reference samples for the classification of the preceding timestep. The IR-MAD transformation identifies suitable invariant pixels for radiometric normalization with a multivariate change detection algorithm [41] that iteratively increases the weight of non-change observations, and thereby, enhances the discrimination of substantial spectral changes between images from different time periods [63,64]. For our analysis, only the median composites of spectral bands and indices (in total 11 bands) were used in the change detection procedure, and a threshold of 0.25 was chosen for the chi-square image based on the IR-MAD algorithm to extract the unchanged area.
Finally, based on baseline classification and change detection, we applied a stepwise classification framework to former years by selecting new training samples using spectral Remote Sens. 2021, 13, 1694 7 of 17 information from unchanged pixels. Due to the changes between the timesteps and the subsequent decrease of the number of unchanged pixels, the number of samples decreased with each timestep. We assessed the accuracy of the final maps based on confusion matrices [65]. For the classification reference, we split the data randomly proportionally for training and validation. From these, we derived overall accuracy (OA), as well as producer's (PA) and user's (UA) accuracy.
The baseline map of 2017 was compared to the two available products at national level [66,67].

Intensity Analysis and Comparison with Population Data
To further understand the underlying processes of land-use change and urban growth and to evaluate whether the growth was sustainable, intensity analysis was applied. In general, the intensity analysis evaluates the predicted uniform rate of change compared to the observed rate of change among categories [68,69], characterizing the process of land use change. It also evaluates the land use change by setting the stability. Although the stable land use conversion may not always result in sustainable development, it can still be an indicator and basis for it [70].
Furthermore, the population data at a district level were used to compare the rates of urban growth and population change and calculate the ratio of land consumption rate (LCR) to population growth rate (PGR). It is based on the ratio of Land Consumption Rate and Population Growth Rate (LCRPGR) (Equations (1)-(3)) [8,51,71], where Urb t and Urb (t+N) is the areal extent of the land consumed at the initial reference year t and at the final year t + N; Pop t and Pop (t+N) input the total population of the spatial unit at the initial reference year and at the final reference year, y is the number of years between t and t + N.

Land-Use Classification and Changes between 1985-2017
The classified maps differentiated between ten classes: Arable land, grassland, deciduous forest, coniferous forest, mixed forest, water, open-pit mines as well as low, medium, and high levels of imperviousness. The baseline classification of 2017 as well as the levels of imperviousness based on RF regression are shown in Figure 3. The resulting maps had acceptable quality, exceeding 75% of overall accuracy. Table 2 presents the OA, UA and PA for all time steps and UA and PA for the three impervious classes as the main target of the study in particular. The baseline map of 2017 was compared to the two available products at national level [66,67]. The maps had high agreement, especially for the classes of high-and mid-levels of imperviousness with the interim land cover product for class build up [66]. Whereas, all imperviousness classes had high agreement with class artificial land in the case of a second product [67]. of the study in particular. The baseline map of 2017 was compared to the two available products at national level [66,67]. The maps had high agreement, especially for the classes of high-and mid-levels of imperviousness with the interim land cover product for class build up [66]. Whereas, all imperviousness classes had high agreement with class artificial land in the case of a second product [67].  We observed an accuracy decrease for the retrospective classification. Several factors can cause this accuracy decrease. First, due to the change detection, the number of unchanged pixels that can be used for training and validation decreases over the years. In the case of our study the total number of reference data decreased by 50% from 2017 to 1985. At the same time, we observe a large drop in accuracies for the time step of 1990, which has the biggest interval for change detection. Furthermore, the resolution of Landsat (30 m), which is higher than the size of many urban objects, can be the cause of mixed  We observed an accuracy decrease for the retrospective classification. Several factors can cause this accuracy decrease. First, due to the change detection, the number of unchanged pixels that can be used for training and validation decreases over the years. In the case of our study the total number of reference data decreased by 50% from 2017 to 1985. At the same time, we observe a large drop in accuracies for the time step of 1990, which has the biggest interval for change detection. Furthermore, the resolution of Landsat (30 m), which is higher than the size of many urban objects, can be the cause of mixed pixels and result in misclassifications [15]. We observed this issue, especially for low imperviousness levels, which exhibited lower PA and UA.
In contrast, the high impervious class was the most accurately classified class, with accuracies exceeding 80%. Another reason for the drop in accuracy can be the reduced data availability before the launch of Landsat-7 in 1999. Among the input variables, NDVI and NDWI, TC brightness minima and maxima, as well as green, red, NIR, and SWIR median and 90 th percentile were among the essential input bands for classification. The combination of these metrics provided information on three imperviousness classes, helping to minimize misclassifications by minimizing seasonal effects and changes.
For 32 years, urban expansion resulted in large-scale increases in other land-use types transforming into a different level of imperviousness. The transitions of land-use are shown in Figure 4. The year 2005 was selected as a middle timestep, as it followed the implementation of the national sustainability strategy in 2002. It can be seen that arable land, forests, and grassland are the primary land resources for urban expansion. Transitions cause the increase in the low imperviousness class from arable land, grassland, and small Remote Sens. 2021, 13, 1694 9 of 17 changes in the forest classes. Simultaneously, a high level of dynamics can be observed between the three imperviousness classes, indicating a transition from low to middle and high level of imperviousness. The analysis of the maps shows that impervious surfaces comprised about 552,464 ha in 1985. Three decades later, impervious surfaces increased by more than 167,000 ha in 2017. The multiyear classification can contribute to quantification of the changes in areas where time series of reference data is not available. Specifically, according to the study of the Federal Institute for Research on Building, Urban Affairs and Spatial Development concluded that an accurate monitoring of land consumption during the conversion period from ALK/ALB to ALKIS, i.e., before and after 2016, is not less than difficult [72]. Three alternative databases were compared but none of them fulfilled the criteria of a possible replacement. ATKIS and LBM-DE (Digital land cover model for Germany) data are often outdated. Accordingly, the suggested approach provides insights on land consumption in NRW with a high resolution and temporal consistence. This makes the data of high value for researchers and the official statistical offices.  Because the degree and dynamics of imperviousness can be used as critical indicators for assessing the sustainability of land-use changes, we analyzed the intensity level of the changes ( Table 3). The interval level of intensity shows the size and speed of change across different time intervals. The category level of the intensity analysis gives information about which land use category is relatively dormant or active during a particular time interval. During our study period, the 2005-2017 interval was fast compared to uniform growth, while the earlier 1985-2005 interval was slower than uniform growth.  Because the degree and dynamics of imperviousness can be used as critical indicators for assessing the sustainability of land-use changes, we analyzed the intensity level of the changes ( Table 3). The interval level of intensity shows the size and speed of change across different time intervals. The category level of the intensity analysis gives information about which land use category is relatively dormant or active during a particular time interval. During our study period, the 2005-2017 interval was fast compared to uniform growth, while the earlier 1985-2005 interval was slower than uniform growth. At the category change intensity of grassland, water, forests, as well as low and high imperviousness exhibited greater than the average intensity during 1985-2005, indicating active change. The changing intensity of arable land was less than the average intensity, indicating rather dormant change. Between 2005 and 2017, open-pit mine and the middle imperviousness class exhibited greater than average change intensity. The rate at which other land use categories are converted to imperviousness is higher than the transition rate of imperviousness classes to other land use categories.
Based on the results of our intensity analysis, between 1985 and 2005 the high and middle imperviousness classes were targeted from the low imperviousness class. In the case of 2005-2017, transitions from low to high impervious level were largely avoided.
The gain intensity of impervious areas ( Figure 5) in the first interval was higher than the annual uniform intensity of 1.71 for the low and the high imperviousness classes, thus making them an active land use category. During this period, other classes, such as arable land or open-pit mine were dormant, with gain intensity lower than uniform. For the second interval, compared to the uniform change rate of 2.72, the middle imperviousness level was an active class. Coupled with the analysis at a community level (Figure 6), these metrics can reveal the patterns of urban change. When considering the net urban growth, the more rural eastern region in NRW had the highest relative gains during the last 30 years ( Figure 6). The gain intensity of impervious areas ( Figure 5) in the first interval was higher than the annual uniform intensity of 1.71 for the low and the high imperviousness classes, thus making them an active land use category. During this period, other classes, such as arable land or open-pit mine were dormant, with gain intensity lower than uniform. For the second interval, compared to the uniform change rate of 2.72, the middle imperviousness level was an active class. Coupled with the analysis at a community level (Figure 6), these metrics can reveal the patterns of urban change. When considering the net urban growth, the more rural eastern region in NRW had the highest relative gains during the last 30 years ( Figure 6).    Figure 7 shows the relationship between the LCR and PGR trends presented in the LCRPGR. LCRPGR estimates for the period 1985-2005 are positive in most municipalities. This shows parallel increase of both population growth and soil consumption rate. In various areas, the ratio varies between zero and one. In these municipalities, a simultaneous increase in the population growth rate can be observed. In the period 2005-2017, the proportion of municipalities that presented a negative estimate of LCRPGR was higher, showing the soil consumption increase combined with urban population reduction. This indicates that land use changes and population growth are not always in parallel, underlining the detachment of population growth as the main underlying cause of urban growth. Municipalities with high LCRPGR can be found in the peri-urban surroundings of the Rhine-Ruhr metropolitan region and in the rural areas in the Northern part. As a result, 1985-2005 was a period of rural growth (according to the population pattern in Figure 7a), while 2005-2017 was a period of urbanization (i.e., growth of urban population, rural exodus, Figure 7b).  Figure 7 shows the relationship between the LCR and PGR trends presented in the LCRPGR. LCRPGR estimates for the period 1985-2005 are positive in most municipalities. This shows parallel increase of both population growth and soil consumption rate. In various areas, the ratio varies between zero and one. In these municipalities, a simultaneous increase in the population growth rate can be observed. In the period 2005-2017, the proportion of municipalities that presented a negative estimate of LCRPGR was higher, showing the soil consumption increase combined with urban population reduction. This indicates that land use changes and population growth are not always in parallel, underlining the detachment of population growth as the main underlying cause of urban growth. Municipalities with high LCRPGR can be found in the peri-urban surroundings of the Rhine-Ruhr metropolitan region and in the rural areas in the Northern part. As a result, 1985-2005 was a period of rural growth (according to the population pattern in    Figure 8 shows LCR and PGR, extracted for municipalities, RRMA, independent cities and non-metropolitan regions. Land consumption mainly occurs in non-independent cities, i.e., rural areas, and has accelerated in some rural municipalities. Population growth was slower in RRMA than the rest in 1985-2005 but is equal in 2005-2019. However, PGR has generally decreased. The average land consumption has lowered slightly for both RRMA and other municipalities.

Land Consumption and Population Dynamics
Using three imperviousness classes, distinct patterns have shown the densification processes and conversions from low and middle levels to high imperviousness from 1985 to 2017 (Figure 9). This is an advantage compared to one single urban class in the classification (Figure 6), which can omit patterns of the densification processes that can act as important indicators of sustainable growth. In general, densification is preferred over expansion and the development of new building areas since forests and agricultural areas can be preserved. Nevertheless, continuous increases in densification hold challenges, such as in logistics and urban climate [40]. This is the type of urban growth, which is the hardest to monitor with binary classification. The observation of the class low imperviousness and the transformation of that class to higher imperviousness levels provide the chance to observe densification of metropolitan regions, but also observe modification of land use for further ecosystem-oriented analyses. One of the main limitation of the study is the inaccuracies for the low imperviousness class, but they are still acceptable, considering the heterogeneity of the class and the spatial resolution of the dataset used. Nevertheless, our study allows land consumption information to be updated at a spatial scale adequate for land use planning at an administrative level. As opposed to existing datasets, our results reveal spatial patterns of densification and intensification within urban agglomerations. The approach considered the specific German situation, where the often conflicting interests in land use planning lead to small-scale changes and distinct patterns of densification in urban centers, and extensive land use conversion in the urban fringe. These processes are neither adequately addressed in existing global products nor in other approaches of urban change monitoring. As opposed to most published papers that often focus on urban areas with high growth rates, our study covers a heterogeneous region with predominantly small-scale urbanization, resulting in densification rather than more areal growth. We provide a methodology that allows us to detect and quantify these specific changes, which can have fundamental implications for socio-economic processes. As a result, the applied methodology contributes to urban mapping by offering a reliable way to account for changes with limited data and further extends the approach for SDG indicator estimation. Using three imperviousness classes, distinct patterns have shown the densification processes and conversions from low and middle levels to high imperviousness from 1985 to 2017 (Figure 9). This is an advantage compared to one single urban class in the classification (Figure 6), which can omit patterns of the densification processes that can act as important indicators of sustainable growth. In general, densification is preferred over expansion and the development of new building areas since forests and agricultural areas can be preserved. Nevertheless, continuous increases in densification hold challenges, such as in logistics and urban climate [40]. This is the type of urban growth, which is the hardest to monitor with binary classification. The observation of the class low imperviousness and the transformation of that class to higher imperviousness levels provide the chance to observe densification of metropolitan regions, but also observe modification of land use for further ecosystem-oriented analyses. One of the main limitation of the study is the inaccuracies for the low imperviousness class, but they are still acceptable, considering the heterogeneity of the class and the spatial resolution of the dataset used. Nevertheless, our study allows land consumption information to be updated at a spatial scale adequate for land use planning at an administrative level. As opposed to existing datasets, our results reveal spatial patterns of densification and intensification within urban agglomerations. The approach considered the specific German situation, where the often conflicting interests in land use planning lead to small-scale changes and distinct patterns of densification in urban centers, and extensive land use conversion in the urban fringe. These processes are neither adequately addressed in existing global products nor in other approaches of urban change monitoring. As opposed to most published papers that often focus on urban areas with high growth rates, our study covers a heterogeneous region with predominantly small-scale urbanization, resulting in densification rather than more Although the applied approach showed potential for urban growth analyses, further investigations can be improved by integrating other datasets, such as night time lights. The addition of SAR data from TanDEM-X and Sentinel-1 can further delineate the urban Although the applied approach showed potential for urban growth analyses, further investigations can be improved by integrating other datasets, such as night time lights. The addition of SAR data from TanDEM-X and Sentinel-1 can further delineate the urban structures and accurately reveal densification processes [73]. This study focused on the classification of imagery at a regional scale. In future studies, Landsat data can be complemented with Sentinel-2 and applied at a larger scale to derive LU/LC information at both high spatial and thematic resolution. Future research should focus on method transferability on multiple levels. The approach can be implemented and applied over the larger scale with the availability of the change detection algorithm on Google Earth Engine [74] and available Landsat time series.

Conclusions
This study presents an approach for automated retrospective LU/LC classification with the use of the spectral information of unchanged pixels. The method's performance indicates that it is a promising approach for quantifying the changes in land use, including the derivation of several levels of impervious surfaces, with the use of reference data only for the baseline timestep. The analysis indicates a continuously growing imperviousness area during 1985-2017 in our study area. Further analysis of intensity showed the importance of identifying three imperviousness levels, which can better describe the densification processes of the urban areas in the study site. Using the approach, land consumption information was updated at a spatial scale, adequate for land use planning at an administrative level. By conducting an intensity analysis and the integration of population data, the spatial heterogeneity of urban expansion and population growth was assessed. We observed that the urban expansion rates considerably increased population growth rates in some of the regions in NRW, indicating more expansive urban growth patterns. Land consumption has accelerated in at least some of the rural municipalities. Furthermore, it is much lower in the RRMA as compared to the remaining municipalities. The analysis of LCRPGR showed that, especially for 2005-2017, land use change and population growth were not simultaneous. The former relationship of population growth and urban growth can be dismissed for the study area in Western Germany. This is an important finding as it proves that when it comes to SDG 11.3.1 urban, as well as rural areas are not maintaining an efficient land use policy, so that a sustainable settlement development until 2030 is hard to be seen in the Federal State of North Rhine-Westphalia. Especially the periurban regions are experiencing an inefficient settlement growth. The study highlights the applicability of earth observation data for accurately quantifying urban change for further sustainable targeted planning practices.
Author Contributions: Conceptualization, A.R., G.G. and F.T.; funding acquisition and supervision A.R. and C.J.; formal analysis, G.G., F.T., C.O. and A.R.; data curation, C.O., S.S. and B.T. All authors have read and agreed to the published version of the manuscript.

Funding:
The project "Town and Country in the Flow-Network for the Creation of a sustainable Climate Landscape" (KlimNet) was funded with means of the Federal Ministry for the Environment, Nature Conservation and Nuclear Safety following a decision of the German Bundestag (funding code 67DAS098 ABC) within the German Strategy for Adaptation to Climate Change (Deutsche Anpassungsstrategie an den Klimawandel, DAS).

Data Availability Statement:
The raw data used for this study (Landsat surface reflectance) can be acquired via USGS (https://earthexplorer.usgs.gov/ (accessed on 1 March 2021)). The reference data are subject to permission of third party (German Authoritative Topographic-Cartographic Information System). Additional demographic data were accessed from State Office Information and Technology North Rhine-Westphalia https://www.it.nrw/ (accessed on 1 March 2021).