Refugee Camp Monitoring and Environmental Change Assessment of Kutupalong, Bangladesh, Based on Radar Imagery of Sentinel-1 and ALOS-2

: Approximately one million refugees of the Rohingya minority population in Myanmar crossed the border to Bangladesh on 25 August 2017, seeking shelter from systematic oppression and persecution. This led to a dramatic expansion of the Kutupalong refugee camp within a couple of months and a decrease of vegetation in the surrounding forests. As many humanitarian organizations demand frameworks for camp monitoring and environmental impact analysis, this study suggests a workﬂow based on spaceborne radar imagery to measure the expansion of settlements and the decrease of forests. Eleven image pairs of Sentinel-1 and ALOS-2, as well as a digital elevation model, were used for a supervised land cover classiﬁcation. These were trained on automatically-derived reference areas retrieved from multispectral images to reduce required user input and increase transferability. Results show an overall decrease of vegetation of 1500 hectares, of which 20% were used to expand the camp and 80% were deforested, which matches ﬁndings from other studies of this case. The time-series analysis reduced the impact of seasonal variations on the results, and accuracies between 88% and 95% were achieved. The most important input variables for the classiﬁcation were vegetation indices based on synthetic aperture radar (SAR) backscatter intensity, but topographic parameters also played a role.


Introduction
Remote sensing has become increasingly popular in the field of humanitarian action because it is an independent and reliable source of information and allows both a quick response to emergencies and the monitoring of gradual changes [1]. This is of particular importance when observations in the field area are not possible because of limited budget, legal restrictions, or security reasons [2][3][4]. Observation of specific places from space is not only crucial for decision making involving cases of humanitarian response to natural disasters or emergencies [5,6], but it also helps to develop a better general understanding of an area and how trends and temporal dynamics have shaped the spatial patterns of the present [7,8]. This temporal aspect is often more valuable than the level of spatial detail of the images alone, as it puts the information into a context and increases the interpretability of the retrieved information [9].
Institutionalized in the International Charter on Space and Major Disasters, a global collaboration of space agencies and disaster response units provides satellite images and value-added products to relief organizations to assist their operations [10]. However, other institutions and platforms also make operational use of satellite images for humanitarian assistance, such as UNOSAT (United Nations To increase the potential of the radar images' information, missions of Sentinel-1 (S1, ESA) and ALOS-2 (A2, JAXA) are combined in this study. This allows the measuring of backscatter intensity in all four polarizations (horizontal-horizontal (HH), horizontal-vertical (HV), vertical-horizontal (VH) vertical-vertical (VV)). Moreover, the different wavelengths of the sensors can be exploited more effectively-the C-band of Sentinel-1 operates at a comparably short wavelength (5.5 cm) which is strongly sensitive to the roughness of canopies and moisture of un-vegetated surfaces [65]. In turn, the large wavelength of ALOS-2 (22.9 cm) is able to penetrate surfaces to a certain extent and can therefore be used to describe structures under canopies, such as stems or branches [66]. This was already confirmed by Inoue et al., who report that the C-band is a more suitable canopy and grassland descriptor, while the L-band has more informative value regarding the biomass of vegetation [67]. Consequently, the combination of both sensors allows for the ability to analyze backscatter mechanisms of both surfaces and covered structures. This complementary use was already successfully demonstrated for biomass estimation [68,69], deforestation monitoring [70,71], the mapping of flooded vegetation [72], and rice or palm oil cultivation [73,74]. However, it has to be denoted that, while Sentinel-1 data can be freely accessed within the Copernicus Programme, ALOS-2 data is commercially distributed by the Japan Aerospace Exploration Agency (JAXA). This aspect will be discussed in section 3.3.
As indicated in the studies presented in the previous section, the study area shows high phenological variations over the year. To be able to take these into consideration and prevent the deduction of false conclusions, SAR image pairs of S1 and A2 for eleven dates between 2016 and March 2019 (Table 1) were selected. As shown in Table 1, most of the pairs had a temporal baseline of three days or less. Furthermore, only images with the same flight direction (ascending or descending orbit) were selected to minimize the impact of topography within one investigated date. To increase the potential of the radar images' information, missions of Sentinel-1 (S1, ESA) and ALOS-2 (A2, JAXA) are combined in this study. This allows the measuring of backscatter intensity in all four polarizations (horizontal-horizontal (HH), horizontal-vertical (HV), vertical-horizontal (VH) vertical-vertical (VV)). Moreover, the different wavelengths of the sensors can be exploited more effectively-the C-band of Sentinel-1 operates at a comparably short wavelength (5.5 cm) which is strongly sensitive to the roughness of canopies and moisture of un-vegetated surfaces [65]. In turn, the large wavelength of ALOS-2 (22.9 cm) is able to penetrate surfaces to a certain extent and can therefore be used to describe structures under canopies, such as stems or branches [66]. This was already confirmed by Inoue et al., who report that the C-band is a more suitable canopy and grassland descriptor, while the L-band has more informative value regarding the biomass of vegetation [67]. Consequently, the combination of both sensors allows for the ability to analyze backscatter mechanisms of both surfaces and covered structures. This complementary use was already successfully demonstrated for biomass estimation [68,69], deforestation monitoring [70,71], the mapping of flooded vegetation [72], and rice or palm oil cultivation [73,74]. However, it has to be denoted that, while Sentinel-1 data can be freely accessed within the Copernicus Programme, ALOS-2 data is commercially distributed by the Japan Aerospace Exploration Agency (JAXA). This aspect will be discussed in Section 3.3.
As indicated in the studies presented in the previous section, the study area shows high phenological variations over the year. To be able to take these into consideration and prevent Remote Sens. 2019, 11,2047 6 of 34 the deduction of false conclusions, SAR image pairs of S1 and A2 for eleven dates between 2016 and March 2019 (Table 1) were selected. As shown in Table 1, most of the pairs had a temporal baseline of three days or less. Furthermore, only images with the same flight direction (ascending or descending orbit) were selected to minimize the impact of topography within one investigated date. To quantify the decrease of forest areas during the investigated period, a supervised classification was chosen (see Section 2.2). However, this required the input of training areas whose collection was either not possible in the field or biased by visual interpretation of existing images (see challenges outlined in Section 1). For this reason, optical data of Landsat ETM+ (L7), Landsat OLI (L8), and Sentinel-2 (S2) were selected on the basis of their acquisition date and amount of cloud cover ( Table 2). Even if some images had gaps because of the L7 SLC-off incident since 2003 [75], were partly covered by clouds, or did not fully cover the study area, they were still usable to retrieve training areas as described in the following section.

Overall Study Design
The overall workflow of the study is displayed in Figure 2. The single components will be explained in more detail in the subsequent sections, but the general approach is shortly outlined as follows. Spectral indices were derived on the optical images to serve as an input for a semi-automated pre-classification of the four target classes for each date (Section 2.2.2). The following classes were selected for the study area: (1) Vegetation-around 40% of the study area is covered by vegetation, mainly consisting of mixed tropical evergreen and semi-evergreen forest dominated by tall Garjan (Dipterocarpus spp.) trees [76]. Peaks of the hills are predominantly covered by sun grass (Imperata cylindrica), herbs, shrubs, and brush woods [55], but their share in the study area is low because these elevations are only reached further in the north and east (as indicated in Figure 1). (2) Open lands (abbreviated as "open" in the following sections)-includes all areas outside the forests which are only sparsely covered by grassy plants, barren lands, or degraded forests affected by shifting cultivation as identified by Giri and Shrestha [77]. (3) Water-the course and floodplains of the Naf River in the southeast of the study area. This also includes the frequently waterlogged marshes in this valley. (4) Informal settlements (abbreviated as "camp" in the following sections)-the study area only contains the refugee camp of Kutupalong and its extension areas. This class represents all areas that are covered by solid humanitarian facilities and simple dwellings of lighter construction materials.
Remote Sens. 2019, 11, x FOR PEER REVIEW 7 of 34 cylindrica), herbs, shrubs, and brush woods [55], but their share in the study area is low because these elevations are only reached further in the north and east (as indicated in Figure 1). 2) Open lands (abbreviated as "open" in the following sections)-includes all areas outside the forests which are only sparsely covered by grassy plants, barren lands, or degraded forests affected by shifting cultivation as identified by Giri and Shrestha [77]. 3) Water-the course and floodplains of the Naf River in the southeast of the study area. This also includes the frequently waterlogged marshes in this valley. 4) Informal settlements (abbreviated as "camp" in the following sections)-the study area only contains the refugee camp of Kutupalong and its extension areas. This class represents all areas that are covered by solid humanitarian facilities and simple dwellings of lighter construction materials. Radar images of S1 and A2 were combined to compute a large feature space for the supervised classification (section 2.2.4.). As many of the target classes are also dependent from topographic conditions, a digital terrain analysis was conducted to provide information on landforms at several scales (section 2.2.3.). The training areas were then applied to the input features to generate a land Radar images of S1 and A2 were combined to compute a large feature space for the supervised classification (Section 2.2.4). As many of the target classes are also dependent from topographic Remote Sens. 2019, 11, 2047 8 of 34 conditions, a digital terrain analysis was conducted to provide information on landforms at several scales (Section 2.2.3). The training areas were then applied to the input features to generate a land cover prediction map for all eleven dates to identify the changes between 2016 and 2019 (Section 3.1). An accuracy assessment was performed based on validation samples which were independently extracted from the pre-classification (Section 3.2). Finally, the importance of all input features was evaluated to estimate their contribution to the resulting land cover classification (Section 3.3).

Processing of Optical Images and Retrieval of Training Areas
As shown in Table 2, some of the optical images closest to the acquisition dates of the SAR image pairs did not cover the full study area and unavoidably contained cloud cover. They were, however, still usable to derive training areas for the latter classification (see the next section). In the first step, cloud masking was conducted on the images based on IdePix (version 2.2, [78]) for Sentinel-2 data and the fmask algorithm (version 4.0, [79]) for Landsat products. In the second step, radiometric calibration of the images to retrieve bottom of atmosphere (BOA) surface reflectance using the calibration constants provided for Landsat by Chander et al. [80] and the sen2cor processor (version 2.5.5 [81]) for Sentinel-2. All image products were harmonized to a common spatial resolution of 10 m using bilinear resampling. Compared to nearest neighbor resampling, this technique alters the original pixel values and leads to a smoother-appearing image that can lead to problems in subsequent pattern recognition techniques [82]. As the optical images are used in this study to derive radiometric indices for the identification of small representative areas of extensive land cover classes, the smoothing might even contribute to larger homogenous areas.
The following radiometric indices were computed as a feature space for the automated extraction of training areas for all investigated dates: • Normalized Difference Vegetation Index (NDVI [83], Equation (1))-This is the most frequently used vegetation index and is based on the red and near infra-red (nIR) band. It highlights the vitality of vegetation based on their chlorophyll content and was used in this study to distinguish vegetated from open or degraded areas.

•
Enhanced Vegetation Index (EVI [84], Equation (2))-This index also addresses the identification of vegetation quality, but it was additionally selected for this study because it is said to have a higher sensitivity in high biomass areas. It is based on red, blue, and nIR, as well as additional constants-C1 and C2 are used to estimate the aerosol resistance, while G is an empirical gain factor, and L allows for the adjustment of the impact of soil.

•
Modified Normalized Difference Water Index (MNDWI [85], Equation (3))-This index is based on a ratio between the reflectance in the green band and the middle infra-red (mIR) band to highlight areas that are covered by water. In this study, it was used to delineate the training areas for water.

•
Built-up Index (BUI [86], Equation (4))-This rather complex index makes use of the red band and the two short-wave infra-red bands (swIR1, swIR2) and was used to highlight built-up surfaces for the delineation of the training areas of the refugee camp.

•
Index-based Built-up Index (IBI [87], Equation (5))-This index consists of the three previously mentioned indices and is said to be more sensitive to built-up areas of different densities, which is important in this setting as parts of the refugee camp consist of different dwellings and construction materials. [31].

•
Normalized Difference Build-up Index (NDBI [88], Equation (6))-This rather simple index uses the near and short-wave infra-red bands (nIR, swIR) which also highlight urban areas. It was additionally selected to mitigate the different center wavelengths and bandwidths of Landsat ETM+, Landsat Oli, and Sentinel-2.
To minimize the effect of subjectivity and manual input in the study design, 25 points were digitized for each class, based on visual interpretation of each cloud-masked image. These 100 points were then used as input for a Classification and Regression Tree (CART [89]) to retrieve a coarse pre-classification of the four classes in each image. The trees were pruned to a maximum depth of six levels, wherein the minimum number of instances per leaf was two, and subsets containing less than five members were not split any further. This resulted in training accuracies between 87.3% and 96.9%. As shown by the trees for the single dates, which are displayed in Appendix A, the use of radiometric indices allowed for nodes to be achieved with high class purity within three levels for most of the dates, and six levels at maximum based on a small set of input rasters. The areas resulting from these classifications were then used as training inputs for the supervised classification (Section 2.2.4). Their sizes per analyzed date are summarized in Table A1 (Appendix B).

Processing of Topographic Input Features
As demonstrated by Chust et al., the integration of topographic features can significantly increase the quality of land cover classifications [90]. This was also proven in other studies and can be explained by the fact that the presence or absence of certain classes of land cover or land-use are strongly determined by elevation, slope, or aspect of surfaces [91][92][93]. For example, forest cover ends at certain elevations, water bodies only occur on flat terrain, and windward slopes are often characterized by different vegetation communities and soil types instead of downwind sides [94][95][96]. Furthermore, both optical and radar images underlie topographically-induced radiometric distortions in terms of illumination and orientation [97,98]. As shown in Figure 1, the study area shows pronounced topographical patterns, such as the ridges of over 100 m elevation ranging in the north-south direction within a narrow zone from the coast. To be able to describe the land cover with respect to these topographic conditions and to enhance the feature space of the supervised classification, a digital elevation model (DEM) at a spatial resolution of 30 m of the Shuttle Radar Topography Mission (SRTM [99]) was used to derive the following parameters: elevation, slope, northness, eastness, plan curvature, and profile curvature by means of digital terrain analysis [100], and a classification of morphometric features (peaks, ridges, passes, channels, pits, and planes [101]). In the first step, the data were resampled to a spatial resolution of 10 m using bilinear resampling to match the resolution of the final classification (see the next section). As these parameters are highly scale-dependent, this analysis was performed for window sizes of 1 × 1, 5 × 5, 9 × 9, 15 × 15, 25 × 25, and 49 × 49 pixels, leading to a total number of 42 topographic features.

Processing of Radar Input Features and Image Classification
Radar data of both S1 and A2 were radiometrically calibrated to Beta0 by applying calibration coefficients from the metadata [102,103] to retrieve radar brightness. Using the digital elevation model described in Section 2.2.3, terrain-flattened Gamma0 was derived, as introduced by Small et al. [98], to minimize the impact of the incidence angle on backscatter intensity. A refined Lee filter [104] with a window size of 5 × 5 pixels was applied to all images to reduce speckle effects. This adaptive speckle filter produces smoother images in both high and low backscatter areas and provides a good compromise between image quality and computing time [105]. For each date, image products of S1 and A2 were co-registered to a stack, based on their orbit information and geocoded using range Doppler terrain correction [106], so that data for each date consisted of all four polarizations (HH and HV from ALOS-2; VH and VV from Sentinel-1) at a common spatial resolution of 10 m.
To enhance patterns emerging at several spatial scales, the pre-processed SAR stacks were filtered by median filters with windows sizes of 5 × 5, 9 × 9, 15 × 15, 25 × 25, and 49 × 49 pixels, leading to a total number of 24 radar bands per analyzed date. As all four polarizations were available, the following radar measures were additionally derived for each date: Radar Vegetation Index (RVI [107], Equation (7)), Radar Forest Degradation Index (RDFI [108], Equation (8)), Canopy Structure Index (CSI [109], Equation (9)), average like-polarized magnitude (LK [109], Equation (10)), average cross-polarized magnitude (CS [109], Equation (11)), Volume Scattering Index (VSI [109], Equation (12)), standard deviation (StD, Equation (13)) and three combined indices (CI1, CI2, and CI3, Equation (14)), which are leaned towards the ratio of different polarizations as defined by the Pauli scattering elements [52]. Technically, all these indices can only be derived from fully-polarized radar data taken by one sensor within a single acquisition because they are sensitive to the orientation of objects on the measured surface [53]. This physical consistency is not given in our study, as the images of Sentinel-1 (VV and VH) and ALOS-2 (HH and HV) were not acquired at the same time, nor do they use the same frequency. Their interpretation regarding the initial use, namely the inherent physical characteristics of surfaces, should be made with respective care. Examples of the indices used in this study and how they are affected by the median filter products of different window sizes are given in Figure 3. It is nicely demonstrated how Figure 3A (based on median filter over 5 × 5 pixels) shows a higher level of detail, while Figure 3B (25 × 25 pixels) results in more homogenous areas favoring the classification of large patches of land cover. Figure 3C The feature space for the supervised classification consisted of 126 features per investigated date (84 radar features and 42 constant terrain features). They represent different types of landforms and land cover at multiple spatial scales, but since they are derived from a small amount of base data (four polarizations and one DEM), they are correlated with each other. To identify the features that are most suitable for the description of the study area, a random forest (RF) classifier was selected. As a method from the machine learning domain, it originates from simple classification trees, but additionally introduces a random component by repeatedly selecting random subsets from both the training samples and the input features within a defined number of independent iterations [110]. For each investigated date, 150 decision trees were calculated on the basis of one third of the available training pixels per class (Section 2.2.2), and twelve feature rasters (square root of the total 126 available features, potentially consisting of both radar and topography features). Accordingly, each decision tree used a different subset of training pixels and features to predict the distribution of the target classes (vegetation, open, water, camp) by a CART-like hierarchical ruleset. This led to a collection of 150 uncorrelated decision trees that were then used to determine the final classification of each pixel on the basis of the class majority of all iterations (Section 3.1). This had several advantages: Firstly, redundancies in the input features, as they were used in this approach, had no negative effect on the classification [111]. This is especially valuable in cases where the relevance of these features is not known in the beginning [112]. Additionally, smaller errors in the pre-classification derived from the spectral indices (Section 2.2.2) will not have a large impact by falsely training the classifier [110]. As a result of the permutation of the input features, their contribution to the final classification can be estimated, which not only identifies which features are most suitable to describe the study area, but also at which spatial scale (Section 3.3).
Remote Sens. 2019, 11, x FOR PEER REVIEW 11 of 34 these features is not known in the beginning [112]. Additionally, smaller errors in the preclassification derived from the spectral indices (section 2.2.2.) will not have a large impact by falsely training the classifier [110]. As a result of the permutation of the input features, their contribution to the final classification can be estimated, which not only identifies which features are most suitable to describe the study area, but also at which spatial scale (section 3.3.).

Land Cover Change
The land cover maps for all eleven dates are shown in Figure 4. Most strikingly, the area of the camp and the expansion of the camp after July 2017 is dominant throughout all images. The camp area increased from 36 hectares in early 2016 to over 1850 hectares at the beginning of 2019. The highest increase is observed between July 2017 and February 2018, where the camp area increased tenfold within only a few months (according to the reports, the largest immigration was between August and November 2017 [64]). Consistent with the findings of Hassan et al. [23], who reported an increase of 1219 hectares between December 2016 and December 2017, the findings of this approach seem plausible at first sight (camp growth of 1313 hectares between February 2017 and February 2018).

Land Cover Change
The land cover maps for all eleven dates are shown in Figure 4. Most strikingly, the area of the camp and the expansion of the camp after July 2017 is dominant throughout all images. The camp area increased from 36 hectares in early 2016 to over 1850 hectares at the beginning of 2019. The highest increase is observed between July 2017 and February 2018, where the camp area increased tenfold within only a few months (according to the reports, the largest immigration was between August and November 2017 [64]). Consistent with the findings of Hassan et al. [23], who reported an increase of 1219 hectares between December 2016 and December 2017, the findings of this approach seem plausible at first sight (camp growth of 1313 hectares between February 2017 and February 2018).
were correspondingly classified as open areas. Inland water bodies showed larger spatial variations because of wetlands, temporarily flooded oxbows, and creeks along the Naf River in the southeast of the study area. Seemingly all these landscape features show similar backscatter characteristics and strongly vary between the rainy and dry seasons. However, this is of no concern for the assessment of camp growth and vegetation dynamics because they are rarely covered by forest and also have no suitable building ground [33].   areas-the only image pair found of S1 and A2 at the end of the rainy season was for October 2016. As shown in the figure, the classified vegetation areas dominate the study area, while vegetation coverage showed more differentiated patterns than in all other maps. As later shown in Section 3.2, this was not caused by low classification accuracy, but rather by the strong presence of vegetation, composed not only of trees, but also bushes and shrubs that are fully developed at the end of the rainy season and therefore additionally contribute to volume scattering. The remaining land areas were correspondingly classified as open areas. Inland water bodies showed larger spatial variations because of wetlands, temporarily flooded oxbows, and creeks along the Naf River in the southeast of the study area. Seemingly all these landscape features show similar backscatter characteristics and strongly vary Remote Sens. 2019, 11,2047 13 of 34 between the rainy and dry seasons. However, this is of no concern for the assessment of camp growth and vegetation dynamics because they are rarely covered by forest and also have no suitable building ground [33].
The strong impact of phenology is also underlined in Figure 5A. It demonstrates the variation of vegetation cover during the investigated period in hectares (green line, primary y-axis) with the monthly precipitation for this area as retrieved from the Global Precipitation Climatology Centre (GPCC [113], blue areas, secondary y-axis). Although the vegetation cover might show a slightly negative trend (green line), the dynamics resulting from seasonal rainfall availability superimposed these patterns at this spatial scale. This underlined the indication that seasonal variations, or more precisely speaking, the date of the used images, had a large impact on temporal analyses in this area, and that studies using single images for comparison (see Section 1) potentially underlie a massive phenological bias that can lead to false conclusions. The water bodies in the study area ( Figure 5A, blue line) underlay the variations discussed above, but showed an overall stable development. Clearly visible is the growth of the camp after July 2017 as a result of the refugee influx starting in August 2017. Unfortunately, no further image pairs of S1 and A2 are available before February 2018 for a more detailed temporal delineation of the camp expansion, especially for the period between 25 August and 11 September 2017 when the largest influx was reported [64]. The strong impact of phenology is also underlined in Figure 5A. It demonstrates the variation of vegetation cover during the investigated period in hectares (green line, primary y-axis) with the monthly precipitation for this area as retrieved from the Global Precipitation Climatology Centre (GPCC [113], blue areas, secondary y-axis). Although the vegetation cover might show a slightly negative trend (green line), the dynamics resulting from seasonal rainfall availability superimposed these patterns at this spatial scale. This underlined the indication that seasonal variations, or more precisely speaking, the date of the used images, had a large impact on temporal analyses in this area, and that studies using single images for comparison (see section 1) potentially underlie a massive phenological bias that can lead to false conclusions. The water bodies in the study area ( Figure 5A, blue line) underlay the variations discussed above, but showed an overall stable development. Clearly visible is the growth of the camp after July 2017 as a result of the refugee influx starting in August 2017. Unfortunately, no further image pairs of S1 and A2 are available before February 2018 for a more detailed temporal delineation of the camp expansion, especially for the period between 25 August and 11 September 2017 when the largest influx was reported [64].  Figure 3D).
To highlight the camp growth, but also to show its impact on the direct surrounding environment, Figure 5B shows the progress of the classes vegetation, open, and camp for an area within a 2 kilometer distance from the main development axis of the camp (delineated in Figure 3D as a white dashed line for the axis and a gray dotted line for the area). It includes all portions of the main camp and its main expansion area in their present extent, but also the Balukhali expansion in the south. The decrease of forest areas became more clearly visible, and even though vegetation cover in October 2016 was above-average compared to the other investigated scenes, a strong decrease remained. What is interesting is that forest areas had already decreased in favor of open areas in the beginning of 2017, accordingly, before the camp expansion happened in August 2017. This can be partially explained by a preceding migration movement from Myanmar into the Cox's Bazar district in October 2016, which was denoted as a "significant influx of people" by the Inter Sector Coordination Group [64]. This movement resulted in 74,000 arrivals until February 2017, of which 31,207 were registered in the Kutupalong makeshift settlement (KMS) and reportedly "thousands of newly built huts have appeared in the hills adjacent to KMS" [114]. Thus, it is possible that the deforestation of the area close to Kutupalong already started before the massive influx of Rohingyas in August 2017 as a result of increased demand for construction materials to provide shelter and accommodation for the arrivals.   Figure 3D).
To highlight the camp growth, but also to show its impact on the direct surrounding environment, Figure 5B shows the progress of the classes vegetation, open, and camp for an area within a 2 km distance from the main development axis of the camp (delineated in Figure 3D as a white dashed line for the axis and a gray dotted line for the area). It includes all portions of the main camp and its main expansion area in their present extent, but also the Balukhali expansion in the south. The decrease of forest areas became more clearly visible, and even though vegetation cover in October 2016 was above-average compared to the other investigated scenes, a strong decrease remained. What is interesting is that forest areas had already decreased in favor of open areas in the beginning of 2017, accordingly, before the camp expansion happened in August 2017. This can be partially explained by a preceding migration movement from Myanmar into the Cox's Bazar district in October 2016, which was denoted as a "significant influx of people" by the Inter Sector Coordination Group [64]. This movement resulted in 74,000 arrivals until February 2017, of which 31,207 were registered in the Kutupalong makeshift settlement (KMS) and reportedly "thousands of newly built huts have appeared in the hills adjacent to KMS" [114]. Thus, it is possible that the deforestation of the area close to Kutupalong already started before the massive influx of Rohingyas in August 2017 as a result of increased demand for construction materials to provide shelter and accommodation for the arrivals. This indication is underlined by the drastic increase of open areas around the center of the camp. Temporally delayed, these open areas were then converted to densely-populated camp sites after July 2017.
The decrease of forest areas in the direct surroundings of the camp matches the 2283 hectares reported by Hassan et al. [23] when the dates of July 2016 and February 2018 are compared (decrease of 2100 hectares), excluding the data from October 2016 (extraordinarily high amount of vegetation). As expected, the camp area increased the most between July 2017 and February 2018 (+1314 hectares) when most of the Rohingya crossed the border (see Section 1). A small reduction of camp area is observed between June and July 2018, which can partially be credited to demolitions "damaging thousands of tents" by monsoon rains in early June [115]. Nonetheless, it is most likely that smaller misclassifications played a role as well. They will be discussed in the next section.
Comparing the first and the last image utilized in this study (both acquired in March) a total loss of 1562 hectares of vegetation can be attested, of which 20% are now camp areas and 80% are open areas. This testifies both to the impact of the camp, but also the overall decrease of vegetation during the investigated period due to increased resource demand [33]. The camp area itself increased by 1830 hectares, which is demonstrated in Figure 6 As expected, the camp area increased the most between July 2017 and February 2018 (+1314 hectares) when most of the Rohingya crossed the border (see section 1). A small reduction of camp area is observed between June and July 2018, which can partially be credited to demolitions "damaging thousands of tents" by monsoon rains in early June [115]. Nonetheless, it is most likely that smaller misclassifications played a role as well. They will be discussed in the next section.
Comparing the first and the last image utilized in this study (both acquired in March) a total loss of 1562 hectares of vegetation can be attested, of which 20% are now camp areas and 80% are open areas. This testifies both to the impact of the camp, but also the overall decrease of vegetation during the investigated period due to increased resource demand [33]. The camp area itself increased by 1830 hectares, which is demonstrated in Figure 6  A different form of change visualization is presented in Figure 7. This approach was originally suggested by Langer [116,117], who investigated landscape changes around the Lukole refugee camp in Tanzania based on a time-series of 42 Landsat images. For this study, squares of the size of 300 meters width and height were generated, and the dominant landcover was assessed for each square and date on the basis of the determination of each square's mode. The square was then divided into nine sub-squares the size of 100 meters to which the land cover classifications between July 2016 and January 2019 were assigned. Accordingly, the change of landscapes can be visualized within each square of 300 meters, as demonstrated in Figure 7. The initial area of the camp is characterized by consistently red squares, while the expansion sites south of it are gradually colored red in the lower parts of the squares. Additionally, the areas of conversion from forest to camp are well identifiable in the map on the right side (pink tones). They make up about 630 hectares, which correlates well with the 572 hectares of forest decrease reported by Labib et al. [35], who estimated the conversion A different form of change visualization is presented in Figure 7. This approach was originally suggested by Langer [116,117], who investigated landscape changes around the Lukole refugee camp in Tanzania based on a time-series of 42 Landsat images. For this study, squares of the size of 300 m width and height were generated, and the dominant landcover was assessed for each square and date on the basis of the determination of each square's mode. The square was then divided into nine sub-squares the size of 100 m to which the land cover classifications between July 2016 and January 2019 were assigned. Accordingly, the change of landscapes can be visualized within each square of 300 m, as demonstrated in Figure 7. The initial area of the camp is characterized by consistently red squares, while the expansion sites south of it are gradually colored red in the lower parts of the squares. Additionally, the areas of conversion from forest to camp are well identifiable in the map on the right side (pink tones). They make up about 630 hectares, which correlates well with the 572 hectares of forest decrease reported by Labib et al. [35], who estimated the conversion between March and December 2017. The trends in Figure 7 (right) were calculated based on the first three (July 2016 to February 2017) and last three (June 2018 to January 2019) images of each square-it shows areas where vegetation has permanently changed to camp (light pink) and or to open areas (violet). The latter could also be some indication for transitions to camp areas, especially in the southwest and south of Kutupalong where new settlements are emerging [115]. Most of the forested areas of the Teknaf Wildlife Sanctuary in the south of the study area are characterized by permanent vegetation cover (similar to the areas shown in Figure 6), but two sub-squares indicate the clearance of forests, probably as a result of the population increase, as reported by Imtiaz [34]. . The latter could also be some indication for transitions to camp areas, especially in the southwest and south of Kutupalong where new settlements are emerging [115]. Most of the forested areas of the Teknaf Wildlife Sanctuary in the south of the study area are characterized by permanent vegetation cover (similar to the areas shown in Figure 6), but two sub-squares indicate the clearance of forests, probably as a result of the population increase, as reported by Imtiaz [34].

Accuracy Assessment
Because decisions are to be based on the findings of this study and other studies on the forest change, as outlined in section 1, the accuracy of the generated output maps must be known [118]. It strongly determines how well conclusions can be drawn from the retrieved land cover dynamics and at which points there might be weaknesses. Assessing the accuracy of digital maps was based on independently collected validation points which were sampled according to specific rules. One way to determine the required number of samples N according to the expected accuracy and confidence level is given by Fitzpatrick-Lins [119] (Equation (15)) based on the binomial distribution formula: where Z is the standard normal deviate for a 95% confidence level (here: 2), p is the expected accuracy of the map (here: 90%), q is 100 − p (here: 10), and E is the allowable error (here: 2.5%). This led to a required number of 576 points per analyzed date. As the different land cover classes were not equally distributed over the study area, we applied a stratified random sampling so that the

Accuracy Assessment
Because decisions are to be based on the findings of this study and other studies on the forest change, as outlined in Section 1, the accuracy of the generated output maps must be known [118]. It strongly determines how well conclusions can be drawn from the retrieved land cover dynamics and at which points there might be weaknesses. Assessing the accuracy of digital maps was based on independently collected validation points which were sampled according to specific rules. One way to determine the required number of samples N according to the expected accuracy and confidence level is given by Fitzpatrick-Lins [119] (Equation (15)) based on the binomial distribution formula: where Z is the standard normal deviate for a 95% confidence level (here: 2), p is the expected accuracy of the map (here: 90%), q is 100 − p (here: 10), and E is the allowable error (here: 2.5%). This led to a required number of 576 points per analyzed date. As the different land cover classes were not equally distributed over the study area, we applied a stratified random sampling so that the number of samples represented the spatial share of each class-250 points for vegetation and 100 points for water (in accordance with the 41% of forest cover and 19.75% of floodplain in the greater study area, as reported by Moslehuddin [55]). One hundred and fifty points were sampled within areas classified as open, as the second largest class and the remaining 76 points were sampled within the camp areas, which did not represent their spatial share in the resulting map, but it is granted that enough sampling points lay within this class to provide a robust accuracy assessment. The points were automatically generated within each of the 11 analyzed maps and manually assigned to a class on the basis of visual interpretation of the corresponding optical image ( Producer's accuracy (PA) refers to the commission error and indicates how reliable the classification of a pixel is, while user's accuracy (UA) refers to the omission error and indicates the completeness of the determination of a class [122]. Cohen's kappa additionally determines how likely the agreement between validation samples and classification occurred by chance and is often claimed to be a more robust measure than the simple percent agreement [123]. Table 3 shows that the high vegetation cover in October 2016 was not a result of over-estimation, as it can be attributed to the time of image acquisition which lies at the end of the rainy season with temperatures above 30 • C and precipitation of up to 1000 mm per month, leading to a pronounced vegetation cover [113]. Nonetheless, the PA of vegetation of 93.6 indicated a slight over-estimation of this date, but it was comparable to those of other dates, ranging from 85.2% to 95.9% in total for the vegetation class. Accordingly, the classified vegetation cover can be trusted throughout the study to a high degree, which was additionally confirmed by the kappa values between 0.8 and 0.9. Throughout almost all investigated dates, open areas had the lowest producer's accuracies, which led to the conclusion they were over-estimated in many cases. This was likely at the edges of large connected patches of vegetation, in areas where the vegetation cover is spatially fragmented, or areas which are generally characterized by lower vegetation, such as shrubs or bushes. The lowest overall accuracies observed in this study were achieved for 2016-03, 2017-07, and 2019-03, which were caused by falsely classified open areas (UA of 72.2%, 84.0%, and 74.7%, respectively) and under-estimated vegetation areas (PA of 86.8%, 89.0%, and 85.2%, respectively).
Lastly, the camp areas were classified with very high kappa values and low errors of commission and omission. This means the expansion of the camps and the accompanying decrease of vegetation not only corresponded to the values found in other studies (see Section 3.1) but they were also well described by SAR data in this approach, most probably due to an increase of dihedral scattering caused by double bounce effects at artificial structures.
One drawback to mention here is that the observed accuracy was a result of the study design, starting with the selection of 25 training samples for the pre-classification based on the radiometric indices. Although the fundamental idea of this study is based on the assumption that this number is enough for classifying an entire image based on radar and topographic features, it has to be noted that small changes in these training samples potentially had a large impact on the final results. This was partly because they defined the reference areas (as derived from radiometric indices) for the training of the actual land cover classification, leading to different rulesets within the classification algorithm. Moreover, their small number was sensitive to outliers which did not fully represent the respective class. As already indicated in Section 2.2.4, this risk was minimized by the massive permutation of input features within the random forest classifier, making this classifier robust to outliers and favoring features with a higher significance to distinguish the target classes [110]. Also, as shown by the independently computed overall accuracies, beginning with only 25 visually identified training points only can lead to very good results.
As for the stratified random sampling used in the accuracy assessment, it is suggested by Townshend and Justice that, instead of using single pixels, clusters of connected pixels (e.g., blocks of 3 × 3 pixels) should be analyzed as a more robust measure for the representation land cover classes and their spatial distribution of land cover classes [124]. While this form of sampling clearly reduces the costs for data collection at an early stage of the study [121], it increases the number of samples to be manually labelled nearly 10-fold in the later validation part.
To put the findings of this study into the context of the current scientific debate, we compared the land cover maps of this study to those created by Hassan et al. in 2018 [23]. Figure 8A shows the pre-influx agreement between the results from Hassan et al., based on optical images from December 2016), with the closest date analyzed in this study from February 2017. The overall statistical agreement was only 66.7% (pre-influx, Figure 8A) and 60.8% (post-influx, Figure 8B), which seems low at first sight. However, the comparison reveals some interesting points about the quality of this study that will be discussed as follows. Figure 8 shows that the camp area was nearly the same in both studies, but larger parts, which were classified as vegetation in this study, were open areas in the study of Hassan et al. (Figure 8A bright blue areas). This can only partly be explained by an overestimation of vegetation areas in our study, as indicated by the producer's accuracy of 93.3% (Table 3). More importantly, the temporal delay (here: 48 days) between both observations might already contain phenological effects caused by the rainy season gradually starting in February ( Figure 5A), resulting in more areas considered as healthy vegetation. Lastly, both approaches were based on different input data-while the study of Hassan et al. used multi-spectral satellite data that mostly measured the amount of chlorophyll of a surface using infra-red information, this study made use of radar data, which is sensitive to backscatter mechanisms. That means that areas which had low chlorophyll values in the optical images might still be recognized as vegetation by the radar images because of their branch structures and crown surfaces. camp class from February 2018 because in the image used for training (13.02.2018) the classifier was closer to the radar image (01.02.2019) than the one used by Hassan et al. (15.12.2017). Nonetheless, the ambiguity in the backscatter intensity of solid and light dwellings can become a problem in studies dealing with humanitarian settings, as already discussed by Braun [125], however, more research is required to understand the differences between their radar signatures [126].

Feature Importance
Besides the land cover maps illustrated in Figure 4, the permutation importance (PI) and the Gini decrease (GD), as measured for the contribution of each feature (input raster), were calculated for the random forest classification. PI measured the mean decrease in the model accuracy by systematically comparing the effect of randomly shuffled input variables [127]. Accordingly, the higher a feature's permutation importance, the more important it was for the respective class. GD made use of the Gini impurity, which determined how well a split separates the permuted samples at each node to estimate each variable's importance on the basis of its occurrence throughout all trees [128]. As an index value, it had no upper boundary, but the higher the GD, the larger the importance of the variable for the classification. For an evaluation of the feature importance, the 500 features with the highest PI throughout all investigated dates were selected. Of the features, 82 were based on topographic parameters (section 2.2.3.) and 418 were SAR-based features (section 2.2.4.). This confirms the design of the study, which combined spectral information with topographic characteristics to predict the land cover classes, where the radar parameters played a significantly larger role. As shown in Table 4, the most frequent feature used throughout all the analyzed dates was VH polarization of Sentinel-1, followed by VV and HH. Accordingly, backscatter intensity of the SAR satellites already had a large informative value for the classification of land cover in the study area. However, indices such as the Radar Vegetation Index (RVI) and the Radar Forest Degradation Index (RFDI) strongly contributed to higher prediction accuracies. As expected, information of the cross-polarized channels (VH, HV, cross-polarized magnitude CS) played a larger role because of their sensitivity to volume scattering caused by vegetation. Table 4 also shows that all features based on topography were of secondary importance, most likely because they were stable throughout the whole period, while the seasonal variations could only be described by the changes of SAR The same applies to the post-influx comparison ( Figure 8B), which also shows a considerable discrepancy between vegetated and open areas along their borders. Two further differences can be observed here: Firstly, the camp is larger in the western part in the map of Hassan et al. (Figure 8B dark blue). These areas concentrate at the western and southern borders of the main camp area (part of the Balukhali expansion area), but also of the smaller site in the south, and it is likely that they contain dwellings which are the most recently constructed. As resources of solid materials have become short since August 2017 [32], one explanation for their under-representation in this study could be that these buildings are no longer covered by corrugated metals, or at least to a smaller degree, as a consequence of material shortage and a shift to lighter construction materials, such as wood, leaves, or plastic (all used for covering the roofs [31]). Furthermore, less densely built parts of the camp were prone to under-estimations of this land cover class. This slight under-estimation was also indicated, but not fully explained by the user's accuracies (92.1%) and kappa value (0.9) of the camp class from February 2018 because in the image used for training (13.02.2018) the classifier was closer to the radar image (01.02.2019) than the one used by Hassan et al. (15.12.2017). Nonetheless, the ambiguity in the backscatter intensity of solid and light dwellings can become a problem in studies dealing with humanitarian settings, as already discussed by Braun [125], however, more research is required to understand the differences between their radar signatures [126].

Feature Importance
Besides the land cover maps illustrated in Figure 4, the permutation importance (PI) and the Gini decrease (GD), as measured for the contribution of each feature (input raster), were calculated for the random forest classification. PI measured the mean decrease in the model accuracy by systematically comparing the effect of randomly shuffled input variables [127]. Accordingly, the higher a feature's permutation importance, the more important it was for the respective class. GD made use of the Gini impurity, which determined how well a split separates the permuted samples at each node to estimate each variable's importance on the basis of its occurrence throughout all trees [128]. As an index value, it had no upper boundary, but the higher the GD, the larger the importance of the variable for the classification. For an evaluation of the feature importance, the 500 features with the highest PI throughout all investigated dates were selected. Of the features, 82 were based on topographic parameters (Section 2.2.3) and 418 were SAR-based features (Section 2.2.4). This confirms the design of the study, which combined spectral information with topographic characteristics to predict the land cover classes, where the radar parameters played a significantly larger role. As shown in Table 4, the most frequent feature used throughout all the analyzed dates was VH polarization of Sentinel-1, followed by VV and HH. Accordingly, backscatter intensity of the SAR satellites already had a large informative value for the classification of land cover in the study area. However, indices such as the Radar Vegetation Index (RVI) and the Radar Forest Degradation Index (RFDI) strongly contributed to higher prediction accuracies. As expected, information of the cross-polarized channels (VH, HV, cross-polarized magnitude CS) played a larger role because of their sensitivity to volume scattering caused by vegetation. Table 4 also shows that all features based on topography were of secondary importance, most likely because they were stable throughout the whole period, while the seasonal variations could only be described by the changes of SAR backscatter. Still, the elevation and exposition of areas seemed to play a role in the study area, most probably for water bodies and permanent forest cover. Table 4. Frequency of the 500 most important features used for the classification. DEM = digital elevation model.

Feature Origin Frequency
VH SAR (S1) 85 VV SAR (S1) 51 HH SAR (A2) 43 Radar Vegetation Index (RVI) SAR (S1 + A2) 40 Cross-polarized magnitude (CS) SAR (S1 + A2) 38 Average like-polarized magnitude (LK) SAR (S1 + A2) 26 Combined index (CI1) SAR (S1 + A2) 23 HV SAR (A2) 21 Standard deviation (StDev) SAR (S1 + A2) 20 Radar Forest Degradation Index (RFDI) SAR (S1 + A2) 20 Canopy Structure Index (CSI) SAR (S1 + A2) 19 Volume Scattering Index (VSI) SAR (S1 + A2) 18  Elevation  DEM  18  Aspect: northness  DEM  18  Morphometric feature  DEM  17  Slope  DEM  17  Combined index (CI3) SAR (S1 + A2) 14 Profile curvature DEM 6 Aspect: eastness DEM 6 In the context of the features used for the classification and their contribution to the final result, it must also be discussed that the indices derived from combinations of both SAR sensors were potentially distorted by the relatively stronger level of backscatter from the L-band compared to the C-band. Accordingly, these indices might be dominated by the contribution of HH and HV (ALOS-2). In this study, the magnitude of input features probably diluted this effect, but the actual difference of both sensors should be investigated in further studies, probably raising the need for a proper mutual normalization of both sensors.
The values in Table 5 show the permutation importance per class, namely how much worse (in percent of the achieved training accuracy) a classifier would perform if the named feature had random values. The table shows that the RFDI and the VSI were the most important variables to discriminate vegetation from open areas, but also that topographical features played a role here, as the forests are restricted to certain heights and are absent on steep slopes [55]. Also, as indicated above, the combined index CI2, potentially highlighting complex surfaces, such as crowns and other volumes, was among the five most important variables for the identification vegetation. The water class was mainly determined by areas with low elevation and small slopes, but also very low values in the HH polarization. Lastly, the camp was characterized by high HH values, but also by CI1, as was probably retrieved from the many horizontal structures. The high importance of the standard deviation among all four polarizations (StDev) again underlined the complex backscatter characteristics of settlements of light construction materials. Aside from the thematic information content of the input data for the classification, all features were evaluated regarding the window size, which was used for their computing. This refers to the scale of the digital terrain analysis (Section 2.2.3) and the pixel size of the median filter applied to all SAR data (Section 2.2.4). Figure 9 shows that large window sizes were selected more frequently during the random forest classification compared with small ones. Of the 500 most important features as listed in Table 4, 124 were of the largest size of 49 pixels (corresponding to a spatial scale of 490 m). Generally, large window sizes have higher prediction values than smaller ones, most probably because they are more suitable to predict large consistent areas. This was already discovered by Braun et al. [48], who reported a dominance of large-scale features for land cover assessments around the refugee camp of Dagahaley in Kenya. As another possible explanation, products of smaller window sizes are often characterized by topographic variations or speckle effects, which do not represent the spatial distribution of the expected land cover. In our study, features with window sizes of five pixels slightly exceeded the linear trend. Considering the spatial resolution of the used raster, this corresponded to a size of 50 m, which apparently played a more important role in the study area. Consequently, it might be sufficient for future approaches to derive filter products of window sizes of 5, 25, and 49 pixels to keep computing time low, while results of similar prediction quality can be expected.  Aside from the thematic information content of the input data for the classification, all features were evaluated regarding the window size, which was used for their computing. This refers to the scale of the digital terrain analysis (section 2.2.3.) and the pixel size of the median filter applied to all SAR data (section 2.2.4.). Figure 9 shows that large window sizes were selected more frequently during the random forest classification compared with small ones. Of the 500 most important features as listed in Table 4, 124 were of the largest size of 49 pixels (corresponding to a spatial scale of 490 meters). Generally, large window sizes have higher prediction values than smaller ones, most probably because they are more suitable to predict large consistent areas. This was already discovered by Braun et al. [48], who reported a dominance of large-scale features for land cover assessments around the refugee camp of Dagahaley in Kenya. As another possible explanation, products of smaller window sizes are often characterized by topographic variations or speckle effects, which do not represent the spatial distribution of the expected land cover. In our study, features with window sizes of five pixels slightly exceeded the linear trend. Considering the spatial resolution of the used raster, this corresponded to a size of 50 meters, which apparently played a more important role in the study area. Consequently, it might be sufficient for future approaches to derive filter products of window sizes of 5, 25, and 49 pixels to keep computing time low, while results of similar prediction quality can be expected. Lastly, the importance of specific features was evaluated by calculating the average PI and GD for each combination of feature and window size throughout all 11 investigated dates. Table 6 and  Table 7 summarize the 10 combinations of features and window sizes with the highest PI and GD scores. Interestingly, elevation 49 had the highest average PI and GD scores, but only occurred in three of the 11 analyzed dates. This indicated a strong correlation between some of the automatically selected samples, as described in section 2.2.2., with topographical features. Interestingly, the Volume Lastly, the importance of specific features was evaluated by calculating the average PI and GD for each combination of feature and window size throughout all 11 investigated dates. Tables 6 and 7 summarize the 10 combinations of features and window sizes with the highest PI and GD scores. Interestingly, elevation 49 had the highest average PI and GD scores, but only occurred in three of the 11 analyzed dates. This indicated a strong correlation between some of the automatically selected samples, as described in Section 2.2.2, with topographical features. Interestingly, the Volume Scattering Index (VSI) and the Radar Forest Degradation Index (RFDI) were among the most important features (Table 4) and showed generally large PI and GD scores (Tables 6 and 7), but CI2 was not among these lists. Its role for vegetation is indicated in Table 5, however, attesting a high significance for the vegetation class. We conclude that VSI and RFDI may be the more important variables for the entire study area as a whole, as they also provide an explanation for the spatial distribution of the other classes (open, water, and camp). Furthermore, the information content of the combined indices (CI) was not directly representative of the scattering mechanisms, as was predictable by a canonical polarimetric matrix derived from fully-polarimetric sensors, as discussed in Section 2.2.4. The high PI scores of slope (sizes 25 and 49) could result from their high importance for the classification of water bodies, which only occur in flat areas. Of all original SAR backscatter features, the co-polarization HH of ALOS-2 had the highest frequency and the highest average variable importance. Furthermore, it was the only variable of this type among the 10 with the highest average scores based on window sizes of 5 and 15 pixels. This confirms the idea that even if most of the important variables were derived from large window sizes, a combination of features derived at different scales would lead to the most promising result. Despite the many indications on variable importance, we did not observe a prevalence of one sensor for the final result. Moreover, the combined measures (e.g., RFDI, VSI) were the ones that systematically achieved the highest scores. For the discrimination of forests from grasslands, the contribution of the sensor's different wavelengths would have become more evident. This was already demonstrated by Braun et al. [129], who combined short and long radar signals of active and passive sensors to derive regional maps of above-ground biomass in Senegal. By the time this article was written, the only satellites capturing radar images at the L-band were ALOS-2 and SAOCOM-1A, which are commercially distributed by JAXA and CONAE respectively [130]. This means that, currently, no L-band data is freely available (besides archived data of ALOS PALSAR between 2007 and 2011). This can be an obstacle for humanitarian organizations that do not have a budget to purchase satellite data. The ALOS-2 products used in this study cost ¥240,000 (roughly 2000 Euro) per scene [131]. Accordingly, time-series analyses such as this one can become expensive, and the transferability of this approach is limited in this respect. However, many organizations from the humanitarian domain have established agreements with operators and distributors so to get data at a discount, or even for free [132], or have increased their budget for digital technologies [133,134]. Lastly, new satellites will emerge in the upcoming years that will provide images acquired at the L-band (NISAR [2021], TanDEM-L [2023]), and even P-band (BIOMASS [2012]), entirely for free [130]. Against this background, it is important that methods already exist that address the special challenges of sensor combinations within humanitarian settings.
To conclude, the random forest classifier effectively handled the partially-redundant feature space and selected the ones which were most suitable towards predicting the land cover classes on the basis of the automatically derived training samples. This makes this approach largely operational, as it is no longer dependent upon an extensive collection of training samples by experts, nor from the manual selection of input features for the classification.

Conclusions and Outlook
This study showed how radar data of different frequencies can be integrated for land cover assessments and the monitoring of environmental conditions in situations of humanitarian crisis. It complements existing studies that, because of the limited availability of optical data, only roughly estimate the decrease of forests and the growth of camps. This study also introduces a high level of automation, as user input is only needed at the beginning for the selection of a few representative reference samples. As demonstrated in this study, a land cover classification with accuracies above 90% can be achieved by the selection of only 25 sample points per class. This number was proven to be a good trade-off between scientific accuracy and transparency, but also applicability and reproducibility. Applying radiometric indices on optical data to construct simple CARTS was used to identify all areas which fulfil similar criteria and therefore were suitable for training the classifier, even if the optical data were obscured by clouds, gaps, or did not fully cover the area of interest. Especially in cases of emergencies, this brings enormous benefits, as the workflow can be transparently transferred and automatically applied to other data sources or different areas. The only decision to be made by the users is the selection of target classes and the identification of a small amount of representative training areas.
It is yet to be tested if other settings, for example, those based on different sensor combinations or located within other ecosystems, in achieving the same level of accuracy as observed in this study.
The use of microwaves additionally increased the quality of the land cover assessments, as they allowed a higher number of input images, especially in the rainy seasons. This helps to consider phenologic variations of the images and reduces their potential bias impacting the results as the images are analyzed as parts of a time series instead of a simple pre-and post-event change detection based on a single image pair. The benefit of time-series analyses over simple binary approaches has been demonstrated in numerous studies and should therefore also be integrated into monitoring frameworks of the humanitarian domain [7,47,[135][136][137]. Furthermore, radar images can increase the quality of land cover classifications because they measure physical surface characteristics instead of spectral reflectance. This makes them especially valuable in areas of high vegetation cover where higher vegetation and grasslands have similar spectral characteristics. Among other features, this was accomplished by the computation of combined indices (CI) based on polarizations of both sensors. Although they were not based on the complex information of fully-polarimetric imagery, their role for the classification of vegetation and camp areas was demonstrated by their high feature importance-vegetation was found to be characterized best by the Volume Scattering Index (VSI) and the Radar Forest Degradation Index (RFDI), and the combined index 2 (CI2) seemed to be more sensitive towards the volumes and complex structures of crowns. In turn, the combined index 1 (CI) element was found to play a large role for the determination of camp areas where horizontal objects caused higher direct backscatter intensity ( Table 5).
The study also showed that topography played a role in the identification of classes, especially in cases where large areas were covered by homogenous classes, such as water bodies and large evergreen forest areas. Topographic features turned out to be of largest information content when retrieved at large window sizes of 25 or 49 pixels, which correspond to distances of 250 to 500 m. It appeared that these were spatial scales that represented the connection between land cover and landforms best. This was underlined by the fact that 16% of the highest features were DEM-based, and were also slope and aspect features. However, these factors were constant throughout all analyzed dates and therefore only contributed to certain areas, while the seasonal variations were highlighted by changes in backscatter intensities and its derivations.
All classified dates showed high classification accuracies, especially for water and camp areas. A comparison with a reference study [23] showed that the camp areas might have been slightly under-estimated at the edges of the expansion areas where the dwelling density is low. This, as well as the high impact of seasonal variations, has to be considered when land cover maps are evaluated for the monitoring of camp expansion and forest degradation. Furthermore, the benefit of a higher amount of input images became evident here. Although the temporal difference between optical images used for the selection of training areas and the radar images used for classification was kept small, this delay potentially introduces errors in the training step initially, which significantly reduces the quality of the later classification. In an ideal case, optical and radar images are acquired at the same date. With the increasing availability of both radar (e.g., the upcoming Radarsat constellation mission [138] plus many of the planned SmallSAT missions [139]) and multi-spectral sensors (e.g., Planet Scope [140]), this will be more likely in the near future. Nonetheless, a minimum number of image pairs over the investigated period is required, as well as a budget to cover the acquisition costs of commercial satellite data.
Of course, the high classification accuracies are also a product of a low number of target classes. In this study, a discrimination of vegetation, open, camp, and water was sufficient, but in ecosystems with more complex land cover or different seasonal dynamics, a higher number of target classes (e.g., grassland, shrubs, but also solid urban areas vs. camp dwellings of light construction materials) should be considered. The limits regarding the number of classes and usable image pairs over a longer period are to be tested in further studies.
In a humanitarian context, the integration of SAR data not only gives higher flexibility regarding the number of images and the definition of suitable acquisition dates, but also allows for the identification of landcover types which are related to characteristics of surfaces (size, form, roughness, material, and scattering mechanisms). This makes them a valuable source of information for future emergencies, given that enough studies exist that demonstrate their applicability and transferability. More research is needed to investigate the benefits of radar data in humanitarian settings.
The presented insights on the development of vegetation patterns and the growth of the Kutupalong refugee camp are of higher detail than in all previously published studies. In particular, the incorporation of data from all seasons helps to understand the phenological variations and puts the overall developments into context. The suggested approach can be automated to a very large degree so that it can feed new data into monitoring frameworks, such as those conducted by the UNDP [56], to ensure the protection of the Teknaf Wildlife Sanctuary and the livelihood of the people living in this area. Edward Cahill for the language editing of this article and both anonymous reviewers for their extensive comments and constructive suggestions. We acknowledge support by the Open Access Publishing Fund of University of Tübingen.

Conflicts of Interest:
The authors declare no conflict of interest.

Appendix A
The trees show the rulesets based on the radiometric indices to assign classes to the sampling points based on thresholding. Each node consists of a base class (middle) and its current representation within this node (indicated in percent and in absolute numbers at the bottom). The feature variable and the threshold used for the split is given at the top of the subsequent node. As shown in the figures, most of the classes in the lowest level showed purities above 85%. Depending on their purity, the nodes are colored in hues of green (vegetation), blue (open), orange (water), and red (camp). Please note that camp areas were not identified for March 2016, July 2016, February 2017, and July 2017 because of low spatial presence and cloud cover. For these dates, training areas had to be delineated manually on the basis of interpretations by experts from partner organizations from the humanitarian domain with regional expertise in the study area.

Appendix A
The trees show the rulesets based on the radiometric indices to assign classes to the sampling points based on thresholding. Each node consists of a base class (middle) and its current representation within this node (indicated in percent and in absolute numbers at the bottom). The feature variable and the threshold used for the split is given at the top of the subsequent node. As shown in the figures, most of the classes in the lowest level showed purities above 85%. Depending on their purity, the nodes are colored in hues of green (vegetation), blue (open), orange (water), and red (camp). Please note that camp areas were not identified for March 2016, July 2016, February 2017, and July 2017 because of low spatial presence and cloud cover. For these dates, training areas had to be delineated manually on the basis of interpretations by experts from partner organizations from the humanitarian domain with regional expertise in the study area.

Appendix A
The trees show the rulesets based on the radiometric indices to assign classes to the sampling points based on thresholding. Each node consists of a base class (middle) and its current representation within this node (indicated in percent and in absolute numbers at the bottom). The feature variable and the threshold used for the split is given at the top of the subsequent node. As shown in the figures, most of the classes in the lowest level showed purities above 85%. Depending on their purity, the nodes are colored in hues of green (vegetation), blue (open), orange (water), and red (camp). Please note that camp areas were not identified for March 2016, July 2016, February 2017, and July 2017 because of low spatial presence and cloud cover. For these dates, training areas had to be delineated manually on the basis of interpretations by experts from partner organizations from the humanitarian domain with regional expertise in the study area.

Appendix A
The trees show the rulesets based on the radiometric indices to assign classes to the sampling points based on thresholding. Each node consists of a base class (middle) and its current representation within this node (indicated in percent and in absolute numbers at the bottom). The feature variable and the threshold used for the split is given at the top of the subsequent node. As shown in the figures, most of the classes in the lowest level showed purities above 85%. Depending on their purity, the nodes are colored in hues of green (vegetation), blue (open), orange (water), and red (camp). Please note that camp areas were not identified for March 2016, July 2016, February 2017, and July 2017 because of low spatial presence and cloud cover. For these dates, training areas had to be delineated manually on the basis of interpretations by experts from partner organizations from the humanitarian domain with regional expertise in the study area.                                    Appendix B

Appendix C
Confusion matrix between reference (columns) and classified (rows) land cover for all investigated dates. V = vegetation, O = open, W = water, C = camp.