Remote Sensing Based Methodology for the Quick Update of Population Exposed to Natural Hazards

: The assessment of the number of people exposed to natural hazards, especially in countries with strong urban growth, is difﬁcult to be updated at the same rate as land use develops. This paper presents a remote sensing based procedure for quick updating the assessment of the population exposed to natural hazards. A relationship between satellite nightlights intensity and urbanization density from global available cartography is ﬁrst assessed when all data are available. This can be used to extrapolate urbanization data at different time steps, updating exposure each time new nightlights intensity maps are available. As reliability test for the proposed methodology, the number of people exposed to riverine ﬂood in Italy is assessed, deriving a probabilistic relationship between DMSP nightlights intensity and urbanization density from GUF database for the year 2011. People exposed to riverine ﬂood are assessed crossing the population distributed on the derived urbanization density with ﬂood hazard zones provided by ISPRA. The validation on reliable exposures derived from ISTAT data shows good agreement. The possibility to update exposure maps with higher refresh rate makes this approach particularly suitable for applications in developing countries, where exposure may change at sub-yearly scale.


Introduction
The aim of this paper is to present a general, simple and remote sensing based procedure for the quick update of the population exposed to natural hazards in rural and urban areas. This procedure explores the potential of the combination of the high refresh rate of satellite nightlights intensity products with the high precision of urban cartography. The final product is a map of exposure which can be easily and quickly updated, with the refresh rate typical of remote sensing products, to be used in developed but especially in developing countries, where exposure may change at sub-yearly scale.
The impact of natural disasters such as floods, earthquakes or cyclones is growing significantly in recent decades [1]. This can be related to the growth of all components of the risk equation: hazard impacted by climate change on one hand and exposure and vulnerability affected by the rapid economic growth and changes in land use of many areas of the world on the other. While the evaluation of the hazard component of the risk equation is a quite well established and consolidated procedure (see, e.g., [2] for floods and [3] for earthquakes) and global hazard maps are produced for several hazards (see, e.g., [4] and [5] for floods and [6] for earthquakes), the assessment of both the built-up areas and the number of people exposed to natural risks is a complex process, which needs information that should be updated often, especially in areas with strong economic and urban development (see e.g. [7]). Nevertheless, this information is not easily available and even more difficult to update at the same rate of land use development. Currently available satellite based exposure products, such as the GUF database [8][9][10] or the ESM [11], are very detailed and reliable, designed for use even on a fine scale, e.g. for demographic monitoring, analysis on the spatial structure of built-up areas, etc. [9]. Nevertheless, the level of processing required to release a final reliable version of these products lasts for several years, making their refresh rates significantly lower than the urban growth rate of many areas of the world. It is worth mentioning here the provisions of the Sendai Framework for Disaster Risk Reduction [12], which requires a frequent update (e.g. yearly) of information on the population or assets exposed at national or provincial scale, to be used for the monitoring of some of its indicators, such as "Number of directly affected people attributed to disasters, per 100,000 population" (indicator B-1) or "Direct economic loss attributed to disasters in relation to global gross domestic product" (indicator C-1).
To provide a solution to the need of highly refresh rate exposure data, it would be possible to link the cartography of exposure to other easily available information that can be updated with a higher refresh rate. Satellite-based nightlights have been used by various authors to map global urban extent and its dynamic [13][14][15][16][17][18] or as proxy for population density [19][20][21], economic activity [22][23][24], and other. Results provided by Bagan and Yamagata [20] show that correlation between population density and nightlights intensity becomes significant for values of DMSP nightlights intensity above 61-62 (97-98 % of the maximum intensity). Ceola et al. [25] use the Version 4 DMSP-OLS Nighttime Lights Time Series to derive information about the exposure of population to the flood risk [26], linking the variation of the light intensity from 1992 to 2012 to the proximity of a river network, identified using the USGS HydroSHEDS database [27]. However, the temporal evolution of the exposure does not only depend on the intensity of the light in a defined area around the watercourses, but it is also linked to the evolution of the extension of the urbanized area (not entirely represented by the pixels surrounding rivers only) and to its intersection with the areas exposed to the flood hazard.
This paper presents and validate an experimental methodology to reconstruct the urbanization density and the exposure of built-up areas and people, to be used for aggregated (e.g. municipal, district, country level) risk and direct damage estimations. The method is based on the identification of a relationship between nightlights intensity derived from satellite sensors and urbanization density derived from global available cartography, starting from the assumption that the latter can be considered as a proxy to quantify exposure in hazard-prone areas. This relation can be assessed at a specific time when both nightlights intensity and urbanization maps are available, and then can be potentially used to extrapolate the urbanization at different time steps, when only nightlights are available. In this way, the cartography of exposure can be quickly updated each time a new nightlights intensity map is available.
To test the methodology, we propose a reconstruction of the spatial distribution of the population exposed to riverine flood in Italy based on a probabilistic relationship between the intensity of the nightlights provided by the DMSP database and the density of urbanization estimated from the GUF database. The probabilistic relation is obtained by fitting a calibration curve on the GUF and DMSP data in three different regions, for the reference year 2011. The obtained curve is then reliably used to redistribute at finer scales the regional population provided by ISTAT, outside the calibration regions, for the whole Italy. As a final step, people exposed to riverine flood is evaluated crossing the obtained distributed population with riverine flood hazard zones provided by ISPRA. The validation of the results is made at municipal scale on the reference exposed population derived from the independent database provided by ISTAT. It shows a good agreement.
The paper is organized as follows: in Paragraph 2 the methodology and the databases used for the analysis are described. Paragraph 3 reports the results obtained applying the proposed methodology for the assessment of the number of people exposed to riverine flood hazard in Italy, together with a validation of the obtained output at municipal scale. Paragraphs 4 and 5 report the discussion of the results and the related conclusions, respectively.

Materials and Methods
The proposed methodology allows to produce updated maps of built-up areas and population at country scale using readily available satellite-based information on small portions of a country (calibration regions). These maps, if combined with hazard maps, can be used for aggregated 3 of 15 (e.g. municipal, district, country level) risk and direct damage estimations. In order to apply the methodology nightlights intensity data, urbanization maps and a population estimation at least at national scale are required.
More specifically, the present approach starts from the identification of a relationship between intensity of nightlights derived from satellite sensors and urbanization density derived from global available urbanization maps, starting from the assumption that the latter can be considered as a proxy to quantify exposure in hazard-prone areas. This relation, named as urban density calibration curve, can be assessed for one or more calibration regions inside the country or region under study and derived at a specific time, when both nightlights intensity and urbanization maps are available. Each time a new nightlights intensity map is available, the resulting calibration curve/s (derived inside the calibration region/s) can be applied to derive an updated urban density map at country scale. On this map it is then possible to redistribute the mostly updated population estimation (available at least at national scale) in order to evaluate the exposed population.

Algorithm
The steps to be performed in order to obtain an updated population distribution map are reported hereafter: 1. Input pre-processing: Nightlights intensity and built-up area data are regridded on a reference grid and normalized between 0 and 1. Resolution of the reference grid is the same of the nightlights intensity layer, which is usually the dataset provided at the lower resolution. The choice to regrid the built-up data on the nightlights grid allows to better calibrate the relationship between them generating a single urban density value for each cell of the nightlights intensity map. In fact, by regridding, the built-up values are transformed in urban density values, measured as percentage of the pixel of the new grid occupied by built-up areas. 2. Curves calibration: Urban density curves are calibrated in one or more regions of the case study area, evaluating the relationship between normalized nightlights intensity and urban density data (obtained from previous step) available for the same time period. The choice of the number of curves to calibrate strongly depends on the distribution of the levels of urban development inside the study area. If the study area has a homogeneous urban development a unique curve can be derived. 3. Urban density update: If new nightlights intensity data are available, the corresponding map is regridded and normalized. The curve is then applied to this new nightlights intensity map to estimate updated urban density for the whole case study area. Therefore, the urban density for each pixel i of the map is evaluated as: Where UD i is the urban density for pixel i, L i is the normalized light intensity for the same pixel and f is the calibration curve which better fits the experimental data, maximizing the R 2 measure. 4. Population distribution: Updated urban density is used to redistribute new available population data provided at least at national scale, applying the following expression: Where i is the pixel in which population is redistributed; j is the area for which population data are available (e.g.: country/region/district); P i,j is the redistributed population in pixel i; P j is the total population in area j; UD i,j is the updated urban density UD i for pixel i belonging to area j. 5. Exposure assessment: Distributed population is crossed with hazard maps and population exposed is estimated.
A graphical representation of the methodological steps above described is reported in Figure 1.

Figure 1.
Graphical representation of the proposed algorithm for the estimation of population exposed to natural hazards, which reports the main procedural steps, input data, intermediate products and the final output.

Test and validation
In order to test the performance of the methodology, an application for the assessment of the number of people exposed to riverine flood in Italy is performed. Input data used for this application are described into detail in Section 2.3 and the obtained results are presented and discussed into detail in Paragraph 3 and Paragraph 4 respectively. In this test, the probabilistic relation between nightlights and urbanization density is derived by fitting a calibration curve for the reference year 2011. With the purpose of testing the goodness of the obtained relationship, the calibration curve evaluated for three specific regions (Lombardy, Tuscany and Sicily) is applied on the same satellite-based nightlights map of 2011 producing an estimated urban density map at national scale. This new urban density map is used to redistribute population data provided at regional scale by ISTAT, applying (2). Estimated distributed population is crossed with hazard maps provided by ISPRA and the estimated exposed population derived.
The choice of using the same reference year for calibration and testing, without realizing an update with more recent nightlight data, has been guided by the fact that mostly updated information on population exposed to riverine flood at municipal scale for Italy (to be used as independent validation data) are available for year 2011.
Indeed, validation is performed comparing the evaluated exposed population obtained from testing with reference exposed population. Reference exposed population is obtained redistributing population provided by ISTAT at the municipal scale on the original urban density map derived from GUF and then crossing the obtained distribution with the hazard zones layer provided by ISPRA. The validation is performed at municipal scale, which is the scale at which reference data are provided. For this reason, both reference and evaluated exposed population from testing are aggregated at the municipal scale to perform the comparison.
A graphical representation of the steps performed for testing and validation described above are reported in Figure 2. Graphical representation of the workflows performed to test and validate the proposed methodology, applied for the assessment of the number of people exposed to riverine flood in Italy.

Input data
In order to test the methodology for riverine flood exposure assessment in Italy, two remote sensing based datasets have been identified as input data for curves calibration: the DMSP night-time light series as nightlights intensity dataset, an the GUF as built-up area map. Furthermore, in order to estimate population exposed to natural hazards, regional population data provided by ISTAT from 2011 census and riverine flood hazard zones provided by ISPRA are used. Both ISTAT and ISPRA datasets are provided as vector layers.
The layer used for validation is obtained using the same riverine flood hazard zones by ISPRA, the urban density directly obtained from GUF dataset and population data provided by ISTAT from 2011 census aggregated at municipal scale. The list of all the datasets used for test and validation are reported in in Table 1. Each of them is then described into detail in the following subsections.  [26]. Data are provided by 6 different satellite platforms for a total of 34 scenes. Nightlights data are expressed as yearly averaged digital number (DN) values ranging between 0 to 63, where 0 represents complete darkness and 63 bright areas. Data are provided in raster GeoTiff format with spatial resolution of 30 arc sec (0.00833 • , nearly 1 km at the equator). The spatial extension of nightlight images is between 75 • N and 65 • S and 180 • W and 180 • E. The database contains three different data types: Average Visible, Stable Lights, and Cloud Free Coverages. Stable Lights contains the lights from cities, towns, and other sites with persistent lighting, including gas flares. Ephemeral events, such as fires have been masked. Then the background noise was identified and replaced with zero values. Areas with zero cloud-free observations are represented by the value 255. The latter was considered as the most suitable information to describe the level of exposure, and then used in the analysis, as it avoids including temporary changes in lighting. The DN derived was subsequently normalized between 0 and 1 to ease comparison with the urbanization percentage.

Global Urban Footprint
The Global Urban Footprint is a world-wide map of human settlements at very high resolution produced by analyzing a global coverage of TerraSAR-X and TanDEM-X images collected in the context of the TanDEM-X mission launched by DLR in 2007 [9]. The global map has been produced starting from around 180.000 single TerraSAR-X/TanDEM-X image products for the reference year 2011, which have been processed by DLRs operational Urban Footprint Processor (UFP), followed by an automated post-processing approach. The final product is a binary, thematic raster dataset in GeoTiff format, with reports value 255 for built-up areas and 0 for non-built-up areas, available at a geometric resolution of 0.4 arc seconds ( 12 m, near the equator) [8,28].

Italian riverine flood hazard zones
Riverine flood hazard maps have been collected by ISPRA, in accordance with the European Floods Directive (2007/60/EC) which requires each Member State to map Flood Hazard and Flood Risk for significant flood risk areas. The riverine flood hazard maps produced by ISPRA report the delimitation of the geographic areas that could be affected by floods for three probability scenarios: low probability (P1), medium probability (P2) and high probability (P3). The identification and delimitation of flood prone areas for different hydraulic hazard scenarios has been carried out by applying two different approaches: (i) an hydro-geomorphological and historical-inventory based methodology, which combines LiDAR data, aerial photos, satellite images, geological maps and historical records of areas affected by past flood events; (ii) an hydrological-hydraulic analysis, in which flood prone areas are identified determining the hydrographs for different return periods and applying a numerical Preprints (www.preprints.org) | NOT PEER-REVIEWED | Posted: 21 October 2020 doi:10.20944/preprints202010.0425.v1 hydraulic model [29]. Riverine flood hazard maps for the entire Italian territory are provide as a vector layer. The release used in this paper refers to year 2017.

Italian census of population
Population data used for testing and validation have been extracted from the data warehouse of the "15th Population and Housing Census" carried out by ISTAT [30]. The warehouse contains a wealth of information, disaggregated to the sub municipal level or aggregated at municipal, regional or national scale, on the demographic and social structure of the population usually resident in Italy and on Italy's housing stock. The census reference date is 9 October 2011. The census variable which have been used for this paper is the final demographic data on the population resident in Italy, used at regional scale for testing and at municipal scale for validation.

Results
As mentioned in Section 2.2, the performance of the proposed methodology is tested and validated in a case study in which the number of people exposed to riverine flood in Italy is assessed. The workflow applied for testing and validation is reported in Figure 2. The main results obtained at each step of the workflow are reported hereafter.

Test of the methodology for the assessment of the number of people exposed to riverine flood hazard in Italy
In this test, DMSP night-time lights and GUF for the reference year 2011 are used in three calibration regions to obtain a statistical relationship between them. The derived curve is then applied on the same DMSP map of 2011, producing an estimated urban density map at national scale to be used to redistribute population data provided at regional scale by ISTAT. Estimated distributed population is then crossed with ISPRA riverine flood hazard maps to estimate exposure. The outputs of each of these passages are described into detail in the following subsections.

Input pre-processing
The DMSP night-time light series and the GUF for the reference year 2011 are regridded on a reference grid (with the same resolution of the DMSP layer, which is around 1 km) and normalized between 0 and 1. While regridded at 1km, the GUF data are transformed in a urban density value, measured as percentage of the pixel of the new grid occupied by built-up area pixels from GUF. The obtained normalized products are shown in Figure 3. The comparison between the two maps of Figure 3 allows to observe that, despite the significant difference in the resolution, the two products show a similar pattern. More specifically, all products are able to well capture big cities and high density urban areas. Nevertheless, the DMSP night-time light tends to overestimate the light intensity around these high density urban areas because of the glare phenomenon. Other night-time light products with better resolution than DMSP, if available, may reduce the extension of this noise.

Curve calibration
The probabilistic relationship between normalized nightlights and urbanization density data is evaluated in three different calibration regions, potentially represented by different levels of urban development: Lombardy, Tuscany and Sicily. The curve that best represents the relationship between the average urbanization values for a given light intensity is a 6 degree polynomial, forced to assign 100% urbanization to 100% light intensity. Nightlights intensity is first divided into 20 classes with equal width (5%) and the distribution of urban density in each class, for each calibration region, is estimated. The fitting is then performed on the 50th quantile of the urbanization density of the three calibration regions. More specifically, the identified curve is described by the following equation: which can be applied to compute UD i according to (1). The values assumed by the coefficients c n , with n = 1, . . . , 9, are reported in Table 2. This calibration curve is able to fit very well experimental data, reaching a R 2 value up to 0.94. The visual comparison between the 50th quantiles for the three calibration regions and the fitting curve is reported in Figure 4.

Urban density estimation
The calibration curve described by (3) is applied using as input the normalized DMSP night-time light (the same layer used to calibrate the curve) producing an estimated urban density map at national scale. As well shown in Figure 5, the comparison between estimated and original urban density shows very consistent patterns. Effects of light diffusion (i.e. glare) are stronger in big urban areas, such as Milan (squared in red in Figure 5). The availability of higher resolution satellite light intensity data may potentially reduce the effect.

Population distribution
The estimated urban density map is used to redistribute at finer scale population data provided at regional scale by ISTAT, by applying (2). The resulting population distribution is reported in Figure 6. . Estimated population distribution obtained redistributing regional population data from ISTAT on the urban density map derived from DMSP, detailed for for North-West Italy (a) and distributed for whole Italy (b).

Exposure assessment
Estimated population distribution is crossed with riverine flood hazard maps provided by ISPRA in order to obtain an estimation of people exposed to flood. ISPRA hazard maps, which report the delimitation of the areas that could be affected by floods according to low (P1), medium (P2) and high Preprints (www.preprints.org) | NOT PEER-REVIEWED | Posted: 21 October 2020 doi:10.20944/preprints202010.0425.v1 (P3) probability scenarios, allows to estimate the number of inhabitants per pixel located in hazard zones P1, P2 and P3 respectively. The population distribution layer has been resampled from DMSP resolution (1 km) up to 100 m to facilitate the overlapping with the detail of the hazard zones, which are provided as a vector layer. The ISPRA riverine flood hazard maps for the three probability scenarios and the number of inhabitants per pixel located in hazard zone P1, as an example, are illustrated in Figure 7.

Validation against reference exposed population
According to the procedure described in Section 2.2, validation of the resulting exposed population is performed through a comparison with reference exposed population, obtained redistributing population provided by ISTAT at the municipal scale on the original urban density map derived from GUF and then crossing the obtained distribution with the hazard zones layer provided by ISPRA. Since reference population data are provided at municipal scale, both reference and evaluated exposed population are compared at the municipal scale. The visual result of the comparison is presented in Figure 8. The correlation between estimated and actual exposed population per municipality for whole Italy is shown in two scatter-plots (Figure 9), which evaluate the correlation for P1 (low probability) and P3 (high probability) hazard zones respectively. Figure 9. Scatter-plots that show the correlation between estimated and actual exposed population at municipal scale in Italy for P1 (a) and P3 (b) hazard zones.
Data presented in Figure 9 proof a good correlation between estimated and reference exposed population at municipal scale, reaching a value of R 2 equal to 0.88 for exposed people in P1 (low probability) hazard zones and equal to 0.85 for exposed people in P3 (high probability) hazard zones.

Discussion
The test and validation of the methodology has allowed to evaluate its goodness for future applications, in which it would be possible to use updated data about nightlights intensity and population estimation at national or regional scale in order to quickly perform an update of population exposed to hazards. Indeed, the validation with independent data shows good agreement between estimated and reference population exposed to riverine flood hazard (Figure 8), demonstrating a strong consistency in the statistical relationship expressed by (3). Furthermore, the shape of the obtained correlation curve (Figure 4) shows an agreement with similar studies presented in literature, such as the correlation curve between population density and nightlights intensity presented by Bagan and Yamagata [20] for Japan. In their study, the correlation becomes significant for values of DMSP nightlights intensity above 61-62 (97-98 % of the maximum intensity). This result is absolutely coherent with the choice to force the calibration curve described by (3) to drawn through (1,1), in order to ensure the maximum urbanization density in correspondence to the maximum nightlights intensity.
Scatter-plots reported in Figure 9 proof a very strong correlation between estimated and reference exposed population at municipal scale, specifically for exposure evaluated in P1 (low probability) hazard zones. In this case, correlation reaches a value of R 2 around 0.88. Nevertheless, for P3 (high probability) hazard zones the correlation slightly decreases (R 2 around 0.86) due to the increasing inconsistencies between nightlights data resolution used as main input to update population distribution maps and the detail of the hazard zones. Indeed, higher probability flooding scenarios are characterized by a smaller and more scattered extent such as P3 hazard zones, which are very thin buffers around small rivers. Therefore, the proposed methodology provides better results on large flood hazard areas, which are more consistent with the light intensity map resolution. The discrepancy between results obtained with P1 and P3 hazard zones is not so considerable at national scale (the difference between the two R 2 is very small) because of compensation. Nevertheless, at finer scale this discrepancy is significantly highlighted. In Figure 10, the correlation between estimated and actual exposed population per municipality for Lombardy region is shown in two scatter-plots, which evaluate the correlation for P1 (low probability) and P3 (high probability) hazard zones, respectively. Preprints (www.preprints.org) | NOT PEER-REVIEWED | Posted: 21 October 2020 doi:10.20944/preprints202010.0425.v1 Figure 10. Scatter-plots that show the correlation between estimated and actual exposed population at municipal scale in Lombardy region for P1 (a) and P3 (b) hazard zones.
In this case, the correlation reaches a value of R 2 around 0.82 for P1 (low probability) and significantly decreases for P3 (high probability) hazard zones (R 2 around 0.50) due to the increasing inconsistencies between nightlights data resolution and the detail of the P3 hazard zones. These inconsistencies can be well observed in Figure 11, where derived and reference population distributions for the city of Milan are overlapped by P3 hazard zones. More specifically, it is possible to observe that the hazard zones are very narrow along the river network, having a cross dimension significantly lower than the resolution of the derived population distribution, obtained from DMSP nightlight data. If hazard maps with a lower detail would be used, such as seismic hazard maps, this issue would be potentially be avoided.
The goodness of the relationship described by (3) can be appreciated also comparing the estimated and the original urban density maps, both reported in Figure 5. The comparison between the two maps shows very consistent patterns. Nevertheless, light diffusion is stronger in big urban areas, such as Milan (squared in red in Figure 5). Indeed, the DMSP night-time light used as input to produce the estimated urban density map, tends to overestimate the light intensity around big urban areas because of the glare phenomenon. The availability of higher resolution satellite light intensity data may potentially reduce the glare effect. Furthermore, products with better resolution than DMSP may also increase the correlation between nightlights intensity and urban density, improving the overall performance of the proposed methodology.

Conclusions
The up-to-date assessment of the built-up area and the number of people exposed to natural hazards represents a key component to perform a reliable risk assessment and to implement successful risk reduction strategies. The work presented in this paper has allowed to identify, test and validate a general and simple remote sensing based procedure for the quick update of the population exposed to natural hazards in rural and urban areas. The methodology is based on the identification of a relationship between nightlights intensity and urbanization density, which can be assessed in some calibration regions and then applied to extrapolate the urbanization for the whole country at different time steps, when only nightlights are available.
The validation with independent data shows good agreement between estimated and reference population exposed to riverine flood hazard, making this approach suitable for future applications. Nevertheless, there are some factors which can influence the reliability of the results, delimiting the field of application of this approach. A first factor is the impact of the resolution of nightlights intensity data on the estimation of population distribution around big urban areas because of the glare phenomenon. A second limitation is represented by the fact that the methodology provides better results on large hazard areas (e.g. flood hazard area P1), which extent is sufficiently larger than the low resolution of the light intensity map, making unnecessary a more detailed spatial distribution of the population. Nightlights products with better resolution than DMSP data may overcome all these two limitations, increasing significantly the performance of the methodology.
Despite the current limitations, the presented approach performs well at the municipal and supra-municipal scale, producing updated exposure maps which can be used to evaluate aggregated risk and direct damage estimations. The possibility to easily and quickly update exposure maps with an higher refresh rate compared to available exposure cartography, makes this approach particularly suitable for applications in developing countries, where exposure may change at sub-yearly scale.

Conflicts of Interest:
The authors declare no conflict of interest.

Abbreviations
The following abbreviations are used in this manuscript: