Multispectral Transform and Spline Interpolation for Mapping Flood Damages Paolo Villa, Student Member, IEEE, and Marco Gianinetto Remote Sensing Laboratory of the DIIAR Department (Hydraulics, Environmental, Surveying Engineering) Politecnico di Milano University, Piazza Leonardo da Vinci 32, 20133 Milano, ITALY email: {paolo.villa; marco.gianinetto}@polimi.it Abstract— In this paper are described some enhancements for a straightforward method recently developed by the authors for evaluating post flood damages using Landsat TM/ETM+ data integrated with Digital Terrain Models (DTMs) and based on the Principal Components Transform (PCT). In particular, the main updates refer to the computational scheme in deriving the flood extension map with the use of the Minimum Noise Fraction (MNF) transformation and in estimating the peak water height and the volumes involved through the use of spline interpolators. Compared with the PCT-based method, the enhanced technique results in improved mapping accuracy. Keywords-Minimum Noise Fraction (MNF) Transform; Spline Interpolation; Remote Sensing; Flood Damages
I. INTRODUCTION Remote Sensing (RS) techniques are well known as the most efficient tool for monitoring the environment at global scale, thanks to the wide availability of data and their easy use. On the contrary, as regards RS applications at regional and local scale, researchers have not yet proved in a completely satisfactory way the usefulness and competitiveness of satellite-based methods compared with ground measures and aerial surveys [1], [2]. Above all, in the field of natural hazards RS potentialities are huge and not yet totally explored. Recent works of the authors have focused on the development of an efficient image processing chain for analyzing and mapping damages due to floods using remotely sensed satellite images, both at local and regional scales [3]. The basic idea is to highlight the environment changes by processing, through the use of a change detection algorithm, imagery collected before and after the flood, together with a Digital Terrain Model (DTM). In previous works [3], [4], the Spectral-Temporal Principal Component Analysis (STPCA) has been used for delimiting the flooded areas, obtaining an Overall Accuracy (OA) with ground truth data between 85% and 90%, while an accuracy greater than 80% has been obtained in the water level and volume estimation through the use of a second order polynomial interpolation. This paper presents some minor improvements of the processing technique regarding the flood map generation and the water level estimation. With respect to the first issue, a new data transformation is introduced, replacing the PCA with the Minimum Noise Fraction (MNF) transformation, while the
peak water height chart is obtained through different interpolation methods, among which spline showed the best results. A case study is presented for the flood that struck the Piemonte Region (North-West part of Italy) in 1994. Four of the most important Italian rivers flow through the Piemonte Region: the Po (460 m3/s of mean annual discharge), as well as several of its tributaries, among which there are the Dora Batlea, the Dora Riparia and the Tanaro rivers (129 m3/s of mean annual discharge alone). During the first days of November 1994 the Piemonte Region was struck by the worst flood event in the century: in four days (from November 4 to November 7) the accumulated precipitations reached the record values of more than 100 mm of rainfall over the whole region and more than 200 mm in the upper part of the Tanaro river basin, with a maximum hourly intensity of 55 mm/h. The consequent flood caused 44 casualties and over 2,000 people homeless. All the inhabited centers placed along the Tanaro river have been inundated and reported serious damages to buildings and infrastructures [4]. In particular, the area investigated in this study is the valley of the Tanaro river, between the cities of Asti and Alessandria, one of the most heavily hit by the flooding. II.
DATASET DESCRIPTION
Two Landsat-5/TM scenes were used for the analysis of post-flood damages, one preceding the flood event and one subsequent. The closest pre flood cloud free image was collected on October 16, 1994 and the first useful post flood image on January 5, 1995, two months after the end of the inundation. Such a temporal lag is the most challenging side of this analysis: applying a methodology for flood detection when flooding waters have completely flowed away. The choice of the Landsat/TM data was considered a good compromise between spatial resolution (30 meters ground resolution), spectral content (6 reflective bands: visible, near-infrared and shortwave-infrared) and area covered by the images (about 35,000 squared kilometers for each scene). Among the further available data there was the DTM of the Piemonte Region supplied in ASCII format by the Cartographic Office of Piemonte Regional Government. Below are listed in details the main characteristics of the dataset used:
0-7803-9510-7/06/$20.00 (c) 2006 IEEE
•
•
Two 30-meters Landsat-5/TM (WRS2 path 194, row 28/29), centered at 44°53’ North latitude and 8°19’ East longitude: o
October 16, 1994 (pre flood image)
o
January 5, 1995 (post flood image);
50-meters gridded DTM derived from the digitalization of existing 1:10,000 scale digital cartography, centered at 44°51’ North latitude and 8°16’ East longitude, with 1 m vertical resolution and 2.5 m Root Mean Square Error (RMSE) vertical accuracy.
III. FLOOD MAPPING The first part of this analysis consisted in the investigation of the maximum extent of the flooding through the generation of a flood map of the area affected by the event at issue. For achieving such an aim it is necessary to separate water related features (linked to inundated terrain), from non-water related ones, through the examination of the spectral differences between the two cited classes of land cover. Pre flood and post flood Landsat-5/TM scenes were first georeferenced in the UTM-WGS84 F32N projection using a first order polynomial warping, being an already geocoded Landsat-7/ETM+ image taken on April 2003 the base for registration. At-sensor radiance data were atmospherically corrected using a low resolution MODTRAN4 model (15 cm-1) combined with aerosol retrieval based on band reflectance ratios and with adjacency correction of path radiance [5], [6], this way drawing scaled ground reflectance data. A synthetic 12-bands file was then created, including first the six reflective bands of the pre flood scene (from band TM1 to band TM5, plus band TM7) followed by the six homologous bands of the post flood scene. To this Spectral-Temporal merge image it was applied the MNF Transform. The MNF is a linear transformation which turns multivariate data with different Signal-to-Noise Ratio (SNR) into a new set of uncorrelated variables, rescaling the noise component to obtain MNF bands laid in decreasing order of SNR, that is to say with increasing image quality [7], [8]. This implementation of the MNF transform is directly deduced from the STPCA, a ChangeDetection technique widely used in RS applications [9], resulting in a Spectral-Temporal MNF (STMNF) operation. After a first of visual inspection of the STMNF results, the STMNF band 2, among which with the higher SNR, was chosen to be representative of a good identification of water related land cover, such as flood water, wet soil and deposited sediments. A threshold was then applied to STMNF band 2 (4%) extremely unlikely to be covered by water, in order to obtain an inundation map of the Tanaro river between the cities of Asti and Alessandria. Both the conditions on the MNF threshold and the terrain slope threshold have to be satisfied at the same time (logical AND).
Figure 1. Flood map derived with STMNF transform and DTM filtering, showing the accordance with ground truth data (both omission and commission errors).
Eventually, the flood map was refined with classical segmentation and clumping techniques to boost the spatial coherency and homogeneity of the final mapping. In Fig. 1 is shown the flood map derived from the processing of the STMNF band 2 compared with ground truth data. Both the errors in omission and commission are shown. IV. WATER LEVEL COMPUTATION Using the flood extension map produced (Fig. 1), the maximum heights reached by the water during the event were then derived by means of interpolating the water level gathered at the border of the flooded area (spatial density of 0.7 points/km2). Along the borderline separating the flooded and non flooded areas a set of terrain heights have been collected through the superimposition of the DTM on the flood map, considering the value of terrain height at the border as coincident with the maximum level reached by the flooding waters. This randomly collected water level data were then interpolated to derive the surface approximating the peak water level. In such a way it was possible to estimate the water height for all the flooded area as the difference between the water level (computed from the interpolation of sparse data) and the terrain level (derived from the DTM). Various interpolation algorithms were tested: 1) linear regression; 2) polynomial regression; 3) and bicubic splines as local component field interpolation [10], [11]. In detail, the four interpolators used in the analysis are listed below and the corresponding water height results are showed in Fig. 2: i.
Linear Regression (LR): direct approximation of the water level with a plane surface (R2=0.94);
ii.
Third Order Polynomial Regression (P3R): direct approximation of water level with a third order polynomial surface (R2=0.97);
iii.
Composed Spline Interpolation (SP2): second order polynomial regression (R2=0.96) for the approximation of the regional component of the
0-7803-9510-7/06/$20.00 (c) 2006 IEEE
heights field and modelling of the residual component through loose bicubic natural splines with Tychonov regularization on second derivative, with the following parameters:
iv.
o
Regularization parameter (µ) = 0.5;
o
Low spatial density: mean mutual distance among radial basis functions centers of about 1,000 meters (416 splines functions over the dominion).
Composed Spline Interpolation (SP3): third order polynomial regression (R2=0.97) for the approximation of the regional component of the heights field and modelling of the residual component through dense bicubic natural splines with Tychonov
(a) Linear Regression Interpolation (LR)
regularization on second following parameters:
derivative,
with
the
o
Regularization parameter (µ) = 0.1;
o
Medium spatial density: mean mutual distance among radial basis functions centers of about 500 meters (1,581 splines functions over the dominion). V.
RESULTS AND ACCURACY
The mapping accuracy was tested using as ground truth data in-situ hydrometric levels taken during the days subsequent the rainfall and an inundation map produced from aerial and local survey, both supplied by the Cartographic Office of Piemonte Regional Government.
(b) Third Order Polynomial Interpolation (P3R)
(c) Second Order Polynomial and Bicubic Splines (SP2)
(d) Third Order Polynomial and Bicubic Splines (SP3)
Figure 2. Water height charts showing the peak level reached during the flood, derived with four different interpolation methods: (a) Linear Regression (LR); (b) Polynomial Regression (P3R); (c) Composed Interpolation with second order polynomial regression for regional component and loose bicubic splines interpolation with regularization for local part (SP2), (d) Composed Interpolation with third order polynomial regression for regional component and dense bicubic splines interpolation with regularization for local part (SP3).
0-7803-9510-7/06/$20.00 (c) 2006 IEEE
The flood map generated from the Landsat/TM data showed an OA of 97.2%, with a Kappa coefficient of 0.82. The best results for the water level retrieval, reached with the second order polynomial and bicubic spline interpolations, (SP2), showed an absolute residual mean value of 0.79 m. In Fig. 3 are showed the statistics for the water height maps in terms of mean value and standard deviation for the absolute residuals computed over 55 Ground Control Points (GCPs), in which peak water levels were known from local survey. Finally, the derived estimated total volume of water and sediments carried by the flooding, on a regional scale showed an agreement of 83% with in-situ measurements, with an overestimation of about 17% of the total.
ACKNOWLEDGMENT Authors wish to thanks the Italian National Research Council (CNR-IREA, Milano) for supplying the Landsat7/ETM+ image taken over Piemonte (Italy) on April 24, 2003 (path 194, row 28/29) that was used for georeferencing the Landsat TM-5/images. Thanks also to the the Cartographic Office of Piemonte Regional Government for supplying the DTM of the Tanaro basin and the in-situ ground truth data. Special thanks to Dr. Mirko Reguzzoni (Politecnico di Milano University-Polo Regionale di Como) for supplying the software (GeoSPLINTER v.1.0) for the spline interpolation. REFERENCES
CONCLUSIONS The method presented has proved to be very effective in deriving inundation maps and local water heights. This makes the methodology a very useful tool to derive a fast mapping of inundated areas and to obtain an accurate chart of the peak level reached by waters during the flood. The use of the STMNF instead of the STPCA and the use of a second order polynomial regression for regional component with a loose bicubic spline interpolation with regularization for local part (SP2), instead of a simpler second order polynomial regression, lead to and overall accuracy increment in flood detection. Enrichment with GIS datasets of local infrastructures (e.g. linear and punctual infrastructures) and land use could help during the post flood decisional phase in affected areas, resulting in a rapid response method of analysis and a decisionsupport tool for governments, local administrative units and also insurance companies.
Figure 3. Histogram and line graph for accuracy displaying of water height interpolation results: the best performance is for the second order polynomial and loose bicubic splines interpolator (SP2). LR = Linear Regression; P3R = Third Order Polynomial Regression; SP2 & SP3 = Composed Spline Interpolation.
[1]
Li Jianping and Zhang Bai, “A GIS-based study on the model for evaluation of direct submerging damage of flood disaster”, Geoscience and Remote Sensing Symposium, 2005. IGARSS '05. Proceedings. 2005 IEEE International, vol. 3, pp.:1807-1810, 25-29 July 2005. [2] H. Bach, F. Appel, K. Fellah, P. de Fraipont, “Application of flood monitoring from satellite for insurances,” Geoscience and Remote Sensing Symposium, 2005. IGARSS '05. Proceedings. 2005 IEEE International, vol. 1, pp.:63-66, 25-29 July 2005. [3] M. Gianinetto and P. Villa, “Monsoon Flooding Response: A MultiScale Approach To Water-Extent Change Detection”, ISPRS Mid-term Symposium 2006 “Remote Sensing: From Pixels to Processes”. The International Archive of the Photogrammetry, Remote Sensing and Spatial Information Sciences, unpaginated CD-ROM, 8-11 May 2006. [4] M. Gianinetto, P. Villa, and G. Lechi, “Post-flood damage evaluation using Landsat TM and ETM+ data integrated with DEM”, IEEE Transactions on Geoscience and Remote Sensing, vol. 44, no. 1, pp. 236-243, January 2006. [5] M.W. Matthew, S.M. Adler-Golden, A. Berk, S.C. Richtsmeier, R.Y. Levine, L.S. Bernstein, P.K. Acharya, G.P. Anderson, G.W. Felde, M.P. Hoke, A. Ratkowski, H. Burke, R.D. Kaiser, and D.P. Miller, "Status of Atmospheric Correction Using a MODTRAN4-based Algorithm," SPIE Proceeding, Algorithms for Multispectral, Hyperspectral, and Ultraspectral Imagery VI, vol. 4049, pp. 199-207, August 2000. [6] A. Berk , G. P. Anderson, L. S.Bernstein, P. K. Acharya, H. Dothe, M. W. Matthew, S. M. Adler-Golden, J. H. Chetwynd, Jr., S. C. Richtsmeier, B. Pukall, C. L. Allred, L. S. Jeong, and M. L. Hoke, "MODTRAN4 Radiative Transfer Modeling for Atmospheric Correction," SPIE Proceeding, Optical Spectroscopic Techniques and Instrumentation for Atmospheric and Space Research III, vol. 3756, pp.348-353, October 1999. [7] A. A. Green, M. Berman, P. Switzer, and M. D. Craig, “A transformation for ordering multispectral data in terms of image quality with implications for noise removal,” IEEE Transactions on Geoscience and Remote Sensing, vol. 26, no. 1, pp. 65-74, January 1988. [8] T. M. Tu, H. C. Shyu, Y. S. Sun, and C. H. Lee, “Determination of data dimensionality in hyperspectral imagery – PNAPCA,” Multidimensional Systems and Signal Processing, Kluwer Academic Publisher, vol. 10, pp. 255-273, 1999. [9] R. A. Schowengerdt, Remote Sensing-Models and methods for image processing. Second Ed. Academic Press, 1997. [10] M. A. Brovelli, M. Cannata, “Digital terrain model reconstruction in urban areas from airborne laser scanning data: the method and the example of the town of Pavia (northern Italy),” IAPRS Proceedings, Commission II, vol XXXIV, part 2, pp.43-48, 20-23 August 2002. [11] H. Mitasova, L. Mitàs, “Interpolation by regularized spline with tension: I. Theory and implementation,” Mathematical Geology, vol. 25, no. 6, pp.641-655, 1993.
0-7803-9510-7/06/$20.00 (c) 2006 IEEE