Impacts of the Tropical Cyclone Idai in Mozambique: A Multi-Temporal Landsat Satellite Imagery Analysis

: The Central Region of Mozambique (Sofala Province) bordering on the active cyclone area of the southwestern Indian Ocean has been particularly affected by climate hazards. The Cyclone Idai, which hit the region in March 2019 with strong winds causing extensive ﬂooding and a massive loss of life, was the strongest recorded tropical cyclone in the Southern Hemisphere. The aim of this study was to use pre- and post-cyclone Idai Landsat satellite images to analyze temporal changes in Land Use and Land Cover (LULC) across the Sofala Province. Speciﬁcally, we aimed—(i) to quantify and map the changes in LULC between 2012 and 2019; (ii) to investigate the correlation between the distance to Idai’s trajectory and the degree of vegetation damage, and (iii) to determine the damage caused by Idai on different LULC. We used Landsat 7 and 8 images (with 30 m resolution) taken during the month of April for the 8-year period. The April Average Normalized Difference Vegetation Index (NDVI) over the aforementioned period (2012–2018, pre-cyclone) was compared with the values of April 2019 (post-cyclone). The results showed a decreasing trend of the productivity (NDVI 0.5 to 0.8) and an abrupt decrease after the cyclone. The most devastated land use classes were dense vegetation (decreased by 59%), followed by wetland vegetation ( − 57%) and shrub land ( − 56%). The least damaged areas were barren land ( − 23%), barren vegetation ( − 27%), and grassland and dambos ( − 27%). The Northeastern, Central and Southern regions of Sofala were the most devastated areas. The Pearson Correlation Coefﬁcient between the relative vegetation change activity after Idai (NDVI%) and the distance to Idai’s trajectory was 0.95 (R-square 0.91), suggesting a strong positive linear correlation. Our study also indicated that the LULC type (vegetation physiognomy) might have inﬂuenced the degree of LULC damage. This study provides new insights for the management and conservation of natural habitats threatened by climate hazards and human factors and might accelerate ongoing recovery processes in the Sofala Province.


Introduction
Tropical cyclones are among the most devastating natural disasters owing to their great potential for loss of human life, significant economic decline and severe environmental damage [1][2][3]. The Southwestern Indian Ocean is one of the main tropical cyclone areas 2012 to 2019; (ii) to investigate how the distance to Idai's trajectory related to vegetation damage, and (iii) to determine the degree of damage caused by the cyclone on different LULC classes. The findings of this study are meant to assist managers of natural resources to design and implement efficient strategies and practices in order to safeguard natural ecosystems and their services such as provisioning (e.g., food, and timber), regulating (e.g., climate regulation, and water purification), supporting (e.g., nutrient cycling, and soil formation), and cultural (recreation, and spiritual) [30].

Study Area
Sofala is a coastal province in Central Mozambique, occupying a surface area of approximately 68,018 km 2 (about 8.5% of the country). It borders on the Indian Ocean to the east, on the Zambezi River to the north, and the Save River to the south (Figure 1a). The Central part of Sofala is intersected by the Púnguè and the Buzi River, and the province is characterized by an inter-connected hydrological system of creeks, swamps, and lakes. All the region's main rivers discharge into the Indian Ocean, providing a suitable habitat for the growth and establishment of a great variety animal and plant species. According to the 2017 census, the total population of Sofala province was estimated at 2,259,248 inhabitants, most of them living in rural areas [31] and largely dependent on unstable natural resources (threatened by climate hazards) for their subsistence. The main economic activities of the local population are slash and burn agriculture, fishing, raising cattle and commerce. The Province of Sofala is characterized by a tropical climate with rainy season (summer) running from November to March, and a dry season (winter) from April to October; southeasterly trade winds are predominant. The annual average temperature is 25 °C and average rainfall amounts to approximately 1300 mm/year [32]. According to Marzoli [33], LULC classes in Sofala are primarily composed of vast native forests (48.81%), wetlands (19%), agricultural lands (5.57%), urban areas (0.16%), barren areas (0.5%), with water and other vegetation formations (height < 5 m) accounting for 25.96%. Sofala's coastline is characterized by the so-called swamp coast [34] and it has one of the largest mangrove areas (932 km 2 ) in Mozambique, second only to the Zambezia province (1219 km 2 ) [33]. Sofala is the province most prone to floods and cyclones (occurring in summer) which greatly affect LULC [9,11].    (Figure 1b). It brought strong winds (180-220 km/h) and torrential rain (more than 200 mm in 24 h), flooding more than 3000 km 2 of agricultural land. It displaced around 400,000 inhabitants, provoking over 600 deaths, while severely damaging LULC classes [12]. The winds receded as the cyclone moved inland, affecting provinces bordering on Sofala (Manica, Zambézia, Tete, and Inhambane) before moving to neighboring Zimbabwe.

Data
We used Landsat 7 Enhanced Thematic Mapper Plus (ETM+) and Landsat 8 Operational Land Imager (OLI) images from April of each year, from 2012 to 2015, and 2016 to 2019, respectively. Landsat 7 ETM+ and Landsat 8 OLI data both had a 30 m resolution and were downloaded from the United States Geological Survey (USGS) gateway (https://earthexplorer.usgs.gov/). Table 1 shows the list of satellite images used, acquisition date, and path row for the study area. The LULC and NDVI analyses can be influenced by seasonal differences (because of changes in atmospheric conditions and photosynthetic activities) [35].

Methods
This study was carried out in four main steps-(1) satellite data pre-processing, which included cloud detection, atmospheric correction, and geometric correction; (2) calculation of the vegetation indicators, including the Normalized Difference Vegetation Index (NDVI), the difference in NDVI (∆NDVI), and the relative change in vegetation productivity after Idai (NDVI%); (3) producing the LULC map for April 2019; and (4) producing a distance bands map obtained using buffer analysis. These steps are detailed in Figure 2. The distance to the cyclone's trajectory was subjected to statistical analyses to evaluate its relationship with vegetation damage using a linear regression model with a single explanatory variable.

Methods
This study was carried out in four main steps-(1) satellite data pre-processing, which included cloud detection, atmospheric correction, and geometric correction; (2) calculation of the vegetation indicators, including the Normalized Difference Vegetation Index (NDVI), the difference in NDVI (∆NDVI), and the relative change in vegetation productivity after Idai (NDVI%); (3) producing the LULC map for April 2019; and (4) producing a distance bands map obtained using buffer analysis. These steps are detailed in Figure 2. The distance to the cyclone's trajectory was subjected to statistical analyses to evaluate its relationship with vegetation damage using a linear regression model with a single explanatory variable. Initially, the Landsat 7 and 8 satellite data were geometrically corrected and orthorectified ( Figure 2) using the Landsat packages ("georef" and "geoshift") in RStudio Software [36][37][38]. All satellite images were georeferenced through the Universal Transverse Mercator (UTM) coordinate system. We were aware of the limitation of Landsat 7 ETM+ because of the Scan Line Corrector (SLC) failure [39,40]. To fill the gaps (Not available "NA" values) caused by that error, we applied the Landsat 7 Scan Line Corrector (SLC)off Gap function, using the traditional Local Linear Histogram Matching method-LLHM [37,41]. The ERDAS image processing software (version 8.7) was applied to perform Scan line error correction with the help of band-specific gap mask files, made with the Landsat 7 Level-1 data product. These mask files help to classify the location of every pixel affected by the original data gaps in the primary SLC-off scene [36]. For this purpose, the DV values (0-255) were converted into top-of-atmosphere (TOA) radiance (at-satellite radiance values) using the parameters provided in metadata files [38,42,43]. We also applied absolute atmospheric compensation techniques (Dark Object and Modified Dark Object Subtraction Method) to identify and remove clouds, aerosols, and cirri [36,44]. Initially, the Landsat 7 and 8 satellite data were geometrically corrected and orthorectified ( Figure 2) using the Landsat packages ("georef" and "geoshift") in RStudio Software [36][37][38]. All satellite images were georeferenced through the Universal Transverse Mercator (UTM) coordinate system. We were aware of the limitation of Landsat 7 ETM+ because of the Scan Line Corrector (SLC) failure [39,40]. To fill the gaps (Not available "NA" values) caused by that error, we applied the Landsat 7 Scan Line Corrector (SLC)-off Gap function, using the traditional Local Linear Histogram Matching method-LLHM [37,41]. The ERDAS image processing software (version 8.7) was applied to perform Scan line error correction with the help of band-specific gap mask files, made with the Landsat 7 Level-1 data product. These mask files help to classify the location of every pixel affected by the original data gaps in the primary SLC-off scene [36]. For this purpose, the DV values (0-255) were converted into top-of-atmosphere (TOA) radiance (at-satellite radiance values) using the parameters provided in metadata files [38,42,43]. We also applied absolute atmospheric compensation techniques (Dark Object and Modified Dark Object Subtraction Method) to identify and remove clouds, aerosols, and cirri [36,44].
The pseudo-invariant features (PIF) function was used to examine the homogeneity of reflectance values of LULC in the images, and then corrected by using major axis regression. The Landsat 7 and 8 data were radiometrically and atmospherically corrected by using an atmospheric simulation model available in the Landsat and R package (RStoolBox) [38,45].

Land Use and Land Cover Classification
The satellite images for 2019 were classified into nine LULC classes-dense vegetation, shrub land, grassland and dambos, agricultural land, wetland vegetation, barren vegetation, barren land, built-up areas, and water bodies ( Table 2). This was performed by applying the Random Forest (RF) classifier. This technique has been recently used for mapping plant species and landscape due to its processing speed, improved accuracy, and reliable classification outputs [46,47]. The processing chain for the RF classification algorithm optimized the proximities between data points [48]: For the LULC classification, we used "randomForest" packages [49] in the open-source RStudio software version 1.3.1073 [50]. Furthermore, the packages Classes and Methods for Spatial Data (Sp) [51], raster [52], and Raster Geospatial Data Abstraction Library (Rgdal) [53] were used to process and visualize the spatial data. Waterbodies River, open water, lakes, streams, estuaries and ponds

Accuracy Assessment
The accuracy assessment is a widely applied technique to quantify how close the result of land cover classification is to the reference image. In this study, the accuracy was assessed by comparing the results of the obtained LULC map (classified image) with reference Google Earth images retrieved from the Google Earth Engine (GEE) [54] for the year 2019. From the classified LULC, seventy-five random point samples were produced for the study period. They were visually distinguished from GEE by comparison with the classified map. A confusion matrix was constructed to assess accuracy [48]. The quantitative accuracy assessment of the classified map was obtained by calculating the kappa coefficient (k) using ERDAS Imagine (version 8.7) [55]. The model's accuracy is classified according to k values as follows-poor (k < 0.2), moderate (0.4 < k < 0.6) and near perfect (k > 0.8) (see Jensen [35] and Louarn et al. [47], for further detail).

Vegetation Indices
The Normalized Difference Vegetation Index (NDVI) is generally used to measure vegetation density and its health status (level of photosynthetic activity) and is less affected by topographic factors and illumination [24]. The NDVI is calculated as follows: where NIR and RED are the spectral reflectances corresponding to the fourth (0.77-0.90 µm) and third (0.63-0.69 µm) Landsat ETM+ bands, respectively [24]. For Landsat 8 OLI, NIR is the fifth band (0.85-0.88 µm) and RED is the fourth band (0.64-0.67 µm) [23]. Normally, NDVI ranges between −1.0 and 1.0 with vegetation land covers ranging from 0.0 to 1.0. The difference in NDVI (hereafter referred to as ∆NDVI) can show the change in LULC, while a negative ∆NDVI represents the vegetation damage caused by Cyclone Idai. It is calculated by subtracting the NDVI image of one date (after) from that of another (before) using map algebra, which is a cell-by-cell process [24,56]. In this study, the Average NDVI in April over the 7-year period (2012-2018) represents the situation before Cyclone Idai (hereafter referred as pre-cyclone) whereas the NDVI of April 2019 represents the postcyclone situation (hereafter referred as post-cyclone). The ∆NDVI is calculated using the following equation: where NDVIpre and NDVIpost are NDVI before and after cyclone, respectively. We also calculated the relative vegetation change activity after Idai (NDVI%) by using the equation below: Please note that a higher NDVI% indicates lower damage. We used Sp [51], raster [52], rts [57], and Rgdal [53] packages in R studio for calculation of NDVI and the changes in vegetation productivity.

Distance to the Cyclone Trajectory
A distance map based on the Cyclone Idai's track was produced with distance bands calculated using ArcGIS (ESRI, 2020); 10 multiple buffers at the specified distance of 25 km from the best track line were obtained [23].

Correlation Analysis
We used the Pearson Correlation Coefficient to assess the relationship between vegetation damage and the distance to the cyclone's trajectory. This relationship was measured by means of a least squares estimator from the linear regression model with a single influencing/explanatory variable [23,24].  (Table 3) shows a high accuracy level, indicating a good agreement between the LULC map for April 2019 and ground-truth based on GEE images [35,47]. Table 3 shows that the Producer's accuracy (i.e., how often are real features on the ground correctly shown on the LULC map) and user's accuracy (i.e., how often the class on the LULC map will actually be present on the ground) presented different values for each class, with a much higher accuracy for barren vegetation (PA: 92.5; UA: 94.1) [58].     Figure 4 illustrates the spatial distribution of pre-cyclone (left) and post-cyclone (right) NDVI. The comparison by optical remote sensing shows the loss of green leaves on the postcyclone imagery, as indicated by the increase of no productivity and/or low productivity pixels [23,24]. Significant NDVI alterations were observed across the entire Province, with the largest LULC changes detected in the Northeast (NE), Central (C), and Southern (S) regions.

NDVI
Remote Sens. 2021, 13, x FOR PEER REVIEW 9 of 17 Figure 4 illustrates the spatial distribution of pre-cyclone (left) and post-cyclone (right) NDVI. The comparison by optical remote sensing shows the loss of green leaves on the post-cyclone imagery, as indicated by the increase of no productivity and/or low productivity pixels [23,24]. Significant NDVI alterations were observed across the entire Province, with the largest LULC changes detected in the Northeast (NE), Central (C), and Southern (S) regions.  Table 4 shows the land cover classes affected by the Cyclone Idai. All post-cyclone images were acquired within three or four weeks of Idai's passage through the Sofala Province, i.e., before the damaged vegetation had recovered. The ∆NDVI ranged from    Table 4 shows the land cover classes affected by the Cyclone Idai. All post-cyclone images were acquired within three or four weeks of Idai's passage through the Sofala Province, i.e., before the damaged vegetation had recovered. The ∆NDVI ranged from −0.07 to −0.46; the highest damage was found for dense vegetation (−58.9%), wetland vegetation (−57.4%), and shrub land (−55.5%); the least damage was observed in barren land (−22.5%), barren vegetation (−26.9%), and grassland and dambos (−27.1%).  Figure 5 clearly shows an abrupt decrease of the area of highly productive vegetation classes (0.5-0.8) in April 2019, associated with a substantial growth of the low productivity vegetation area. Overall, almost all the highly productive vegetation experienced a severe decrease (over 99%) while the areas with low productivity vegetation increased by 22-80% as shown in the inset table of Figure 5.

NDVI
Remote Sens. 2021, 13, x FOR PEER REVIEW 10 of 17 Figure 5 clearly shows an abrupt decrease of the area of highly productive vegetation classes (0.5-0.8) in April 2019, associated with a substantial growth of the low productivity vegetation area. Overall, almost all the highly productive vegetation experienced a severe decrease (over 99%) while the areas with low productivity vegetation increased by 22-80% as shown in the inset table of Figure 5.

Influence of the Distance to the Cyclone Trajectory
The distance to the Cyclone Idai trajectory map is shown in Figure 6 for Sofala Province. The NDVI% was calculated for each distance class (Table 5). NDVI% increased with the distance from the cyclone trajectory (Figure 7 and Table 5).
At distances shorter than 75 km from the cyclone trajectory, land degradation was higher than 50% in NDVI. Even at distances over 200 km, NDVI was still above 20% ( Figure 7 and Table 5).

Influence of the Distance to the Cyclone Trajectory
The distance to the Cyclone Idai trajectory map is shown in Figure 6 for Sofala Province. The NDVI% was calculated for each distance class (Table 5). NDVI% increased with the distance from the cyclone trajectory (Figure 7 and Table 5).    At distances shorter than 75 km from the cyclone trajectory, land degradation was higher than 50% in NDVI. Even at distances over 200 km, NDVI was still above 20% ( Figure 7 and Table 5).

Figure 7.
Pearson's correlation between mean NDVI% and distance to the Idai's trajectory.

The Changes in NDVI
The vegetation productivity changes from 2012 to 2019 were quantified ( Figure 5). In general, the areas with no productivity and/or low productivity tended to increase with time and an abrupt increase was observed in April 2019 (post-cyclone), whereas the high productivity pixels showed an opposite trend. The loss of vegetation productivity with time is a great concern globally and, most particularly, in low-income countries like Mozambique prone to climatic hazards. Mozambican forests and other natural ecosystems are under greater anthropogenic pressure (e.g., slash and burn agriculture, timber for building local infrastructures, firewood, and charcoal production, and urban expansion) while being adversely affected by natural factors (e.g., climate hazards, rising sea levels and frequent flooding) [11,59,60]. The category four tropical Cyclone Idai was the most devastating and deadliest ever recorded in the Southern Hemisphere. It brought strong winds and torrential rains, causing extensive flooding and leaving many communities in the Sofala Province submerged under 10 m of water [61]. It also caused a loss of 99.5% of high-productivity vegetation and an increment of up to 80.5% of the low-and moderate-productivity areas ( Figure 5).
NDVI clearly changed over the whole Sofala Province after Idai's passage. The Northeastern, Central, and Southern regions (Figure 4) appear to be the most devastated areas. At "NE" a high concentration of the pixels of no productivity and/or low productivity region is clearly visible, where our results show the predominance of dense vegetation, agricultural area and water courses with regularly flooded wood and/or herbaceous species (Figures 3 and 4) [62]. Our findings are in line with Couto et al. [63] and Ramsar [64] who described this Northeastern region as part of the Ramsar site of the Zambezi Delta and of the Marromeu Complex, which are designated Wetlands of International Importance (conservation hotspot). This region is characterized by a broad flat alluvial area with a rich mosaic of grassland and dambos, woodland, thicket, large water swamps, and an extensive area covered in mangrove. The devastation in Central Sofala Province is

The Changes in NDVI
The vegetation productivity changes from 2012 to 2019 were quantified ( Figure 5). In general, the areas with no productivity and/or low productivity tended to increase with time and an abrupt increase was observed in April 2019 (post-cyclone), whereas the high productivity pixels showed an opposite trend. The loss of vegetation productivity with time is a great concern globally and, most particularly, in low-income countries like Mozambique prone to climatic hazards. Mozambican forests and other natural ecosystems are under greater anthropogenic pressure (e.g., slash and burn agriculture, timber for building local infrastructures, firewood, and charcoal production, and urban expansion) while being adversely affected by natural factors (e.g., climate hazards, rising sea levels and frequent flooding) [11,59,60]. The category four tropical Cyclone Idai was the most devastating and deadliest ever recorded in the Southern Hemisphere. It brought strong winds and torrential rains, causing extensive flooding and leaving many communities in the Sofala Province submerged under 10 m of water [61]. It also caused a loss of 99.5% of high-productivity vegetation and an increment of up to 80.5% of the low-and moderate-productivity areas ( Figure 5).
NDVI clearly changed over the whole Sofala Province after Idai's passage. The Northeastern, Central, and Southern regions (Figure 4) appear to be the most devastated areas. At "NE" a high concentration of the pixels of no productivity and/or low productivity region is clearly visible, where our results show the predominance of dense vegetation, agricultural area and water courses with regularly flooded wood and/or herbaceous species (Figures 3 and 4) [62]. Our findings are in line with Couto et al. [63] and Ramsar [64] who described this Northeastern region as part of the Ramsar site of the Zambezi Delta and of the Marromeu Complex, which are designated Wetlands of International Importance (conservation hotspot). This region is characterized by a broad flat alluvial area with a rich mosaic of grassland and dambos, woodland, thicket, large water swamps, and an extensive area covered in mangrove. The devastation in Central Sofala Province is closely related to Idai's trajectory and the influence of the Púnguè and Buzi rivers and Urema Lake whose basins are locally considered as wetlands [65,66]. In this region we found a mixed spatial distribution of the land cover, which mainly includes dense vegetation, agricultural land, and wetland along the coastline. The Southern region is characterized by abundant shrub land and a mixed vegetation distribution pattern including wetland vegetation, dense vegetation, and agriculture land which were all strongly affected by Cyclone Idai. To the south, the Save River, constituting the southern border of the Sofala Province, forms an estuary that is also considered wetland [63]. Li et al. [62] reported particularly significant levels of land degradation in the south-eastern region of the province, which could explain the more reduced agricultural land coverage we observed in this region.

Effects of Idai on Different LULC
The changes in the different land-cover classes were calculated and compared with published literature. Dense vegetation (decreased by 58.9%), wetland vegetation (−57.4%) and shrub land (−55.5%) were among the most devastated vegetation classes, while barren land (−22.5%), barren vegetation (−26.9%) and grassland and dambos (−27.1%) were the least affected. From a biological perspective, our results show that the intensity of damage to different vegetation classes is closely related to their physiognomy. Dense vegetation, wetland vegetation, and shrub land are more vulnerable to physical damage by a tropical cyclone than the barren areas and herbaceous or sparse vegetation because of the higher tree layer and greater tree canopy cover. Numerous studies reported that forest stands with older and bigger trees (height, diameter and canopy) are more prone to wind damage than open, short and sparse forest stands [67][68][69]. An intense tropical cyclone like Idai typically causes great defoliation, branch stripping, bole snapping and uprooting [70][71][72]. Other individual tree characteristics, such as less dense wood, poor health status, and growth in more exposed sites (e.g., coastal Sofala, as shown by Cabral et al. [11] and Charrua et al. [10]), contribute to increasing its vulnerability to cyclone-related damage [70,73]. Other factors, such as the size of each LULC area, need to be deeply investigated. Rossi et al. [74] reported greater damage in larger forests in Northern Nicaragua as a result of hurricane Felix (2007). Generally, the spatial distribution of severe post-Cyclone Idai damage in Sofala closely matched the spatial distribution of the most devastated LULC. The damage behavior shown in different LULC in our study is consistent with other findings [23,24,70,73].

The Influence of Distance on Damage
Our results demonstrate a strong positive correlation between distance and damage. The greatest damage in NDVI was found at distances shorter than 75 km (NDVI% from −52.29 to −55.07, Table 5), in accordance with previous studies [16,23,75]. The Pearson Correlation Coefficient between NDVI% and distance was 0.95 and R-square was 0.91, indicating a strong positive linear correlation (Figure 7). The spatial distribution of highly damaged regions across the Sofala Province ( Figure 4) and the distance to Idai's trajectory suggest that the degree of damage does not entirely depend on the latter and that other factors may interfere, like the LULC type (vegetation physiognomy) [74,76].
Although this was not a ground data-based study because of the well-known difficulties of field research (e.g., time-consuming, limited availability of human and logistic resources), our results are realistic and they were satisfactorily verified with the local damage reports [12]. In fact, Landsat image analyses provided reliable results, and have been widely and successfully used to map LULC damages from cyclones [16,24,27,77]. This study can be employed as a future reference in related studies. Our approach clearly identified changes between pre-and post-cyclone LULC in Sofala; moreover, we found a strong correlation between the degree of damage and the distance from Idai's trajectory. Nevertheless, further studies are needed to quantify the local/regional changes that did not fit the linear model.

Conclusions
This study focuses on an analysis of the damage caused to the LULC by Cyclone Idai, which hit Sofala Province in Mozambique in March 2019. For that purpose, we used Landsat 7 and Landsat 8 images taken during the month of April over an 8-year period (2012-2019). All of the LULC classified areas showed a decrease of NDVI after the cyclone. The greatest damage was found in dense vegetation, wetland vegetation, and shrub land; barren land, barren vegetation, and grassland and dambos showed the lowest relative damage. The most heavily affected regions were the Northeast, Central and Southern Sofala. The distance to Idai's trajectory greatly influenced LULC damage levels-the greater the distance, the lower the damage. Besides the distance to Idai's trajectory, other factors such as LULC type (vegetation physiognomy) may have played an important role in terms of vegetation damage. The information here provided is relevant for LULC managers and all stakeholders to take appropriate measures for better planning and future management of the territory. In addition, our findings may help to speed up ongoing recovery processes that were activated in the wake of tropical Cyclone Idai.

Data Availability Statement:
Publicly available datasets were analyzed in this study. This data can be found here: https://earthexplorer.usgs.gov.