Multi-spectral imaging: a tool for the study of the area burned in wildfires

Multi-spectral images provide powerful information to assess the effects of wildfires. Data can be processed, combined and analyzed in such a way that can contribute to suggest a protocol of action to help recover devastated areas. The objective of this paper is to describe methods to evaluate the surface burned in a forest fire. Combining the information recorded by the visible and infrared sensors of the Landsat 8 satellite, we can estimate the temperature of the surface after the wildfire. Finally, by comparing images recorded before and after the fire, we can estimate the surface burned.


Introduction
Climate change is one of the most serious threats that is facing the planet. Nowadays, largescale forest fires are becoming more frequent, especially in dry areas with vegetation that can be found in Eastern Australia, California, or the Mediterranean countries, among others [1][2][3]. Moreover, even the Amazon rainforest can burn because of forest degradation and deforestation [4].
At the beginning of the summer of 2019, Catalonia suffered a severe heat wave. High temperatures, strong winds and the absence of rain led to the appearance of forest fires. On June 26, a forest fire was declared in the Ribera de l'Ebre county, located 100 km west of Tarragona, burning an area of approximately 50 km 2 [5]. The fire was declared extinguished on July 7. Landsat 8, a satellite operated by the US Geological Survey and NASA, can record multispectral images in 11 bands ranging from the visible to the long-wave infrared [6]. Nowadays, multi-spectral information plays a key role in the analysis of the effects of fire in forests. The objective of this paper is to describe a method to evaluate the surface burned in a forest fire. Starting from satellite images, we estimate the at-sensor temperature using a black body approach. Then, we assess the severity of the burned land using analysis methods on images recorded before and after the fire. Finally, we identify the damaged area using mathematical morphology operations. The values of the estimated area are compared with those reported by official sources.
Multi-spectral imaging is routinely used in the analysis of wildfires (see, for instance, [7] and references therein). Also, this topic is analyzed in textbooks on remote sensing [8] and several practical tutorials [9, 10]. Other authors have considered topics that partially overlap with our work. Interestingly, Caridades and collaborators described educational experiences on basic image processing techniques and its application to the estimation of the area of a surface [11,12]. Moreover, multiple experiments involving infrared cameras have been described; we refer to [13][14][15] to cite a few.
The paper is organized as follows: in section 2, we describe the image set required to perform the calculations. In section 3 we describe the procedure to perform geometrical changes on the images in such a way that no zero values appear. In section 4, we estimate the at-sensor temperature using a black-body approach then, the differential normalized burn ratio is introduced in section 5. This index is useful to classify the severity of the damaged area. Finally, the burned area is calculated in section 6. We present our conclusions in section 7. Moreover, a notebook with the Python code that calculates the images of the paper is available (see supplementary materials).

Motivation
Digital image processing is a basic technique [16] that is widely used in experimental physics. At the Universitat de Barcelona, Image Processing and Computer Vision (IPCV) is a 3 ECTS credit subject intended for junior or senior students in Physics 1 . At this point, students have taken two courses (6 ECTS each) in coding techniques and computational physics and are skilled in Python. Depending on their academical profile (fundamental or experimental physics) they have some extra experience with Fortran or Labview.
IPCV is a problem-based introductory course focused on algorithm implementation. Around 30 highly motivated students take these credits every year. The students attend to lectures, computer labs and develop a computational project of their choice. They perform a bibliographical search, select a proper topic, develop the code and write a journal-like report. The present paper is an extended version of an original student project.

Landsat 8 image set
Satellite images are extremely useful for developing analysis tools that entitle us to understand the expansion of fires. Landsat 8 is an Earth observation satellite, launched on 11th February 2013. Landsat 8 comprises two sensors: the thermal infrared sensor (TIRS) and the multispectral operational land imager (OLI) with 11 different channels available.
Landsat 8 datasets can be downloaded from the USGS Earth Explorer platform [17]. They comprise 11 geo-referenced images (corresponding to the bands described in table 1), and two . Reproduced from [17]. Image stated to be in the public domain.
plain-text files of metadata. In this work we use images from bands 4, 5, 7 and 10, taken on August 14, 2019 and June 24, 2018.

Image preprocessing
The size of Landsat 8 images is 7861 × 7981 pixels with a 16 bit pixel depth; each pixel describes an area of 30 × 30 m 2 . They are delivered as 125 MB TIFF files. This large amount of information cannot be easily handled and for this reason, we decided to reduce the number of pixels by a factor of 4 (i.e. 3930 × 3990). Figure 1 shows Landsat images recorded on August 14, 2019. Moreover, the original images are geometrically manipulated to remove pixels with values equal to zero, that can distort the data range during the processing stage. However, this procedure introduces some resampling artifacts: in particular, the images display a slightly but noticeable shift. This effect can be compensated, as explained in section 5. The use of this approach is interesting from a pedagogical point of view, because important image processing techniques such as affine transform and 2D-correlation can be introduced. However, it is not required to modify the images since the zero values can be masked. Homography matrix H relates the position of the pixels of the image before and after an arbitrary geometrical transformation.
where (x, y, 1) and (x , y , 1) refer to the pixel coordinates before and after the transformation [16]. Note that the third component (with value equal to 1) is used for mathematical convenience. Usually, matrix H is normalized in such a way that h 33 = 1.
In what follows, the size of the image is N (rows) × M (columns); the coordinates of the top-left and the bottom-right corners are (0, 0) and (N − 1, M − 1) respectively (see figure 2). Then, the coordinates (r, c) of the four vertices v 1 , v 2 , v 3 and v 4 , are: • v 1 with coordinates (0, c 1 ): column c 1 is determined as the column of the first pixel with non-zero value in the first row of the image. Then, the lengths of each side of the image can be calculated as where . stands for the norm. Note that in the present case, l 1 ≈ l 3 and l 2 ≈ l 4 and thus, parallelism is preserved. Accordingly, the geometrical problem (equation (1)) becomes an affine . Adapted from [17]. Image stated to be in the public domain.
Then, vertices v 1 , v 2 , v 3 , and v 4 will be transformed to a new set of vertices p 1 , p 2 , p 3 , p 4 with coordinates: The resulting images are displayed in figure 3. The overlaid rectangle indicates the area of interest where the fire burned. The physical length of this area is 36 × 36 km 2 .

Examination of thermal imagery
The top of atmosphere spectral radiance (L λ , measured in W m −2 srad −1 μm −1 ) can be calculated using the image recorded by the LWIR1 sensor (band 10), namely I 10 . L λ and I 10 are related by means of where M L = 3.3420 × 10 −4 and A L = 0.1 are band-specific factors that can be found in the metadata file attached to the image set [18]. Note that I 10 is a 16 bit image with values ranging from 0 to 65 535. Now, the spectral radiance is converted to the at-sensor equivalent blackbody temperature T s using a model based on Planck's law: where k B , h and c are the Boltzmann's and Planck's constants, and the speed of light, respectively. T s is determined using effective thermal constants K 1 = 774.8853 and K 2 = 1321.0789 that matches the behavior of the satellite sensors. This information can be found in the metadata file that accompanies the Landsat images: At-sensor temperature, T s . Adapted from [17]. Image stated to be in the public domain.
Note that the at-sensor temperature T s of the burned area is close to 320 K, as shown in figure 4.

Normalized burned ratio (NBR)
The normalized burn ratio (NBR) index is introduced to assess the severity of the damage produced by fire [7]. NBR is defined as: where I 5 and I 7 are Landsat 8 images recorded on bands 5 (NIR) and 7 (SWIR) respectively. Healthy vegetation shows a very high reflectance in the near infrared NIR band and low reflectance in the SWIR. But, recently burned areas show low NIR reflectance and high SWIR reflectance. The difference between the spectral responses of healthy vegetation and burned areas reach a peak in the NIR and the SWIR regions of the spectrum. The differential NBR (dNBR) compares data obtained before and after the wildfire. It is used to map the level of severity of the burned areas [7]: dNBR = NBR pre−fire − NBR post−fire . The severity of a burned area can be assessed according to table 2 [19]: Figures 5(a)-(d) display the calculated NBR distributions. Figures 5(a)-(c) show the NBR for images recorded on 2018-06-24 (pre-fire) and 2019-08-14 (post-fire), and the differential dBNR respectively. Note that the images are slightly misaligned due to the affine transform resampling used to remove those pixels with zero value, and this can lead to some errors in the dNBR calculation. The relative shift among these two image can be detected by calculating the 2D cross-correlation among them. The position of the cross-correlation maximum with respect to the center of the 2D array determines the image shift. Given two arbitrary signals f and g,  . The color scale applies only to cases (c) and (d). Adapted from [17]. Image stated to be in the public domain.
the cross-correlation product is defined as where the symbol * stands for complex conjugate. Taking advantage of the convolution theorem, cross-correlation can be calculated in a computational efficient way as: where F = F{ f } and G = F{g} are the corresponding Fourier transforms of f and g, respectively. Figure 5(d) displays the new version of the dNBR after removing the relative shift among the NBR images (figures 5(a) and (b)). Interestingly, figure 5(d) is much less noisy that figure 5(c).

Wildfire surface estimation
To provide an estimation of the burned area, the dNBR distribution is binarized, using a threshold value equal to 0.100 (low severity damage lower bound, see table 2). The resulting image is shown in figure 6(a). It is noticeable that numerous pixels that do not belong to the burned area exceed the threshold. These pixels can be removed with the help of mathematical morphology. Let B be a binary image and E is the so-called structuring element. The pixels of B are processed one at a time. E is superimposed on one of the pixels of B in such a way that the center of array E determines the pixel of B being processed. Erosion and dilation operators are described as follows: Figure 6. (a) Binary dNBR; (b) binary dNBR after opening; (c) dNBR after opening and small objects removal; (d) dNBR after opening, small objects removal and closing. Adapted from [17]. Image stated to be in the public domain.
(a) The resulting pixel after erosion is 1 if and only if all the pixels of B that are covered by E are 1; otherwise, the pixel is set to 0. (b) The resulting pixel after dilation is 1 if any of the pixels of B that are covered by E are 1; otherwise, the pixel becomes 0.
Opening and closing are combinations of erosion and dilation: opening is the dilation of the erosion; it is used to remove small bright objects. Closing is the erosion of the dilation, and it is useful to remove small holes. A detailed explanation on morphological techniques can be found in [16].
In the present problem, we first used a binary opening operator with a 3 × 3 structuring array E to discard the unconnected pixels ( figure 6(b)). Then, the remaining isolated spots are deleted using the remove small objects function (figure 6(c)). This operator looks for isolated items smaller than a predetermined number of pixels (64 by default). Finally, we use a binary closing operator with a 5 × 5 structuring array E, to fill some minor gaps inside the area of interest [16].
Then, the estimation of the area burned is just the number of white pixels times the corresponding area of one pixel. We assess a burned area of 54.4 km 2 , slightly bigger than the official figure ( 50 km 2 ) [20].

Conclusions
In this paper we described a method to estimate the damage caused by a wildfire. Using Landsat 8 multi-spectral images recorded several days after a wildfire, we implemented a procedure to determine the temperature of the burned area. This calculation involved the combined use of infrared and visible images. Then, the differential NBR is used to compare infrared images recorded before and after the wildfire. Using this metric and with the help of mathematical morphology techniques, we provided a reasonable estimation of the area burned.