Natural Disaster Damage Indices Based on Remotely Sensed Data: An Application to Indonesia

Combining nightlight data as a proxy for economic activity with remote sensing data typically used for natural hazard modeling, this paper constructs novel damage indices at the district level for Indonesia, for different disaster events such as floods, earthquakes, volcanic eruptions and the 2004 Christmas Tsunami. Ex ante, prior to the incidence of a disaster, district-level damage indices could be used to determine the size of the annual fiscal transfers from the central government to the subnational governments. Ex post, or after the incidence of a natural disaster, damage indices are useful for quickly assessing and estimating the damages caused and are especially useful for central and local governments, emergency services, and aid workers so that they can respond efficiently and deploy resources where they are most needed.


Introduction
Quickly assessing and estimating the damage caused after the incidence of a natural disaster is important for both central and local governments, emergency services and aid workers, so that they can respond efficiently and deploy resources where they are most needed. Recently, remote sensing technologies have been used to analyze the impact of disasters, such as hurricanes (Myint et al., 2008;Klemas, 2009), floods (Haq et al., 2012;Wu et al., 2012Wu et al., , 2014Chung et al., 2015), landslides (Nichol et al., 2006), earthquakes (Fu et al., 2005;Yamazaki & Matsuoka, 2007), wildfires (Holden et al., 2005;Roy et al., 2006), volcanoes (Carn et al., 2009;Ferguson et al., 2010) and tsunamis (Römer et al., 2012). These remote sensing techniques are useful for providing quick damage estimates shortly after the disasters giving emergency services a chance to respond quickly and local governments an overview of estimated costs and necessary repairs.
In addition to their usefulness in the aftermath of a disaster, estimates of the potential damage associated with a natural disaster are also useful for policy making prior to the realization of the natural hazard event. In many cases the incidence of a natural hazard event can turn into a natural disaster simply because of inadequate preparation ex-ante. Indonesia, for example, is highly exposed to natural disasters by being situated in one of the worlds most active disaster hot spots, where several types of disasters such as earthquakes, tsunamis, volcanic eruptions, floods, landslides, droughts and forest fires frequently occur. The average annual cost of natural disasters, over the last 10 years, is estimated at 0.3 percent of Indonesian GDP, although the economic impact of such disasters is generally much higher at local or subnational levels (The Global Facility for Disaster Reduction and Recovery, 2011). The high frequency of disasters experienced has important impacts on expenditures by local governments that could be anticipated, at least in part, through upward adjustments in the annual fiscal transfers from the central government to the subnational governments. 1 Such ex-ante adjustments in the level of fiscal transfers would be more useful if they could be based on estimates of the potential damages associated with the incidence of a natural disaster as opposed to estimates of the intensity of the potential natural hazard that might occur. However, although in recent years there has been much progress towards the modeling of the main natural hazards, there continues to be a scarcity of estimates of the damages associated with the incidence of these disasters. The value of damage caused by a natural disaster is typically a complicated function of the size of population living in that area, the level and type of economic activity carried out, the value of the physical infrastructure in place, and the resilience of infrastructure and people's livelihoods to the natural hazards.
This paper fills some of the gaps in the literature by using different remote sensing sources and data on the physical characteristics of the events to construct four damage indices for natural disasters in Indonesia. The indices cover floods, earthquakes, volcanic eruptions and a tsunami, and are all weighted by local economic activity in an area, and then aggregated up to a district level. 2 All data used in the construction of the indices are free and publicly available, making the methods used a potentially very useful alternative for both central and local governments to quickly get a rough estimate of the damages caused by a disaster (either ex-ante or ex-post). 3 Importantly, all of the indices constructed take into account local exposure. Given limited access to highly disaggregated local economic activity data, nightlight intensity derived from satellite imagery has proved to be a good proxy; see, for instance, Henderson et al. (2012), Hodler & Raschky (2014) and Michalopoulos & Papaioannou (2014). By utilizing the grid cells of approximately 1 square kilometer we can break down areas in cities and districts into where they are busiest, and thus take into account not only the local physical characteristics of a natural disaster but also the local economic activity exposed to it.
The paper is structured as follows. Section 2 of the paper discusses in more detail the incidence and types of natural disasters. Section 3 discusses the nightlights data. Sections 4-7 discuss in detail the construction of the four damage indices, while section 8 concludes.

Natural Disasters in Indonesia
Natural disasters are prevalent events across most parts of Indonesia. According to the Indonesian National Disaster Management Authority (BNPB) there were more than 19,000 natural disasters in the period (National Disaster Management Agency, BNPB, 2016, making Indonesia a useful country for any natural disaster analysis. The most frequent disasters are floods and landslides (52 percent), strong winds (21 percent) and fires (15 percent), while the most damaging ones are earthquakes, tsunamis and volcanic eruptions, which all cause major damage to buildings and infrastructure in addition to the human casualties. The deadliest year according to the BNPB data was 2004, where there were more than 167,000 deaths due to natural disasters and 166,671 of them stemming from the tsunami in December 2004.

Floods
The tropical climate of Indonesia often leads to annual floods. The BNPB data registered more than 10,000 incidents of floods or landslides leading to more than 3,500 Source: G.R. Brakenridge (2016) 4 Magnitude is defined as: M = log(D * S * AA), where D is the duration of the flood; S is the severity on a scale consisting of 1 (large event), 1.5 (very large event) and 2 (extreme event); and AA is the size of the affected area. Flood events registered by DFO have mainly been derived from news and governmental sources. 5 The provinces where no large scale centroid was present were Bangka Belitung, Riau Islands, Kalimantan Barat, Yogyakarta, Sulawesi Barat, Kalimantan Utara and Maluku. Note that some of these, like Kalimantan Utara, Kalimantan Barat, Sulawesi Barat and Yogakarta, did most likely experience large scale flood during these years, but that the centroid was in another province. The remaining three provinces consist mainly of smaller islands, so the flooded area will most likely not constitute a large scale flood event.

Earthquakes
Due to Indonesia's location inside the Pacific Ring of Fire, one of the most seismically active areas in the world, it is often struck by earthquakes. BNPB counted almost 400 earthquakes from 2001 to 2015, with the largest number of casualties coming from the tsunami created by a 9.0 earthquake located off the coast of Aceh, otherwise known as the earthquake that caused the 2004 tsunami. Apart from that, there were more than 8,000 registered fatalities due to earthquakes over the same period. Overall, this makes earthquakes the deadliest of the natural disasters that strike Indonesia. Figure 2 shows how common earthquakes are in Indonesia by displaying contour maps 6 of all earthquakes of magnitude 5.0 and above that struck Indonesia from 2004 through 2014. In total, the United States Geological Survey (USGS) registered 261 earthquakes. 7

Volcanic activity
Indonesia has the highest number of active volcanoes in the world, numbering almost 150. Of these, many have had eruptions in both more historical times and after the year 2000. The most famous eruption is probably the explosion of Krakatau in August 1883, when two-thirds of the Krakatau Island erupted and disappeared, killing more than 35,000 people and causing a global mini ice age and weather disruptions for years. BNPB have registered 92 eruptions over our 15-year time period and more than 60 major volcanoes that have had eruptions since 1900. The most recent one is the 2010 Mount Merapi eruption 6 These maps are also known as ShakeMaps, which are produced by USGS. 7 There are 1,002 earthquakes registered by USGS that were of magnitude 5.0 or more that had a point with a PGA of at least 0.05 within Indonesia. Many of these points create little to no damage. The 261 earthquakes mentioned above are quakes that are mostly contained within Indonesia. that killed 324 people and dislocated more than 320,000.
In addition to the BNPB data, during the years 2004 through 2015 the Darwin Volcanic Ash Advisory Centre (DVAAC) had 587 days where they issued a red warning, implying an ongoing or imminent volcanic eruption. The most active volcanoes -measured by number of days with red warnings -are shown in Table 1, with the top 5 volcanoes constituting almost 75 percent of the red warnings. These are also volcanoes that have been in the media, with Merapi already mentioned and Sinabung, which had several eruptions in 2010, 2013 and 2014. These two volcanoes are located close to densely populated areas, with Sinabung located in North Sumatra and Merapi in central Java.

2004 Christmas Tsunami
The Christmas Tsunami in 2004 is the worst singular natural disaster during the modeling period, and one of the worst natural disasters in world history. As seen in the photo in Figure 3 the destruction was absolute in parts of Indonesia. The total death toll across Indonesia and 13 other countries was more than 230,000 people and there were many more missing. In addition, the World Bank (2005) estimated a total economic impact of 4.5 billion US dollars. The official BNPB data for Indonesia estimates 166,671 deaths due to the tsunami.
The cause of the tsunami was an earthquake of magnitude 9.0 150 miles south-south east of Banda Aceh on the morning of 26 December. This quake created a tsunami with waves more than 20 meters high at the highest. Due to the fault line of the earthquake being in a north-south direction, the greatest strength of the tsunami was in an east-west direction (Athukorala & Resosudarmo, 2005). This led to the largest damages being in the northern part of Sumatra, in the province of Aceh, where entire villages were wiped out as seen in the photo of Banda Aceh (Figure 3).

Nightlight Data
Natural disasters are inherent local phenomena in that they either affect only parts of areas and/or affect parts within areas differently. It is thus important to take the local population/asset exposure into account when constructing more aggregate proxies. Arguably one would like to have measures of exposure as spatially disaggregated as possible. For a country like Indonesia, data are usually sparse and at a very aggregated spatial level.
An alternative approach is thus to use nightlights as a proxy for local economic activity. As a matter of fact, nightlights have found widespread use where no other measures are available; see, for instance, Henderson et al. (2012), Hodler & Raschky (2014) and Michalopoulos & Papaioannou (2014). In Henderson et al. (2012), Indonesia is used as an example of using nightlights to capture an economic downturn following the Asian financial crisis in the late 1990s. Their results show that swings in GDP change can generally be captured. Nevertheless one has to account for factors such as cultural differences in light usage, latitude and gas flares. In our case this is unlikely to affect our results since we use nightlights to capture exposure within a country rather than across countries.
The nightlight imagery we employ is provided by the Defense Meteorological Satellite Program (DMSP) satellites. In terms of coverage each DMSP satellite has a 101 minute near-polar orbit at an altitude of about 800km above the surface of the earth, providing global coverage twice per day, at the same local time each day, with a spatial resolution of about 1km near the equator. The resulting images provide the percentage of nightlight occurrences for each pixel per year normalized across satellites to a scale ranging from 0 (no light) to 63 (maximum light). Yearly values were then constructed as simple averages across daily values of grids, and are available from 1992. 8 We use the stable, cloud-free series; see Elvidge et al. (1997).
The data revealed 414,644 cells which had a nightlight value greater than 0 in them at least once during the period 2001-2013. Figure 4 -containing all cells with nightlights in 2012 -shows that the large cities and densely populated areas on Java, Sunda Islands, coastal Kalimantan and Sumatra are fully covered in lights. Inner parts of Kalimantan and most parts of New Guinea are more sparsely lit. The modeling of floods can be done by remote sensing (Brakenridge & Anderson, 2006;Wu et al., 2012;Haq et al., 2012) or through a combination of weather data and GIS systems as for example in Knebl et al. (2005); Asante et al. (2007); Dessu et al. (2016). We utilize the latter, as remote sensing is useful for assessing whether an area is flooded or not, but it is weaker on modeling the intensity of the flood. Moreover, cloud cover generally limits the accurate detection of floods from remote sensing sources.
To model floods we have decided to use the Geospatial Stream Flow Model (GeoSFM) which is a software that is "designed to use remotely sensed meteorological data in data sparse parts of the world" (Artan et al., 2008). GeoSFM was developed by USGS and USAID and is a hydrological modeling tool used to model stream flows across large areas, in particular areas where highly localized data are lacking. It has been used in regions such as the Great Horn of Africa (Asante et al., 2007;Mati et al., 2008;Dessu et al., 2016) and Nepal (Shrestha et al., 2011), with Dessu et al. (2016 finding that the model captures 76% of the monthly average variability, making it useful for flood simulation. The inputs needed to model stream flow for basins are soil-and terrain-based -such as digital elevations models (DEM) and land cover and soil data -and weather-based, such as precipitation and potential evapotranspiration (PET) data. The HYDRO1K data set from USGS, which is a DEM made for hydrological modeling based on the USGS' 30 arc-second DEM of the world, is used as elevation input. The land cover data are the Global Land Cover Characterization (GLCC) data set also from USGS, while the soil data are from the FAO Digital Soil Map of the World.
The daily precipitation data are from the 3-hourly data set from the Tropical Rainfall Measurement Mission Project (TRMM) and the PET data are 6-hourly data from the Global Data Assimilation System (GDAS), both data sets are aggregated up to daily data. The PET data are available from February 2001 and onwards, while the precipitation data are available for the period 1998-2014. Given that we only have nightlight data through 2013, we will focus on floods for the period 2001-2014.
GeoSFM uses the inputs to construct basins based on the terrain and then uses a linear soil moisture accounting routine to model surface runoff and soil moisture based on precipitation and PET. It is worth noting that although a more complex and better non-linear routine is also supported, it does not work well for our more generalized macro-modeling with fairly low resolution data. Finally, GeoSFM models the stream flow for each basin for each day of our time period.
Note that GeoSFM does not model coastal floods, nor does it model flash floods in areas where there are no rivers or streams of a certain length. Figure 5 shows that there are parts of Indonesia and even one province -Riau Islands -which have no basins. Another weakness is that it does not take into consideration the specific terrain within each basin. Floods are generally very localized events and the low resolution of our data makes it impossible to model the intensity of the stream flow within a basin and also causes some river outlets to be slightly inland instead of running all the way to the ocean.

Creation of Index and Results
The first part of constructing the index involves defining when a flood event is happening. In Wu et al. (2012) they propose four runoff based methods to define a flood threshold, and in addition Wu et al. (2014) propose a slightly modified flood threshold definition with a point being flooded when: where R is the routed runoff in millimeters, P 95 is the 95th percentile value and σ is the standard deviation of the routed runoff. Q is the discharge in cubic meters.
We found that with the GeoSFM modeled data, runoff was not a good proxy for flooding, due to it only capturing a limited number of floods. Discharge, Q, was a better proxy, leading to a new -but very similar -equation: Q > P 95 + σ and Q > 10m 3 /s (2)  Kreibich et al. (2009) look at different parameters such as velocity, depth, energy head, stream flow and intensity. They find velocity to be a poor parameter for assessing damage, while water depth and energy head show the best results. Stream flow and intensity are also weak as parameters. Finally, Merz et al. (2010) assess different damage influencing parameters and point to the fact that most "damage influencing factors are neglected in damage modeling, since they are very heterogeneous in space and time, difficult to predict, and there is limited information on their (quantitative) effects". Overall, there is limited support in the literature for a strong correlation between these parameters and damages on anything but a very localized scale.
As for assessing the damage itself, Merz et al. (2010) discuss damage functions and the two main approaches, which involve one empirical approach where damage data are collected after the flood and one synthetic approach where they construct potential what if-scenarios. Once again the assessments rest on very localized data, which we do not have for Indonesia. Overall, it means that we cannot expect anything more than rough estimates. A common denominator for the papers mentioned above is that there is some measurement of intensity. Given that stream flow is an intensity proxy, we have used that to construct a simple measurement for intensity. The equation is: where I b,t is the intensity of the flood in basin b at date t, Q b,t is the stream flow in the same basin at the same time andQ b and σ b are mean and standard deviation of stream flow in b. The intensity is set to zero if the flood threshold -95th percentile plus 1 standard deviation above the average -has not been exceeded. By normalizing, we obtain a measure that is comparable across all regions and that is independent of the absolute river flows. The assumption is that people living close to rivers will be prepared for variations in water levels, and that people living close to rivers with highly variable stream flows are more prepared for these events than people living close to more stable rivers.
To aggregate the flood impact each basin is weighted based on the nightlights in it. The weights per basin, W b,t−1 , used are: where L b,t−1 is the sum of lights in basin b one year, t − 1, before the flood year and L p,t−1 is the same at a province level.
Finally, the weights from (4) are multiplied with the intensity from (3) to get the overall flood impact, FI b,t in that basin on the province: One thing to note here is that for basins that span several provinces or districts, we have assumed the same intensity, but the weight is based on nightlights within each individual province.

Results
The stream flow was simulated for 5,082 consecutive days, from 1 February 2001 to 31 December 2014. 9 Table 2 shows that the top 10 basins with most flood days had close to 200 days of flooding over the 14-year period. As expected, these basins do overlap with some of the busiest flood areas according to the DFO, as shown in Figure 6. The lower part of Table 2 reveals that the driest basin had a mere 12 days of flooding. All 33 provinces with a basin had days that went above our flood threshold set in Equation (2).
All months in our model have flood events, but there are big differences. The range goes from 527 events every March and down to 154 events every August, with the traditional rainy season (November-March) producing the highest number of flood events, whereas the dry season months (June-October) are the driest. Aggregating the numbers for the rainy season, there are 2,215 events every year across the basins, while there are only 995 events every year during the dry season.
Total number of basins that are partly or fully inside Indonesia is 495, and these basins have a total of 55,605 flood events or slightly more than 112 per basin. In other words, the average basin has been flooded for a total of 8 days a year over the 14-year period in question. This is not entirely unexpected given the climate in Indonesia and the way our threshold is made. Also, if we compare with the DFO data where they have 3,808 floods of magnitude 4 or higher through their period from 1985-2016, which converts to almost 123 fairly large scale flood events per year, our model provides a reasonable proxy for events.
Even though the results seem logical on a per basin basis, the time steps in the model are 1 day at a time, which is too slow for the unfolding of a flood event, implying that downstream basins that would normally fill up very quickly will now only be filled up the day after, and then the next basin will be filled two days later and so forth. This means that the amount of days with floods are inflated. We believe that this does not affect our results much, though, as the number of events per province will not affect the end results, since we weigh by affected nightlight and not by number of days of floods.
Despite the numerous floods in Indonesia, they generally do not affect a large percentage of the population, as per Table 3. The mean of nightlights when excluding areas with 0 nightlight is 3.39 percent.
9 For Bali we did it for 5,080 days due to problems with 30 and 31 December 2014.   The DFO flood database is mostly based on news sources, providing an overview of the big floods in Indonesia. To check the database against the GeoSFM model, the focus will be on the largest events of magnitude 6 and above. Given how the DFO data do not give any intensity estimates and focus primarily on displacement numbers and area, while our model is driven by intensity the comparison will only focus on whether GeoSFM results do overlap in time and/or province with the centroid of the DFO floods. Table 4 shows the DFO data on the left side, first column being the start month of the flood, followed by duration, magnitude, the province where the centroid of the flood is, dead and displaced. The right side shows the modeling results where the focus is on duration. The first column under model results shows the number of days for the centroid province, then overall number of days with floods anywhere in Indonesia during the period, then looking at the same island -using that as a proxy for neighboring regions -where one examines total days the island provinces were flooded during the flood and finally how many of the days of the flood duration that a province on the same island was flooded. Generally, the model performs well, in particular on Sumatra and Kalimantan (Borneo), with the example where the 2008 flood was captured for all 25 days in the centroid province. Overall, it shows at least one flooded basin on Kalimantan and Sumatra for 85% of the days the major floods happened. The results are somewhat worse on Java, where only 37% of the days have a flood. A primary reason for this might be that Java is very narrow and hence the streams are short and might not be captured in our model. Java also has larger percentage of land not covered by a basin, ref Figure 5, also due to its narrowness which makes the low resolution landcover data underestimate the size of the basins.

Aggregated numbers
Finally, to aggregate up to a district or province level, we have used a simple method for the total damage experienced per year: where j is the province or district, T is the year, sum of t is all the days for year T , sum of b are all the basins in the province or district and FI b,t is the flood impact from Equation 5. Normally one might use an average flood impact across the year, but by doing this, we capture repeated flood events and areas that experience generally high flooding.
Using the above method, Table 5 provides the aggregated data for all provinces across all years. The impact is fairly even for the most impacted ones, with the impact numbers for the top 10 ranging from 44 to 55. Furthermore, Sumatera Selatan, Lampung and Yogyakarta make up 8 of the top 10 impacted provinces. The overall picture fits with the DFO floods in Figure 1, with the populous provinces in Java, Sumatra and Sulawesi being impacted, whereas the smaller island provinces and parts of Kalimantan are not affected much. For some of the island provinces the numbers are probably underestimated, no basins will have been constructed and modeled there due to the many small islands.
Finally, Table 6 shows the most impacted districts over the years 2001 through 2014. The impact is much larger than for the provinces as one would expect due to the more localized data and impact. The districts are also more geographically spread out than the provinces.    This paper uses the latter method, by utilizing ShakeMaps from USGS, which are automatically generated maps providing several key parameters following an earthquake, such as peak ground acceleration (PGA), peak ground velocity (PGV) and modified Mercalli intensity (MMI). More specifically, the ShakeMaps use data from seismic stations that is interpolated using an algorithm which is similar to kriging. To model the intensity in a given coordinate, the model also takes into account ground conditions and the depth of earthquake. Wald et al. (2005) point to the magnitude and epicenter locationwhich are parameters common for the entire earthquake -that have historically been used to determine how severe earthquakes were, but that the damage pattern is not just dependent on those two parameters, but also on other, more localized parameters that the ShakeMaps use to generate intensity measures. This is exemplified by several earthquakes such as magnitude 6.7 and 6.9 earthquakes in California in 1994 and 1989, respectively, where some areas further away from the epicenters got more damaged than closer areas. The reason why the more localized ShakeMaps with their ground shaking parameters are a better gauge than magnitude and epicenter distance is explained on page 13 of Wald et al. (2005) which states that: "..., although an earthquake has one magnitude and one epicenter, it produces a range of ground shaking levels at sites throughout the region depending on distance from the earthquake, the rock and soil conditions at sites, and variations in the propagation of seismic waves from the earthquake due to complexities in the structure of the Earth's crust." The ShakeMaps are interpolated grids with point coordinates spaced approximately 1.5 kilometers apart (0.0167 degrees). Figure 2 shows contoured maps of these points.
The PGA is a measure of the maximum horizontal ground acceleration as a percentage of gravity, PGV is the maximum horizontal ground speed in centimeters per second and MMI is the perceived intensity of the earthquake, a subjective measure. Figure 7 -which is originally found in Wald et al. (1999) -explains the relationship between the different parameters and the potential damage from different values. The assumption is that damage starts at an MMI level of V and a PGA of 3.9 percent of g. These levels are found for California in Wald et al. (1999), but the relationship has been found for other areas in the US in Atkinson & Kaka (2006) and Atkinson & Kaka (2007) and for places such as Costa Rica (Linkimer, 2007) and Japan, Southern Europe and Western US (Murphy & O'Brien, 1977). It should be noted that the numerical relationship differs from region to region. There are no known papers estimating these values specifically for Indonesia.

Figure 7: ShakeMap Instrumental Intensity Scale Legend
Source: Wald et al. (1999) The different measures are largely interchangeable, and in GeoHazards International and United Nations Centre for Regional Development (2001) report, they use PGA to measure damage, pointing to the fact that PGA, unlike MMI is an objective measure, implying that MMI is not easy to obtain reliably across the globe. Also, for large scale modeling, where it is unfeasible for one to model local conditions precisely, PGA serves as a good proxy for intensity of earthquakes.

Damage Index
To construct the damage index, two types of data will be used; the intensity data -expressed as PGAand building inventory data, to assess what damage one could expect for different intensities.
To take into account the building types in Indonesia, we use information from the USGS building inventory for earthquake assessment, which provides estimates of the proportions of building types observed by country; see Jaiswal & Wald (2008). The data provide the share of 99 different building types within a country separately for urban and rural areas. For Indonesia the building type information was compiled from a World Housing Encyclopedia survey, while the split between urban and rural is from the urban extent map of Center for International Earth Science Information Network -CIESIN -Columbia University et al. (2011). Without any other information available, we use this as an indication of the distribution of building types in Indonesia, but, necessarily, assume that the distribution is homogenous within urban and rural areas.
Fragility curves by building type are derived from the curves constructed by Global Earthquake Safety Initiative project; see GeoHazards International and United Nations Centre for Regional Development (2001). More specifically, buildings are first divided into 9 different types. 10 Each building type itself is then rated according to the quality of the design, the quality of construction, and the quality of materials. Total quality is measured on a scale of zero to seven, depending on the total accumulated points from all three categories. According to the type of building and the total points acquired through the quality classification, each building is then assigned one of nine vulnerability curves which provides estimates of the percentage of building damage for a set of 28 peak ground acceleration intervals.
In order to use these vulnerability curves for Indonesia we first allocated each of the 99 building types given in the USGS building inventory to one of the 9 more aggregate categories of the GESI building classification. However, to assign a building type its quality-specific vulnerability curve we would further need to determine its quality in terms of design, construction, and materials, an aspect for which we unfortunately have no further information. We instead assume that building quality is homogenous across building type in Indonesia and experiment with seven different sets of vulnerability curves, each set under a different quality ratings scenarios (ranging from 0 to 7).  To model estimated damage due to a particular earthquake event the data from the ShakeMaps and GESI are used. Then, one identifies the value of peak ground acceleration that each nightlight cell in Indonesia experiences by matching each earthquake point with its nearest nightlight cell. If the cell is further away than 1.5 kilometers or if it experiences shaking (PGA) of less than 0.05 the value is set to 0. In order to derive at a region j specific earthquake damage index, ED, the following equation is applied: where DR is the damage ratio according to the peak ground acceleration, pga, and the urban-rural qualification of cell i, defined for a set of 8 different building quality q categories. The weight w i is the same as before; the sum of the nightlights of the affected cells over the sum of the total provincial nightlights.

Results
With the above method, we find that 27 of the 34 provinces were damaged by earthquakes at some point in time. 11 Table 7 shows that the big islands Java and Sumatra have the most affected nightlight cells, given how densely populated they are and how much seismic activity is experienced there this is expected. Finally, there were 5,251 cases where the instance hit a nightlight cell that was lit. Table 8 shows that the individual nightlight cell weights are small, as expected, but the impact of having buildings of quality 4 is that within a cell that is hit, on average a bit more than 6 percent of the buildings are destroyed. 12 12 We did the same for buildings of quality 0 (the best) and 7 (the worst), something which led to maximum damage values of 33% for the best buildings and 84% for the worst versus 55% for our base case, showing how the overall damage is highly dependent on building quality information.

Aggregated data
When aggregating, a similar method as in section 4 is used, but now the aggregation is done directly by nightlight cells instead of by basin. The equation is: where t is year, l is all nightlight cells in the province or district j and ED is the damage from equation 7. Table 9 provides the full overview of damage by province, showing the large differences between the provinces and how they vary from year to year, as one can expect with highly randomized events such as earthquakes. Using Yogyakarta -which was only impacted by earthquakes in 2006 -as an example, there was a loss of 0.4 percent of the total building mass, causing damages estimated to be approximately 3.1 billion US dollars in addition to more than 5,000 deaths and tens of thousands of injured and displaced people. Apart from that, provinces on Sumatra make up 6 of the top 10 most damaging years. Even though Indonesia is often hit by severe earthquakes, even in the worst of years they only destroy about 1 percent of the buildings in a province. That being said, 1 percent of total building mass and infrastructure being damaged does constitute a significant portion of local budgets. As another example, the September 2009 earthquake in West Sumatra inflicted damages for an estimated 2.3 billion US dollars, with repair costs and losses of 64 million US dollars on government buildings (Raschky, 2013).
The numbers per district are shown in Table 10, and they suffer much more damage than the provinces, with the most impacted district losing 5 percent of building mass.

Volcano Damage Index
A volcanic eruption consists of ash clouds, pyroclastic flows and lava flows, the latter two which are very difficult to model without extensive local data. Unfortunately, there is little to no academic research that has looked into large scale volcanic modeling for all aspects of eruptions. For modeling ash clouds, Joyce et al. (2009) points to remote sensing through satellite images that detect SO 2 emissions as a potential method.
To construct a damage index for eruptions, we use a two-fold process. First, volcanic ash advisory data are used from Volcanic Ash Advisory Centers (VAAC) to detect eruptions; second, satellite images containing sulphur dioxide data from the OMI/AURA satellite are used to model the intensity of the eruptions. The OMI/AURA images have been utilized by Carn et al. (2009) andFerguson et al. (2010) to model eruption intensity.

Volcano Modeling
There are at least two types of software that are used to model eruptions. Voris (Felpeto et al., 2007), which models ash clouds, lava flows and energy cones, and HYSPLIT from Air Resources Laboratory, which models air pollution dispersion. Voris relies on highly localized data due to the lava flow and energy cone modeling, which one will not have in most cases and that does not fit well for large scale modeling across time and space. HYSPLIT does have batch inputs, but still requires several inputs per eruption, some which are not easily obtainable.
The third solution, which is related to the ash clouds mentioned by Joyce et al. (2009), is based off ash advisory data to determine whether an eruption happened and OMI/AURA satellite data to determine the scale of the eruption. This will not help with modeling lava flows and energy cones, but due to the very localized nature of the flows, there are no good sources or methods to model it for several eruptions from different volcanoes, leaving the ash clouds as a good proxy for all damages.

Ash advisory data
Ash advisories from the Darwin VAAC (DVAAC), which are ash cloud warnings for airplanes, are used to determine whether an eruption is happening. The warnings show relevant data such as volcano name, position, summit height, height of clouds and a color code that reflects the condition of the air/volcano. There are 4 different codes ranking from the normal state, green, to imminent danger of or ongoing volcanic eruption, red. Over the period from 2004 until 2015, the DVAAC issued 12,962 warnings and of these more than 90 percent were either of code red or orange. Data on advisories from code orange or below were not used, due to eruptions of this scale not being large enough to be properly captured by the SO 2 -data. By limiting the data to code red events, there are 1,785 events spread across 587 dates. 13

OMI/AURA Satellite images
To measure the intensity of an eruption, data from the Sulphur Dioxide images of the OMI/Aura project (Krotkov & Li, 2006) are used. These consist of satellite images from October 2004 and onwards. The data have been used to model ash cloud intensity and movement in several articles such as Carn et al. (2009) andFerguson et al. (2010).
The satellite images have a spatial resolution of 13 * 24km and are taken from 80km above ground. The spectral imaging shows the SO 2 vertical column density in Dobson Units and there are 14 or 15 orbits per day, where one orbit covers an area approximately 2,600km wide. A dobson unit is a measure of density, and at sea level the typical concentration in clean air is less than 0.2DU. The images contain 4 values for column density based on the center of mass altitude (CMA), which is a measure of the altitude one assumes the center of the distribution is at. There are 4 different estimates for the vertical column density, ranging from 0.9km above ground to 18km above ground.
For volcanic activities one normally uses a CMA of either 8km or 18km (OMI team, 2012), where the former is a middle tropospheric column (TRM) and is for use in medium eruptions, while the latter is an upper tropospheric and stratospheric column (STL) and is for explosive eruptions. Despite this difference, the data are interchangeable in the sense that one can interpolate from one CMA to the others.
OMI is more sensitive above clouds, which both measures mentioned will normally be. The standard deviation for both measures is as low as 0.1DU over Indonesia. The data for both STL and TRM are very similar and this paper uses the STL-data as that are most useful for the biggest events.

Damage Index
When constructing a damage index based on SO 2 values from ash clouds, one has to set thresholds for distance from the event and from the centroid of the nightlight cell and also a lower sulphur dioxide-value. There are no papers or literature that have estimated any parameter values and there are no usable local data, so the thresholds have been set somewhat arbitrarily.
First off, one wants to set a distance threshold estimating how far away the eruption could cause damage. Note that one wants ground results and not for the aviation industry, since planes can be affected very far away as evidenced by the total stand still of planes across Europe during the 2010 Eyjafjallajökull eruption. We decided to set a very relaxed condition with any point closer than 10 degrees of latitude and longitude included. Figure 9 portrays the plume approximately 7 hours after one of the biggest Merapi eruptions on 4 November 2010, where the plume moved relatively slowly and after 1000km it dissipated at the lower altitudes, which shows that our 10 degree threshold works well.
Secondly, to match the nightlight data with the OMI/Aura data, a maximum distance between a nightlight point and the nearest SO 2 point is set at 50km. The SO 2 points are fairly scattered due to cloud covers, hence to get a more consistent grid of nightlight and SO 2 values we have chosen a distance that is approximately two times the longest side of an OMI cell.
Finally, a minimum SO 2 value in Dobson Units is chosen. According to the Belgian Institute for Space Aeronomy, a typical normal level in air is 0.1DU and a strong eruption is above 10, which is the threshold value chosen.
Once the thresholds have been set, the same nightlight weighting method as for our other indices is applied and then the weights are multiplied with the SO 2 value to get an intensity value. The equation is: where i is the nightlight cell on date t, and w i,T −1 is the previously used weight where i is the nightlight cell, T − 1 is the nightlight strength from the prior year and it is divided by the sum of total nightlights in the province or district.

Results
Applying our constraints, the 587 dates with a red warning have been reduced to 16 days. Of these, 7 are from the 2010 eruption on Merapi, the biggest event during the time period. Table 11 provides the affected nightlight cells by year and volcano and the results are closely correlated with the events of the period. The main one is the Merapi eruption in 2010, Soputan with volcanic explosivity index events of level 2 and 3 in , 2007(Global Volcanism Program, 2013 and Sinabung with several eruptions in 2014, which all fit the model well.  Table 12 refer to province impacts, and the eruptions affected numerous provinces on Java and Sumatra, with Jawa Barat and Jawa Tengah being the most affected with more than 120,000 cells with an SO 2 -value above 10 at some point. This is further underlined by nine of the ten most affected districts in Table 13 being in these two provinces, which are linked to the Merapi eruption in 2010. Apart from that, Sulawesi Utara were affected all the years from 2004 through 2008 by the eruptions on Soputan.
The final table in this section, Table 14, provide descriptives of the SO 2 variable and the nightlight weights, as well as the product of the two. Overall, the mean SO 2 value during these events is almost 20, with a max close to 60, which is 600 times the normal level of 0.1DU SO 2 .   Weights and Intensity multiplied by 1,000

Aggregated data
To aggregate the data, the same method as before is applied, with: where all days t of year T and all nightlight cells i in province or district j are aggregated.
The province overview, shown in Table 16, is caused by the Merapi eruption, with 3 of the top 4 being due to that eruption. Jawa Tengah and Yogyakarta are the two most affected provinces, given their immediate proximity to Merapi. One thing to note is that Jawa Timur, which is east of Merapi, was hardly affected at all.
For districts, the Merapi results are even more pronounced, with all districts in the top 10 being from the 2010 eruption. All the Jawa Tengah districts are west of the volcano. It is somewhat surprising to find that some of the districts in the immediate vicinity of Merapi are not on the list, but this can be due to the timing and quality of the satellite images. Given the time interval between images, the SO 2 clouds could have traveled past the closest districts by the time an image was taken. Regardless, the results are uniform in that all affected districts are neighbors. Overall, the model seems to give a fair picture of when and where the eruptions were at their most intense, although the ground level intensity can be hard to specify.

Tsunami Damage Index
The final disaster damage index constructed is for the 2004 Christmas tsunami. There is little local district level damage data available, so it was decided to use the methodology from Heger (2016), whereby inundation maps are used to construct a district level damage index assuming a uniform damage across all flooded areas.

Damage Index
Despite all the media coverage and attention the 2004 tsunami had, there is not much detailed spatial information readily available for research. Heger (2016) has done some research on the causal effects of the tsunami in his PhD thesis, and we will follow his method closely to model flood impact.
To construct an inundation map of the affected areas, a map based on MODIS satellite pictures from Anderson et al. (2004) is used. The map itself is fairly low resolution, but it provides a good overview of the inundated areas. In terms of the intensity of the flood, there are no existing data on that, but a uniform flood intensity across all flooded areas is assumed, just as in Heger (2016). The resulting map is shown in Figure 10, which shows that a large proportion of the Aceh coastline was struck by the tsunami.
To make this map, the inundation map from DFO was used as a base. Spatial algorithms were then applied to detect the difference in color between inundated and non-inundated areas. This process started with overlaying the base map on a regular shapefile of Indonesia, then detecting the specific color of inundated areas, before constructing a new shapefile where all inundated areas (areas with the same color) have value 1 and all other areas have value 0. , the year after the disaster, there is a strong decline, with 306 and 6,352 lit cells for the inundated areas and Aceh, respectively. The average light intensity has gone down, from 7.39 per cell to 6.33 in the inundated areas and from 6.37 to 5.18 in the province as a whole.

Figure 11: Aceh Nightlights and Tsunami Affected Nightlights
Finally, Table 18 shows the weights which are -again -defined as nightlight in the cell over total nightlight in the province. Although the numbers look very small, by multiplying the means by number of cells, one gets approximately 5.5 percent. Knowing that the census numbers for Aceh in 2000 gave a population of just under 4 million and in 2010 of just under 4.5 million and if one multiplies the population numbers with the affected cells number of 5.5 percent of the total, one gets 221,894 and 249,631, respectively. Given the official numbers of 166,671 dead due to the tsunami, an assumption of total destruction in all inundated areas seem a bit high, given that 166,671 of 230,000 is 72.47 percent. A damage of 75 percent in the inundated cells is chosen, giving a final damage index formula: where TD i is the province weighted damage from nightlight cell i, W i = li

Aggregated data
Aggregating the data is done using the same method as in all previous sections, where the nightlight cells across the province or district is summed up: where all nightlight cells i in province or district j are aggregated.
Since the tsunami only affected one province, it is easy to see the total damage done by it on Aceh. With our assumptions, the tsunami destroyed about 4 percent of the buildings in Aceh. This is clearer once broken down into district level damage. There were 6 districts affected by the tsunami, with Aceh Jaya, Banda Aceh and Pidie all experiencing damage of more than 20 percent. The other 3 affected districts -Aceh Barat, Aceh Besar and Bireuen -experienced damage between 5 and 10 percent due to the tsunami.

Conclusion
With the continuous increase in remote sourcing data, it has gotten much easier and cheaper to monitor and assess the damages from natural disasters. Joyce et al. (2009) and Gillespie et al. (2007) have done an extensive review of how satellite images can be used to map natural disasters, while this paper has contributed with providing new techniques that utilize other remote sourced data such as ShakeMaps and ash advisory data.
Throughout, techniques based on freely available data have been used to construct damage indices for different disaster types. Generally the indices can be used in any area of the world, and if calibrated with local data, they could provide an excellent tool for local governments or stakeholders in early disaster assessments.
The indices can be used to get quick damage estimates and inform where to provide relief, as well as in research such as what the authors have done in Skoufias et al. (2017), where the indices are used to analyze district budget redistributions following natural disasters.
The main caveat is the indices have not been validated against local level damage data. If one had access to high resolution monetary or intensity data, the estimates would be much more precise. Table 20 gives an overview of the different data sources and software used. All disasters have used the DMSP global nightlights data to weight the indices based on economic activity. Recently, the VIIRS nightlight data provide an alternative for assessing economic activity or events, as showcased for GDP in China (Li et al., 2013;Shi et al., 2014) and Africa (Chen & Nordhaus, 2015) and for storms and floods in the US (Cao et al., 2013;Sun et al., 2015). If one is interested in events after 2012, the VIIRS data provide higher spatial resolution and track changes by month instead of yearly.
Alternatives do exist for choice of software if one wants everything to be open source and free. Instead of ArcGIS, packages such as QGIS and RGDAL for R are potential alternatives. R and Python are used more and more for spatial data and can often provide the necessary tools to do the modeling. Instead of statistical software such as Matlab and Stata, R and Python once more provide excellent alternatives.