Spatially distributed rainfall information and its potential for regional landslide early warning systems

Crucial to most landslide early warning system (EWS) is the precise prediction of rainfall in space and time. Researchers are aware of the importance of the spatial variability of rainfall in landslide studies. Commonly, however, it is neglected by implementing simplified approaches (e.g. representative rain gauges for an entire area). With spatially differentiated rainfall information, real-time comparison with rainfall thresholds or the implementation in process-based approaches might form the basis for improved landslide warnings. This study suggests an automated workflow from the hourly, web-based collection of rain gauge data to the generation of spatially differentiated rainfall predictions based on deterministic and geostatistical methods. With kriging usually being a labour-intensive, manual task, a simplified variogram modelling routine was applied for the automated processing of up-to-date point information data. Validation showed quite satisfactory results, yet it also revealed the drawbacks that are associated with univariate geostatistical interpolation techniques which solely rely on rain gauges (e.g. smoothing of data, difficulties in resolving small-scale, highly intermittent rainfall). In the perspective, the potential use of citizen scientific data is highlighted for the improvement of studies on landslide EWS.


Introduction
Rainfall-induced landslides pose a threat to people and infrastructure around the world (e.g. Guzzetti et al. 1999;Crosta and Frattini 2003;van Westen et al. 2008;Günther et al. 2013). Rising temperature (e.g. Gobiet et al. 2014) and rainfall variability (e.g. Casty et al. 2005) in the European Alps make it of great importance to increase efforts in dealing with the consequences of global change on landslide susceptibility. Albeit advancements in the past within the field of geotechnical engineering have led to an increasing in situ damage control in many parts of the world, landslides triggered by heavy rainstorms still cause great losses where no protective structures are available or where they have not been appropriately designed. Indeed, fatal landslide-triggering rainstorm events are occurring in many places around the globe. Recent examples include Messina 2009 (Lombardo et al. 2014) or Rio de Janeiro 2010 (Calvello et al. 2015). Besides those landslides occurring in natural, undisturbed conditions, also established geotechnical mitigation measures are not sufficient and fail in many cases. Although landslides are usually restricted to local sites when investigating single events, they can indeed occur in clusters in the aftermath of a storm event and consequently they can be regarded as a regional phenomenon at specific times (Jaedicke et al. 2014). Floods, on the other hand, are not locally restricted either; however, their spatial occurrence is topographically foreseeable and controllable. This situation shows the necessity of landslide early warning systems (EWS) (Glade and Nadim 2014;Thiebes and Glade 2016).
The observation of intense rainfall events on small spatio-temporal scales is crucial for the development of a landslide EWS (Segoni et al. 2009;Thiebes et al. 2013). The predominant approach in implementing rainfall data into landslide EWS is the employment of empirical rainfall thresholds. This requires the precise knowledge of total precipitation accumulated in a given period, or the rate of precipitation in a period, most commonly measured in millimetres/inches per hour (Guzzetti et al. 2007). A certain threshold is then defined for a given spatial extent by determining the rainfall amount that triggered landslides. Hereby, it is important to differentiate the lowest and highest landslide-triggering rainfall thresholds as the lower and upper limit, below which landslides were never reported and above which landslide were always reported (Gariano et al. 2015).
One major challenge within this approach is, how the recorded landslide-triggering rainfall events were selected. In almost all cases, the precise amount of a landslidetriggering rainfall at a certain location remains unknown. In practice, rain gauges with the closest proximity to a landslide location or which provide the best representation of a certain region are selected for determining a landslide-triggering rainfall event. An indepth consideration of the accurate spatial distribution of rainfall is often neglected for landslide EWS (Thiebes and Glade 2016).
Consequently, the aim of this paper is to provide a useful basis for real-time spatiotemporal rainfall data and to show its potential integration in a regional landslide EWS. Instead of assuming uniform rainfall over a certain area, an automated geostatistical approach is presented. This allows for an approximation of spatially distributed, hourly rainfall predictions in real time based on gauged rainfall data available on the Internet.
2 Review of current approaches 2.1 Landslide early warning systems The UNEP (2012) suggests four key elements which are required for any operational EWS: (a) a comprehensive assessment of the risks, (b) the implementation of monitoring and predicting capabilities, (c) a reliable, synthetic and simple dissemination strategy and (d) development of response strategies combined with the need for raising public awareness and education. Besides financial constraints, however, many operational landslide EWS are discontinued due to the fact that only the purely technical components are considered, and the social aspects such as the elements (c) and (d) (UNEP 2012) are not addressed (Baum and Godt 2010). Not many cases of an integrated approach have been published (e.g. Thiebes 2012;Heil et al. 2014).
Many studies describe the calculation of rainfall thresholds best suitable for a given region (e.g. Glade et al. 2000;Guzzetti et al. 2007;Aleotti 2004). The region under consideration for a certain threshold may range from local scale (Keefer et al. 1987) to global scale . Commonly used thresholds are either attributing rainfall directly or consider additionally the antecedent rainfall. The first one calculates rainfall thresholds for a certain region based on rainfall intensity over a certain time span that caused landslides (e.g. intensity-duration threshold), the latter takes the antecedent rainfall before slope failure into account (e.g. antecedent rainfall threshold), thus considering indirectly also the soil moisture conditions. Literature reports a landslide EWS in Hong Kong implemented in 1977 probably as the first one of its kind (Chan et al. 2003). The first operational landslide EWS in the USA were implemented in the early 1980s (Baum and Godt 2010). A notable mention is the landslide EWS in the San Francisco bay area in 1985 (Keefer et al. 1987) which is said to be the pioneer of modern landslide EWS in combination with rainfall thresholds (Stähli et al. 2015). Ultimately, every successful landslide EWS is a function of expectations that does not only satisfy the requirements of the developers, but also the decision-makers, end-users and the regulations of the prevalent legal system (Thiebes and Glade 2016).
Based on reviewing over 50 EWS for mass movements in Switzerland, Stähli et al. (2015) suggested the following classification scheme for EWS: (a) Alarm systems are directly coupled to sensors that immediately release an alarming upon exceeding a predefined threshold. The accuracy of the prediction is high and the lead time very short (e.g. rockfall). (b) Warning systems detect significant changes in time-dependent factors before an event occurs (e.g. a rainfall threshold). The initial alert is based on predefined thresholds too but ultimately the alarm is released after an expert evaluation. With an extended lead time, warning systems are mainly used for processes with progressive stages of failure (e.g. rock slides, translational and rotational soil slides). (c) Forecasting systems do not use any predefined thresholds, but the level of danger is evaluated in regular intervals. Experts are issuing warnings based on modelling results or sensor data (e.g. degree of avalanche danger for the next day).
Those thresholds involved in landslide EWS can be derived by different approaches. Alarm systems can be based on direct displacement measurements (e.g. Intrieri et al. 2012), while warning systems and forecasting systems obtain their respective thresholds based on rainfall measures (e.g. Capparelli and Tiranti 2010;Aleotti 2004) or deterministic models (e.g. Montrasio et al. 2014;Chien et al. 2015). The most common types applied in landslide EWS are basic rainfall thresholds. Dealing with thresholds always implies a certain amount of a priori knowledge. This knowledge is not readily available in most cases, especially with respect to the precise date and time of landslide occurrence (van Westen 2006(van Westen , 2008 and its respective rainfall conditions that ultimately triggered the failure (Wieczorek and Guzzetti 1999).
When dealing with landslide EWS and having the various data constraints in mind, it is evident that it is impossible to cover all early warning situations (Casagli et al. 2010). Although the calculation of rainfall thresholds is the most common approach for assessing regional landslide early warning situations, they only represent a simplification of the physical processes involved (Reichenbach et al. 1998). This means that in most cases there is more than just this one causative factor (rainfall) involved (Huang et al. 2015).

Observation of rainfall and its spatial representation in landslide EWS
The most important aspect of any landslide EWS is the precise and timely determination of rainfall, no matter whether the approach of the system is based on empirical rainfall thresholds or combined hydrological and slope stability modelling in a landslide forecasting framework. Most researchers are aware of the importance of the spatial variability of rainfall in landslide research. In applications, however, it is commonly not addressed due to the scarcity of available rain gauges (Chiang and Chang 2009). With respect to rainfall thresholds, common approaches to determining the actual rainfall for a single landslide site or a certain region are to utilize single rain gauges near a specific landslide site (e.g. Capparelli and Tiranti 2010) or selected rain gauges as representative locations for a predefined region (e.g. Segoni et al. 2015;Rosi et al. 2015). Aleotti (2004) pointed out that using representative rain gauge locations is dependent on the distance not only from the landslide, but also from other settings such as elevation, aspect or the wind direction. Lagomarsino et al. (2013) present a more advanced approach in their SIGMA warning system by not using a single rainfall threshold for the entire territory of Emilia Romagna (Italy). Instead, they artificially split the region into smaller spatial units, socalled territorial units (TU), which is each characterized by its own threshold. However, each TU is characterized by single rain gauges, thus only offering uniform rainfall for the entire area. With respect to the amount of rain gauges used to characterize aerial rainfall, Guzzetti et al. (2004), for example, use seven rain gauges for an area of ca. 5500 km 2 for their rainfall-induced landslide studies. Similarly, Bathurst et al. (2006) use five for ca. 500 km 2 or Schwab et al. (2008) just one rain gauge for 324 km 2 .
Another seemingly very attractive approach to determining continuous rainfall fields is the utilization of radar technology. Radar technology is capable of offering spatially varying rainfall fields at temporal resolutions of 10 min (and less) and spatial resolutions of around 1 km, which is highly desirable for landslide studies (Chiang and Chang 2009). In addition, there is also vertically differentiated information on precipitation values available. The quantification of rainfall estimates using Doppler radar can provide a real-time comparison with rainfall thresholds to form the basis of landslide warnings (Wieczorek and Guzzetti 1999). Radar data may capture the spatial variation of rainfall fields better than gauged data, especially in mountainous regions where rain gauges are sparsely distributed (Yang et al. 2004;Segond et al. 2007). At the same time, however, terrain obstacles may cause (partial) beam blockage (or beam shielding) when radars are operated in mountainous regions. Beam blockage correction offers ways to enhance rainfall estimates in mountainous regions, but sometimes beam occlusion cannot be prevented in a region of complex terrain (Lang et al. 2009;Anagnostou et al. 2010). On the other hand, radar coverage is not available everywhere (Chappell et al. 2013). What also must be kept in mind is that radar is not capable of measuring precipitation directly. Instead, the rainfall magnitudes are derived from the magnitude of measured reflectivity from the hydrometeors in the atmosphere (Harpold et al. 2017). This procedure requires calibration with rain gauge data in order to adjust the high-resolution spatial pattern with actual measured rainfall data. There is an abundance of literature dealing with this merging process (e.g. Collier et al. 1983;Jewell and Gaussiat 2015;Velasco-Forero et al. 2009;Berndt et al. 2014).
Within the landslide community, the potential of radar data remains often underexploited. There are only few studies implementing spatially continuous rainfall data from radar technology for landslide research (e.g. Crosta and Frattini 2003; Chiang and Chang 2009), even fewer for early warning purposes in combination with numerical weather prediction models and radar data (e.g. Schmidt et al. 2008;Segoni et al. 2009). Although the technology has made great advancements in recent years, the associated uncertainties and generally low success rates (in terms of correctly predicted landslide occurrences) reduce its application within the landslide community. Schmidt et al. (2008) proposed a coupled regional forecasting system in New Zealand based on multiple process-based models (weather forecast, soil hydrology, slope stability). Segoni et al. (2009) proposed a similar approach, yet with the same conclusion that the meteorological uncertainty has the highest influence on the final Factor of Safety map that serves as the basis for landslide early warning. The hydro-meteorological community is far more involved with the utilization of radar data for the deduction of continuous rainfall fields. However, they share the same concerns (Jasper et al. 2002) or even have taken a step further by implementing numerical weather predictions (Cloke and Pappenberger 2009).
Another emergent technology for assessing the spatial extent of rainfall is satellite precipitation data, especially in regions without rain gauges and ground radar coverage. Satellite-based rainfall estimates provide synoptic estimates of the spatial distribution of precipitation events (Chappell et al. 2013). This information is provided at 0.5-3-h intervals at spatial resolutions between 0.07°and 0.25° (Joyce et al. 2004;Kubota et al. 2007;Huffman et al. 2007Huffman et al. , 2010. Tian and Peters-Lidard (2010) questioned the overall accuracy of satellite rainfall products; however, NASA's Tropical Rainfall Measuring Mission (TRMM) has collected precipitation data for over 17 years by now, which is a valuable data archive for global precipitation data (Kirschbaum and Patel 2016). The popularity of satellite-based precipitation data experienced a considerable rise in February 2014 when NASA launched the Global Precipitation Mission (GPM) as a follow-on mission to TRMM (Harpold et al. 2017). GPM provides rainfall and snowfall estimates every 3 h with sensors that are far more advanced and permit better quantification of the physical properties of precipitation particles (Hou et al. 2014). With respect to landslide research, there are only few studies available that utilize satellite-based precipitation data (e.g. Rossi et al. 2012;Kirschbaum et al. 2015), all of which still rely on TRMM data. It will be interesting to observe how near real-time GPM data with higher spatial resolution will be adopted by the landslide community in terms of early warning applications (Stanley et al. 2017).
What has not been covered in detail for landslide studies so far is the wide topic of spatial interpolation for the generation of continuous rainfall fields, which will be addressed in the next chapter as the focus of this study.

Spatial interpolation techniques for the generation of continuous rainfall fields
Spatial prediction is almost always based on samples, but in reality the measurements represent a continuum in space from which the samples have been taken (Oliver and Webster 2014). Historically, spatial prediction was undertaken by purely mathematical approaches that considered only systematic or deterministic variation, but not any error.
Geostatistical prediction, namely kriging, is the logical successor that overcomes most of these drawbacks (Webster and Oliver 2001). Geostatistics aim explicitly at correctly portraying spatial variation (Srivastava 2013), or more precisely, kriging has become a term for several related least-squares methods that provide best linear unbiased predictions (BLUP), in which best is meant in the sense of minimum variance (Oliver and Webster 2014). The atmosphere, or any other environmental feature, is the sum of all kinds of physical, chemical or biological interactions. Although physically determined, they still remain more or less a black box due to its complex interactions that are not fully understood, thus making the variation appear to be random (Oliver and Webster 2014). Consequently, many environmental variables, such as rainfall, can be considered as spatial random variables.
There are many studies available that compare different interpolation methods for assessing spatial rainfall distribution, the majority from the hydrological or hydro-meteorological communities (e.g. Ly et al. 2013;Mair and Fares 2011;Schuurmans et al. 2007;Haberlandt 2007;Goovaerts 2000). Most studies focus on monthly or annual rainfall estimates (Ly et al. 2013). There is only a small number of studies available that uses hourly time steps for the spatial interpolation of rainfall (e.g. Haberlandt 2007;Velasco-Forero et al. 2009;Schiemann et al. 2011;Verworn and Haberlandt 2011). Very few studies are available for landslide applications (e.g. Chiang and Chang 2009). The literature suggests a differentiation between deterministic and geostatistical approaches. The most frequently used deterministic methods for estimating rainfall are Thiessen polygons and inverse distance weighting (IDW) (Ly et al. 2013). The Thiessen polygon method (also Voronoi polygons, Dirichlet tessellation) is one of the earliest and simplest techniques. The region sampled is divided into polygons by perpendicular bisectors between the sampling locations. In each polygon, all points are nearer to its enclosed sampling point than to any other sampling point (Webster and Oliver 2001). For example, Godt et al. (2006) used this technique to characterize rainfall for shallow landsliding in Seattle (USA). The IDW method is rather popular among deterministic spatial interpolation techniques. It is based on inverse functions of distance with the result that unknown locations to the sampling point carry larger weight than those further away. The advantage of weighting by inverse squared distance is the quick diminishing of the relative weights with increasing distance, making the interpolation sensibly local. However, the selection of the weighting function is arbitrary; also there is no indication of error (Webster and Oliver 2001). Chiang and Chang (2009), for example, used IDW to characterize the spatial rainfall distribution for modelling rainfall-induced landslides.
With respect to geostatistical approaches, kriging is the predominant method by connecting mathematical concepts with geoscientific requirements. Kriging is a generalized least-squares regression technique that offers accounting for the spatial dependence between observations (Schuurmans et al. 2007). A huge benefit that comes with this technique is the provision of a measure of certainty. In kriging, the unknown target variable is estimated using a weighted sum of the available point observations. The weights of the data are chosen so that the interpolation is unbiased and the variance is minimized (Ly et al. 2013). One important assumption is that the process under consideration is stationary. This allows researchers to assume that there is the same degree of variation from place to place and that the covariance between two observations only depends on the distance between these observations (Oliver and Webster 2014;Ly et al. 2013). For kriging applications, a graphical summary is used to analyse and understand spatial variation: the variogram (or more correctly: the semivariogram, as it depicts half the variance of the difference of the covariance, but for the sake of simplicity it is mostly called variogram).
The variogram plots variation as a function of distance. This means that rain gauges, for example, being in close proximity will have data values that are more similar to each other than rain gauges further away. Thus, the variogram is a plot of the average squared difference between pairs of data values and thus the central part of any subsequent kriging prediction. Just as the variance, the units of the variogram are the square of the units of measurement (Srivastava 2013). The plotted variogram based on the sampled data is called experimental variogram (Oliver and Webster 2014). There are three key characteristics of a variogram: (a) sill: the plateau that is reached by the variogram estimating the variance of the random process; (b) nugget: the y-intercept for the unresolved variation; and (c) range: distance at which the variogram reaches the sill (Srivastava 2013). To make the variogram applicable for the kriging prediction, a curve has to be fitted to the experimental variogram to neglect any inherent point-to-point erratic fluctuations. This fitted curve has a mathematical expression that describes the variance of random processes with changing distance and guarantees non-negative variances in the predictions (Oliver and Webster 2014). Commonly used correlation functions deriving the theoretical variogram are exponential, spherical, Matérn or Gaußian. In practice, there are different types of kriging. In the geoscientific literature, common prediction techniques are (a) ordinary kriging (OK), (b) kriging with external drift (KED), and (c) ordinary cokriging (OCK). Applied to spatial precipitation pattern, OK uses only rain gauge information, while the other two techniques incorporate sampled secondary information (e.g. weather radar information or elevation) to improve the kriging prediction (Goovaerts 1997). There is an abundance of related geostatistical literature available that readers are referred to for more in-depth information (e.g. Isaaks and Srivastava 1989;Cressie 1993;Goovaerts 1997;Webster and Oliver 2001).

Obtaining real-time weather information from automated web scraping
This study uses web-based, hourly rain gauge data being obtained in an automated data workflow. Hourly time steps were selected due to the fact that longer observation periods would average rainfall intensity, which is detrimental for landslide early warning purposes due to the underestimation of peak (maximum) rainfall (Guzzetti et al. 2007). In this study, hourly time steps are considered as a compromise that still allows for creating spatially distributed rainfall information in near real-time and reflects the short-term rainfall intensity requirements for early warning applications, although shorter time steps would be even more desirable. In the past few years, environmental sciences are witnessing an increasing amount of published information on the Internet (Vitolo et al. 2015), which includes high-quality real-time weather data from various data sources. In many scientific fields, data obtained from the Internet are being used in a different way, ranging from simple data extraction tasks to fully automated data processing workflows. This growing demand leads many operators of databases and servers that exhibit a certain volume of traffic and where well-profiled usage expectations are available, to design publicly available application programming interfaces, so-called APIs. Data APIs are agreed-on programming interfaces providing a structure to download and link large chunks of heterogeneous data. Although such web services are the standard and recommended way to enable external access, such APIs are not always available, especially when the potential user base, and consequently the related demand on data, is minor. In such cases, web scraping is an alternative and valuable method for extracting and combining content from the Internet in a systematic way, even in the absence of an API (Glez-Peña et al. 2013).
Web scraping mimics the human user interaction with a website in a systematic way. The web scraping application is accessing as many websites as desired, parses its content to find and extract relevant information and structures it in a way ready to use for subsequent analyses. Although there are some desktop-based web scraping solutions (Glez-Peña et al. 2013 list a few), the most common approach is to use any suitable programming language to achieve maximum flexibility. There are many third-party, open source libraries available for implementing them in the source code for developing own web scraping applications. Figure 1 shows a rough overview of the entire web scraping process on how hourly rainfall data from various data sources is fetched, parsed and stored so that it can be retrieved for the generation of spatially distributed rainfall information. To establish the web scraping application to obtain hourly real-time weather data, the JavaScript based MeteorJS web framework was used in this study. It allows for rapid prototyping and the establishment of a server/client application. The scraped website access is executed by accessing the document object model (DOM) of the HTML content. This is done with Cheerio, which is based on the jQuery library to simplify client-side scripting. Cheerio (https://atmospherejs. com/fermuch/cheerio) itself is specifically designed for the Node.js runtime environment. With Meteor's http-package, the HTML content was loaded. To ease parsing, the Node.js module html-parser2 (https://www.npmjs.com/package/htmlparser2) was used. To provide the extracted data also for later usage and not just real-time applications, the obtained information is stored in Meteor's default database (MongoDB) in a structured way (Fig. 2). For automated, time-based scheduling of the scraping process, the synced-cron-package (https://atmospherejs.com/percolate/synced-cron) was used. To access the data directly from the database in the source code of our subsequent analyses, the restivus-package (https://atmospherejs.com/nimble/restivus) was used to create RESTful APIs. The Mon-goDB used can be classified as a NoSQL database that stores data as JSON-like (JavaScript Object Notation) documents which makes it a fast and far more flexible way to store scraped data. Fig. 1 Flow chart of the implemented web scraping application to fetch, parse and store real-time rainfall data for subsequent analyses Nat Hazards

Automated spatial interpolation
There are always reliability concerns with respect to web-based data (de Vos et al. 2016). Although no crowdsourced weather data are included in this study at present, there is still the possibility of erroneous or bogus data involved. The quality issues of rain gauge data come from remotely, automated data sampling and the electronic transmission of the data through several ports before being used in an application (Kondragunta and Shrestha 2006). Therefore, real-time rain gauge quality control (QC) is necessary to detect major inconsistencies. Additionally, single-station QC checks alone are not sufficient due to the high spatial variability of rainfall; thus, adjacent stations have to be taken into consideration too. For this purpose, three automated QC/plausibility checks were implemented in this study: (a) a range filter; (b) a spatial consistency filter; and (c) an autocorrelation filter (for the kriging application only). The range filter is a simple check performed on a single observation for a given location for a specific time. Rain gauges with hourly rainfall intensities below 0 mm (physically impossible) or above 25 mm (very unlikely in that area) are excluded for subsequent analyses. The spatial consistency filter uses a distance matrix to calculate the 95th percentile of all available rain gauge intensities within 20 km (Fig. 3a). This distance was chosen to cover a substantial amount of neighbouring rain gauges and due to pronounced spatial autocorrelation that was shown in the variogram within this distance. If a rain gauge contains an intensity greater than the 95th percentile, it is discarded for subsequent analyses, but only if the intensity is higher than 15 mm. This threshold was chosen due to the fact that predominantly high rainfall intensities are important for landslide early warning purposes while the severe weather centre in Austria (Unwetterzentrale) defines heavy rainfall from 17 mm. Therefore, some tolerance was added to that boundary condition. The spatial consistency filter does not serve purely as a plausibility check, but also as a means to deal with outliers for the geostatistical interpolation, although we cannot know whether an unexpectedly large value is a real outlier (i.e. punctual high-intensity rainfall) or not. Outliers cause serious distortions in a variogram and should be removed if they are suspected to belong to a process other than the one interested (Oliver and Webster 2014). Therefore, we exclude such outliers from the data set to model the spatial relationship between the rain gauges, but use the entire data set for the kriging prediction. This way, the variogram modelling is not corrupted by outliers, but the kriging surface still accounts for the extreme values.
The last implementation is the autocorrelation filter. When examining the distribution of environmental parameters, it often tends to be positively skewed towards the smaller values. This is also the case with hourly rainfall where some regions contain a certain amount of rainfall and other regions do not receive any rainfall at all. So there is a strong clustering of natural zero values in the distribution which makes transformation difficult. The comparison between means of observations is more unreliable because the variances are likely to differ from one set of data to another (Webster and Oliver 2001). Therefore, non-transformed data are used but isolated zero values are treated with the autocorrelation filter with regard to variogram modelling. Schuurmans et al. (2007) also point out some problems arising with zero rainfall for kriging estimates. Geostatistics is based on the premise of autocorrelation. In an area where rainfall is recorded, there is a certain autocorrelation between the rain gauges. In an area without rainfall, zero values are negligible for the purpose of this study. The boundary between an area containing rainfall and an area without rainfall is, however, still relevant because stations close together are still autocorrelated to a rather high degree. The autocorrelation filter has, consequently, three conditions that must be satisfied: (a) keep all stations that have at least 0.1 mm rainfall; (b) remove all stations that have 0 mm rainfall and only have stations with 0 mm rainfall within distance; and (c) keep all stations that have 0 mm rainfall, but have stations with at least 0.1 mm rainfall within distance (Fig. 3b). As with the spatial consistency filter, the search distance equals 20 km. Additionally, the filtered data set is used for the variogram modelling, but all sample points are used for the kriging prediction. Naturally, this automated filter approach does not work in all situations, especially when the amount of remaining rain gauges is significantly reduced (Webster and Oliver (1992) point out the effect of sample size on variogram estimation). Therefore, a set of data with and without this filter is used in subsequent analyses and omitted on poorer validation results. Large areas with no rainfall that had their rain gauges Fig. 3 a The spatial consistency filter removes rain gauges from an hourly data set that contain suspiciously high rainfall intensities compared to surrounding rain gauges; b the autocorrelation filter removes rain gauges from an hourly data set that contain zero rainfall and have only rain gauges with zero rainfall in their vicinity. Rain gauges that have zero rainfall but contain rain gauges within its vicinity that exhibit rainfall, are kept due to very apparent spatial dependency. The x and y axes refer to projected UTM coordinates. Triangles indicate retained rain gauges, squares indicate eliminated rain gauges, large circles indicate search radius for the distance matrix, and numbers refer to rainfall amount (mm h -1 ) removed this way have a significantly higher variance consequently. Figure 4a shows an overview of all filters that are applied to the raw web-based rainfall data.

Geostatistical methods
Exploratory data analysis was the initial step to assess how representative and consistent the available rain gauges are distributed in the study area. For geostatistical analyses, areas with a high point density provide more reliable estimates at unsampled locations than areas with only a few rain gauges. To assess whether the average distance from an arbitrary point to the next nearest sample point is significantly short, testing for complete spatial randomness (CSR) might give an indication (Diggle 2003). The function behind estimating CSR is a nearest neighbour distance distribution function G(r). Testing for CSR covers the horizontal domain. To check whether there are significant differences in distribution with respect to elevation, an overlay sampling of the rain gauge locations was performed on a digital elevation model and tested with a nonparametric KS test. Kriging uses the semivariances based on a fitted variogram function. Creating a theoretical variogram is usually a labour-intensive task as many choices must be made (such as finding suitable range, sill, nugget values). There are attempts to automate the procedure, yet it carries risks when fitting a variogram without surveillance. Cressie (1985), for example, suggests a specific variogram where the weighted sum of squares between experimental and theoretical variogram becomes a minimum. However, the sensitivity of the variogram estimation on the interpolation itself is often not quite high. On the other hand, if a precise estimation of the error variance is needed (e.g. for uncertainty assessment), a solid variogram estimation is required (Haberlandt 2011). Therefore, the focus of this study lies in the rapid estimation of a suitable variogram for real-time applications.
The initial step for the variogram modelling is importing the georeferenced rain gauge data (filtered and unfiltered) from the web scraping application which happens every hour as soon as up-to-date information is available. To proceed, usually the variogram parameters range, sill and nugget should be estimated by exploratory analysis or expert knowledge. To automate this procedure, the decision is made to assign some plausible Fig. 4 Flow chart showing the automated workflow from a web-based data generation, the application of various filters to the raw rainfall data to eventually exclude single rain gauges, through b the interpolation procedures for producing hourly real-time rainfall raster maps initial values. Therefore, the initial range is defined as 0.1 time the diagonal of the input shapefile that defines the boundaries of the study area (the bounding box). The initial sill is calculated as the mean of the maximum and median of the semivariance; the initial nugget is defined as the minimum of the semivariance. The next step uses a loop to iterate over multiple permitted variogram models (Spherical, Matérn incl. M. Stein's parameterization, Gaussian, etc.) to select the model with the smallest residual sum of squares (Oliver and Webster 2014). Additionally, all kappa values are tested for the Matérn model. This procedure with its associated values is based on the automap-package in R (Hiemstra et al. 2009). Omnidirectional variograms are created for both filtered and unfiltered input data. In case there is uniformly no rainfall at all or just very little in just a few sampling locations, the variogram is basically a horizontal line, indicating a pure nugget variogram due to the fact that there is no variance between the samples. The prediction in that case would be everywhere the same and thus the mean of the data (Oliver and Webster 2014). It would be difficult to interpret a pure nugget estimate, but the physical justification was given more weight: uniformly zero, or almost zero, rainfall does not matter too much for the underlying research question, which is more interested in larger rainfall intensities for landslide applications. And those are never uniformly distributed.
Here, a univariate approach is applied for spatial prediction, namely ordinary kriging (OK). Adding auxiliary variables as additional predictors might be helpful in some cases, but Ly et al. (2011) have shown, that adding elevation does not improve interpolation accuracy for short time intervals. Also Haberlandt (2007) found that elevation information plays a minor role for hourly rainfall data. When incorporating additional predictors to multivariate approaches, the benefit might be marginal if correlations become too small, as concluded by Goovaerts (2000). Figure 4b shows a flow chart of the automated kriging procedure used in this study. Using OK, hourly rainfall was estimated at unsampled locations on a 1-km square grid across the study area. The 1-km grid spacing used here was set for consistency and subsequent comparison with the radar data that contains the same grid size. Also according to Hengl (2006), who provided some empirical and analytical rules for the selection of suitable grid sizes, the selection of a 1-km grid size is justified based on the amount of rain gauges available. Kriging estimates might take negative values when negative kriging weights are applied. This is undesirable because this can lead to non-physical estimates. Possible solutions to avoid negative estimates are either a posteriori corrections of the kriging weights, as suggested by Deutsch (1996), or simply replacing all negative values with zero (Ly et al. 2013). In this study, the latter approach was used to achieve a physically sound rainfall estimation.

Deterministic methods
With regard to deterministic estimation procedures, the inverse distance weighting (IDW) and Thiessen polygon methods are used. The IDW method gives each rain gauge a weight that is inversely proportional to the distance between that rain gauge and an unknown sample point. Critical user input to this method is a distance parameter (an exponent) that controls the degree of dependence to a rain gauge in closest proximity (Srivastava 2013). A smaller value gives rain gauges further away higher importance, while accordingly larger exponent values assign higher importance to closer rain gauges. IDW is not capable of providing any quantitative indicator of reliability; therefore, an iterative approach was used to test multiple exponent values between 2 and 5, which was considered a plausible range for avoiding biased estimates. The distance exponent with the highest goodness of fit, in terms of the highest coefficient of determination in the validation process, is ultimately used for the final IDW based rainfall estimate (Fig. 4b). Another deterministic method used to compare kriging and IDW estimates is the Thiessen polygon method. This method simply divides the study area into polygons by perpendicular bisectors between the rain gauge locations. Within a polygon, all unknown points are closer to its enclosed rain gauge than to any other rain gauge (Webster and Oliver 2001).

Radar rainfall estimates
In 1965, Austro Control (the Austrian air navigation services provider) started operational service of its first weather radar at the airport Wien-Schwechat. Since then, many improvements have been made until 2011, when the first dual-polarization ground radar was installed. Those new-generation systems enable the transmission of radio signals with both horizontal and vertical polarization, while the conventional Doppler systems only transmitted and received radio waves with single horizontal polarization (Harpold et al. 2017). Since 2001, radar data in Austria are available with a spatial resolution of 1 9 1 9 1 km and a temporal resolution of 5 min. The operational weather radar network in Austria consists of five stations, each with a range of 224 km covering the entire territory of Austria (Kaltenböck 2012). The used 2d-weather radar composites contain 14 quantization steps given in reflectivity (dBZ). Each radar composite has a spatial resolution of 1 9 1 km and contains data from a 5-min scan. Consequently, the 5-min rainfall reflectivity was summed by hour and divided by 12 for the hourly average to match the information from the rain gauges. To convert the reflectivity Z to rainfall intensity R (in mm h -1 ), the empirical Marshall-Palmer equation was used with the relation Z = 200R 1.6 (Lovejoy et al. 2008). For this study, the radar data remained uncalibrated; thus, the converted rainfall rates from the radar data cannot be directly compared with the rainfall rates from the interpolated rain gauge data, but are used as a qualitative means of validation.

Study area
The study area includes the entire federal state of Lower Austria (Niederösterreich) in the north-eastern part of Austria. The size of the study area is 7408 km 2 . The mean annual precipitation rates in Lower Austria for the period 2001-2010 show a gradient from lower rates in the northeast (approximately 500 mm) to higher rates in the southwest (approximately 1600-1700 mm) . Schweigl andHervás (2009) andSchwenk (1992) mention exceptional rainfall and/or snow melt as main triggers for landsliding in Lower Austria. A recently compiled landslide inventory for Lower Austria based on 1-m LiDAR DTM derivatives and orthofoto mapping reveals 13,166 landslides .
The Cretaceous-Early Tertiary Rhenodanubian Flyschzone (RDF) contains approximately 6300 mapped landslides but contributes only to 14% of the territory of Lower Austria. Although there has been extensive work on statistical landslide susceptibility assessment in Lower Austria (e.g. Petschko et al. 2014;Steger et al. 2015Steger et al. , 2017, there are no published rainfall thresholds or process-based modelling approaches available. To test the proposed methodology, a rainfall event from June 2009 that triggered many landslides in the southern parts of Lower Austria was selected. In total, 92 rain gauges were available to apply different spatial interpolation techniques for the automated generation of real-time continuous rainfall fields (Fig. 5). All rain gauges used in this study are contained in officially operated, automatic weather station networks from different weather service providers.

Variogram modelling
Exploratory data analysis was performed as a first step to check how the rain gauges are spatially distributed in the study area with respect to spacing and elevation. Figure 6a shows the CSR plot of the empirical functionĜ r ð Þ against the theoretical expectations G(r), indicating whether the average distance from an arbitrary point to the next nearest sample point is significantly short. The upper and lower simulation envelopes are calculated based on 100 simulations and indicate significance bands (the number of simulations was chosen arbitrarily but was found justified as no significant changes are expected  Fig. 6 a Testing for complete spatial randomness (CSR) to assess the spatial representativeness of the rain gauge distribution (''sampling design'') with respect to geographical space.Ĝ r ð Þ represents the empirical function (continuous line), G(r) the theoretical expectation (dashed line) within its significance bands (''envelopes''); b relative frequency distribution of rain gauges with respect to overall distribution of elevation in the study area Nat Hazards beyond that). Given the fact that there was no sample design involved in the selection of the rain gauge locations and that the average spacing between the rain gauges is around 10 km, the distribution of the sample locations is quite representative with respect to geographical space. With respect to elevation, using a nonparametric KS test showed that the distribution of rain gauges is not significant with respect to the overall distribution of the elevation in the study area (Fig. 6b). Especially higher elevations and elevations in the 400-500 m range are overrepresented by rain gauges.
Variance in every rainfall event differs; thus, also each variogram is different. From the automated fitting procedure, it can be concluded that generally there is quite a large variability in fitted ranges across hourly events. On a multitude of variograms, however, the range is between 30 and 50 km. But there remain many variograms with a very large range leading to spurious autocorrelations, probably caused by large-scale trends extending throughout the study area. From the iterative model fitting procedure, the most common applied model was the Matérn model with Stein's parametrization (resulting from iteratively comparing the smallest residual sum of squares). In some cases, the filtering procedure also leads to a reduction in sill (Fig. 7).

Real-time rainfall interpolation
To demonstrate the different interpolation techniques described in Methods, three consecutive hours of a frontal rainfall event were selected that caused landslides in the southern parts of Lower Austria. 92 rain gauges served as the basis for the spatial interpolation. Figure 8 shows a composite of those 3 h for all rainfall predictions. All predictions were performed on a 1 9 1 km grid. All predicted rainfall fields are in good visual accordance; however, the Thiessen polygons (Fig. 8e) are considered to grant no real advantages over the conventional representative rain gauge approach to characterizing spatially distributed rainfall for a certain area as the transition between the arbitrary polygon boundaries is unrealistically rough. This is even more extreme for very short time intervals and in mountainous regions where recorded rainfall intensities may vary significantly within short distances. Also in just those 3 h, two problems associated with the kriging technique are apparent. When considering the point information from the rain gauges (Fig. 8a), rainfall in areas with high intensities is much lower for the kriging estimates (Fig. 8b, c). Thus, kriging loses variance by smoothing, yet it gives the best estimates from a statistical point of view. The second issue with kriging is the presence of Fig. 7 Exponential variograms for the same time but without data filtering (a) and with filtering (b). In this case, the filtering leads to a slightly decreased sill which indicates lower variance of the residuals at greater distances  Fig. 8b), a single high-intensity rain gauge is present that leads to a strange, bullseye-like structure. This rain gauge was removed with the automated filtering (Fig. 8c), still it seems that a spatially limited high-intensity rain cell causes a big influence in this part of the study area that cannot be properly resolved by the automated variogram modelling.

Performance comparison for hourly interpolation
To estimate the performance of the OK, filtered OK and IDW predictions, two different types of validation were performed. The leave-one-out cross-validation (looCV) removes one rain gauge at a time and recalculates its value from the remaining data. Validation was also performed by splitting the sample size randomly into a training (80%) and a test (20%) subset. The training subset was then used to predict the values of the test subset. For the three consecutive hours shown here, the validation results are presented in Fig. 9a. The problems in the kriging estimates for the local high-intensity rain cell at 01:00 24 June, 2009, are also reflected in the validation results showing much lower coefficients of determination (around 0.65). The next 2 h produced more balanced kriging estimates resulting also in better validation results (around 0.8). For those 2 h, the application of the different filters (range, spatial consistency and autocorrelation filter) also leads to a slight increase in performance. For every rainfall prediction, a standardized residual plot is generated which measures the strength of the difference between observed and predicted values (Fig. 9c). The residuals from the automated rainfall prediction tend to be symmetrically distributed (homoscedastic) and there are no clear patterns in general, which indicates that the automated model prediction is feasible.
Many studies indicate the root-mean-square error (RMSE) as a performance indicator. We refrain from this practice as it does not give justice to the spatially differentiated variances involved in kriging estimates. We found a spatially distributed representation of kriging variances (the estimation error) to be a more suitable tool (Fig. 9b). This map also reveals much larger variances near the boundaries of the study area due to the reduced Fig. 9 a Validation results (coefficient of determination) for the ordinary kriging (OK) and IDW estimates. (looCV = leave-one-out cross-validation; subsample = sub-sampling into training and test subset; OK_f = OK filtered); b kriging variances for an unfiltered (left) and filtered (right) prediction serving as a spatially distributed error estimation; c automatically generated standardized residual plot indicating constant variance Nat Hazards number of rain gauges in close proximity. Additionally to this quantitative validation, a qualitative validation was performed based on the visual comparison with radar rainfall data. To use the radar imagery for a quantitative validation or as an auxiliary variable for multivariate kriging approaches, it needs to be calibrated first. This would be highly desirable for an additional means of quantitative validation because it provides an independent data set from the interpolation results. Comparing the radar data from Fig. 10 with the interpolated rainfall predictions from Fig. 8, the overall picture shows a rather good match. However, the radar technique is capable of capturing more fine structured, intermittent rainfall fields.

Discussion
This study suggests a fully automated workflow from the hourly, web-based collection of rain gauge data to the generation of spatially differentiated rainfall predictions based on deterministic and geostatistical methods. The underlying research question envisages the implementation of those hourly rainfall predictions into a dynamic, combined hydrological and slope stability modelling application with the purpose of estimating landslide failure probabilities in near real time. The entire methodology was implemented solely with open source technology (JavaScript, R, Python, QGIS).
When web data are used, legal and policy issues have to be considered. Legal implications with web data are not always clear in all cases and countries; however, the terms of use should always be respected to limit the permitted data requests within a certain amount of time or, in general, to prevent copyright infringement (Glez-Peña et al. 2013). The usage of meteorological data is better regulated. The Twelfth World Meteorological Congress in Geneva in June 1995 approved a resolution on the international exchange of meteorological data and products (WMO 1996). This resolution stipulates the member states of the World Meteorological Organization the right to use data and products at no costs for noncommercial use. When using web scraping for data generation instead of an API request, a problem might arise in the way how the data are automatically parsed from a website. As soon as the structure of an HTML document is changed, the web scraping application does not work anymore and has to be readjusted, thus requiring constant maintenance. Using these data for geostatistical analyses is usually a very labour-intensive work.
Modelling the variogram is critical for the quality of the kriging estimates; thus, automating this procedure is not so straightforward and required some simplifications in Fig. 10 Uncalibrated radar data for three consecutive hours used for qualitative validation (hence the deliberately omitted legend). The spatial pattern is in good visual accordance with the interpolation rainfall estimates shown in Fig. 8 Nat Hazards terms of defining initial modelling parameters (i.e. sill, range, nugget). Therefore, a stronger emphasis was put on creating a feasible workflow that produced good quality rainfall predictions for further use. Also worth mentioning is that anisotropy was not considered in this study, although rain gauges recording frontal rainfall are likely to be directionally dependent. However, detecting anisotropy in an automated workflow is not straightforward and could only be approached by iteratively calculating multiple directional variograms and perform kriging predictions for all of them accordingly while comparing their respective validation results. On the other hand, not all rainfall situations are directional (e.g. convective rainfall); therefore, omnidirectional variogram models were used in this study. Additionally, this suggested procedure is only feasible for rainfall and no other forms of precipitation, as only direct rain gauge readings are considered. This is acceptable for this study area due to highest rainfall intensities recorded in the summer months which is relevant for landslide initiation. When passing along those automatically generated rainfall predictions to a landslide modelling application, a solid quality indicator is required as a diagnostic tool beforehand in order to guarantee a good performing model. Or in other words: where to set the threshold between good and bad rainfall predictions. One may use the coefficient of determination or, as suggested by Oliver and Webster (2014), the mean-squared deviation ratio (MSDR) which is the mean of the squared errors divided by the corresponding kriging variances. Possible thresholds would be high coefficients of determination or MSDR values close to 1. However, this has to be extensively tested in the landslide application which one is the most suitable with respect to the underlying research question.
Another major issue commonly reported in other studies (e.g. Schuurmans et al. 2007;Kann et al. 2015) is the presence of small-scale, convective heavy rainfall events, mostly occurring in summer or spatially highly intermittent rainfall (e.g. Chappell et al. 2013). Thus, different types of rainstorms may provide different levels of performance for the proposed methodology. Simple ordinary kriging purely relying on distance-based rain gauge information might not be capable of solving this problem with a limited number of sampling locations alone. In that case, using multivariate kriging approaches that incorporate radar data as additional predictors might be more appropriate. Another interesting approach to mimicking true in situ variability on short-scale changes is the utilization of conditional simulation, or how Srivastava (2013) describes it: the spatial version of Monte Carlo procedures. Kriging usually produces a smoothed surface by losing variance and thus underestimating large values and overestimating small values. Also, kriging estimates produce only a single prediction. Conditional simulation, on the other hand, produces many equally likely scenarios (so-called realizations) by using the same variogram, but at the cost of accuracy. Consequently, this probabilistic procedure produces a more realistic picture of small-scale variations, but should be complemented by kriging estimates when spatial variability and error estimates are crucial for the underlying research question (Srivastava 2013).
With regard to regional landslide EWS, this study offers a direct integration for both threshold-based and process-based approaches. The methodology presented in this study is especially suitable for the implementation in warning systems (following the classification of Stähli et al. 2015) that contain predefined thresholds and are mainly used for processes with progressive stages of failure (e.g. rock slides, translational and rotational soil slides). Instead of using uniformly distributed rainfall for an entire region, spatially differentiated rainfall values can be used for a real-time comparison with previously established rainfall thresholds. Similarly, hourly rainfall predictions can also serve as the time-dependent input in process-based approaches that determine how changes in pore-water conditions alter slope stability conditions. Raia et al. (2014) presented a promising approach where the dynamic, infinite-slope-based model TRIGRS (Baum et al. 2008) was used for the forecasting of rainfall-induced shallow landslides over large regions. In their probabilistic modification of the model, however, they assumed constant rainfall intensity to force slope instability for the entire study area. Salciarini et al. (2017) followed a similar probabilistic approach. It would be interesting to observe how spatially distributed rainfall intensities will behave on slope stability conditions when compared to uniformly distributed rainfall input.

Outlook
A major issue in almost all natural sciences is data scarcity. Especially in landslide research proper event data are often lacking. With the rise of Web 2.0 applications, however, there was a large boost in collaborative and fast data acquisition initiatives. Olyazadeh et al. (2016) present a WebGIS Android App for fast data acquisition of landslide hazard, and Klonner et al. (2016) present a review of how volunteered geographical information is collected in natural hazard analysis. Baum et al. (2014) highlight the Report a landslide website operated by the USGS for engaging public in the identification of geological hazards. With respect to weather data, there are already much longer lasting initiatives in operation. Such citizen science initiatives have proven to be highly valuable for supplementing primary instrumented rain gauge networks (Harpold et al. 2017). The National Weather Service Cooperative Observer Program (COOP) in the USA was formally created in 1890 for volunteers to take observations. There are many websites that provide open weather data collected by weather enthusiasts that even offer APIs for direct data integration (e.g. http://openweathermap.com/ or http://wunderground.com/ with their personal weather station network). Those privately operated weather stations that are hosted by such online weather networks have explicit terms of service that facilitate API data usage.
Quality issues might be a big concern with such data; however, the benefit of offering a highly densified rain gauge network providing rainfall data in real-time is not to be underestimated and should clearly be addressed in the near future. This paper suggests an automated workflow that enables the quick integration of additional real-time rainfall data from multiple online sources with either an API or web scraping integration. Therefore, very dense rain gauge networks could be established for providing accurate spatially distributed rainfall predictions for the integration in regional landslide EWS. Future work will incorporate this approach in dynamic, grid-based regional slope stability analysis to evaluate in a real-word hindcast situation if and to what degree spatially distributed rainfall information in landslide research contribute to improving regional landslide EWS. and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made.