DETECTION OF LAND COVER DISPLACEMENTS THROUGH TIME-SERIES ANALYSIS OF MULTISPECTRAL SATELLITE IMAGERY: APPLICATION TO DESERT SAND DUNES

: Detection of changes occurring on Earth surface has become an established practice in remote sensing science since the availability of high-resolution, global coverage, and multitemporal satellite imagery that has constantly increased during the last decades. Meanwhile, the open data policies embraced by some of the principal Earth Observation programs have boosted the spread of such analyses. This asset has attracted also the interest of the private sector which is developing strategies to exploit the potential of these data. In view of the above, we present an experimental procedure to investigate land cover displacements through an application on open multispectral imagery from the Landsat 8 and Sentinel-2 missions for desert sand dunes movements analysis. While most of the change detection techniques focus on locating changes and describing them by means of variation in the pixel spectral responses, the proposed technique aims at describing spatial and temporal patterns of displacements (i.e. directions and magnitude) applying cross-correlation analysis on a multitemporal images stack. Results of the proposed analysis are critical to a number of construction engineering operations that require sand mitigation planning and continuous site monitoring to prevent windblown sand interactions with infrastructures in the desert environment. An overview of the preliminary results is presented together with an extensive discussion on further improvements requested by the procedure. Computational steps leverage exclusively open data and free and open source GIS software thus providing large rooms to empower, replicate and improve thereof.


INTRODUCTION
The emerging availability of high-resolution, global coverage, and multitemporal imagery -captured by the modern satellite missions -has empowered a number of studies focusing on the detection of changes occurring on the Earth's surface. On the one hand, the increased interest in multitemporal imagery analysis is mainly triggered by the rising number of satellites with short revisiting time that enables the acquisition of either long time series or frequent bi-temporal images (Ghamisi et al., 2019). On the other hand, the open distribution policies of both archive data as well as new observations -adopted by some of the principal Earth Observation programs -enable outstanding opportunities for the user community to grow and, in turn, facilitate the spread of case studies exploiting multitemporal imagery (Belward, Skøien, 2015).
Best-known examples of the above are represented by the NASA Landsat (NASA, 2020) and the ESA Sentinel (ESA, 2020) missions which fully embrace the open data principles. By focusing on the aforementioned satellite missions, meaningful outcomes have been achieved for the monitoring of largescale environmental phenomena including e.g. urban sprawl (Wendl et al., 2018), deforestation (Finer et al., 2018), and polar ice degradation (Runge, Grosse, 2019). However, the high spatial and spectral resolution -nowadays provided by the open multispectral satellite imagery -also allows for applications at a finer scale. These assets are further attracting the interest of * Corresponding author the private sector (Davids, Rouyet, 2018) which is developing strategies to exploit the potential of these data (Laefer, 2020).
In view of the above, we present here an experimental procedure to investigate land cover displacements by means of multitemporal analysis of multispectral imagery. The procedure is applied to the assessment of sand dune movements in desert areas. The proposed application is critical to a number of construction engineering operations that require planning of mitigation measures to prevent undesired windblown sand interactions with infrastructures that are built in the desert environment. The interest in applying open satellite imagery for sand mitigation planning is justified by the significantly higher costs of the traditional field surveys that are often inadequate to understand sand dynamics on a large scale and, preliminarily, to identify the locations where such sand mitigation measures are to be applied. This is especially true for wide construction areas such as transport infrastructure sites (Wickramasinghe et al., 2018). A number of authors in the literature have exploited the multitemporal analysis of multispectral imagery for assessing dunes movement (Hugenholtz et al., 2012). However, most of these studies focused on single dune fields, where only a few target dunes were tracked, see e.g. (Al-Mutiry et al., 2016), (Michel et al., 2018), and (Hassoup, 2019). Indeed, the proposed procedure aims at identifying movements at a pixel level but suitable to wide scenes analysis entirely covering the extension of transport infrastructure sites which usually cross multiple dune fields as well as desert land cover textures.
The procedure is tested on two study areas, located respect-ively in Egypt and United Arab Emirates (UAE) deserts, exploiting multitemporal and multispectral satellite imagery from the Landsat 8 and the Sentinel-2 missions. The computational work is carried out exclusively through the use of Geographic Information Systems (GIS) Free and Open Source Software (FOSS). The procedure consists of a single workflow that encompasses input imagery acquisition, temporal stacks generation, land cover displacement (i.e. dunes movement) analysis, and results synthesis. Preliminary applications show a marginal agreement of the results with the dominant wind direction in the areas which is here considered as the principal driving force of dune dynamics (Lancaster, 2013). However, a number of parameters have to be defined for running the procedure. These likely have a critical influence on the outputs and require an extensive tuning that is not carried out within this preliminary work, leaving large room for improvements and precise directions for the future developments of this research.
The paper continues as follows. Section 2 describes in details the study areas and the satellite imagery collection and preprocessing steps. In Section 3, the proposed land cover displacements analysis is presented. A preliminary application of the procedure is presented discussed in Section 4. Conclusions and planned future improvements of the work are outlined in Section 5.

Study areas and data procurement
The selected study areas are located in the Central United Arab Emirates (UAE) desert (∼ 6000 Km 2 ) and in the Central Egypt desert (∼ 2000 Km 2 ), close to the Kharga Oasis (Hassoup, 2019) (Figure 1). The UAE study area is mostly characterized by wide linear crescentic Quaternary dunes and Quaternary eolian sand sheet composed by continental sands and gravels, whereas the Egypt study area includes widespread development of barchans dunes. Additional information on the geomorphology of the two deserts can be found e.g. in (Lancaster, 2013).
The proposed procedure exploits images repetitively acquired by a multispectral satellite camera over the same area at different times. Considered images in this work are both Landsat 8 (11 spectral bands up to 30 m resolution) and Sentinel-2 (13 spectral bands up to 10 m resolution) images acquired in the periods May 2013 -November 2015 (Landsat 8) and December 2015 -December 2018 (Sentinel-2). The choice of the period reflected the availability of medium to high-resolution open imagery according to the first acquisition time of the two satellites over the study areas. This allowed investigating the longest series of dunes movement -at the starting time of the presented analysis -as well as the performances on different resolution images. This is discussed in Sections 4.1. The actual average acquisition time over the two areas was approximately 16 days for the Landsat 8 and 7 days for the Sentinel-2. This allowed for retrieving -at least -one full scene per month covering both areas with limited or absent cloud cover by selecting among the available tiles in the imagery catalogues. Data were downloaded according to the granules overlapping the study areas ( Figure 1) using the open data portals of the providers, i.e. the USGS EarthExplorer for the Landsat 8 imagery (USGS, 2020) and the Copernicus Open Access Hub for the Sentinel-2 imagery (Copernicus, 2020). Landsat 8 images were available from the catalogue as Level-2 data (surface reflectance). Sentinel-2 catalogue offered instead only Level-1C data (top of atmosphere reflectance) for both areas in the considered time period.

Temporal stack generation
Critical to any change detection analysis using multitemporal and multispectral satellite imagery is the stacking operation that is the generation of consistent radiometric and spatially coregistered images on which running the analysis. This is done to ensure pixels having equal spectral information and representing always the same land portion within a multitemporal stack (Tewkesbury et al., 2015). The stacking operation for Landsat 8 involved only spectral bands at 30 m resolution while for Sentinel-2 images, spectral bands at 10 m and 20 m meters resolution were considered through upsampling of the 20 m bands to 10 m. Due to the availability of only Sentinel-2 Level-1 data for the considered areas and time period, a Dark Object Subtraction (DOS) (Chavez Jr, 1988) correction was preliminarily applied to these bands to enhance the imagery comparability in the multitemporal stack by adjusting atmospheric noises. The selected bands were then clipped over the study areas extends and recombined into a series of multipsectral images thus obtain the multitemporal stack. Two stack were produced with a monthly time step respectively for the Landsat 8 and the Sentinel-2 data.
The above operations were fully achieved by coupling GRASS GIS (GRASS Development Team, 2018) with custom Python scripts. The use of GRASS GIS was preferred among other available GIS software for the following reasons. a) the images of a stack produced in GRASS GIS are consistently coregistered in space in an automatic way. This because of the standard routine of the software that constraints raster data processing operations to a fixed computational region. b) the software is fully equipped with optimized input/output operations for raster data that facilitate the manipulation of large data volumes. c) the extended support to Python programming allows for developing of complex processing pipelines allowing to automatize recursive operations on data.

LAND COVER DISPLACEMENTS ANALYSIS
The proposed land cover displacements analysis for the investigation of dunes movement considers iteratively a couple of scenes in the stack acquired at two different times -t (master) and time t+∆t (slave) -over the same geographical area. A set of pre-processed scenes are sorted according to a time interval and a time step at which dunes movement is intended to be assessed. In this study, we considered two different master/slave relationship that we called fixed and moving master. In the fixed master setting, shifts are computed with respect to the initial scene of the time interval in the input stack such as t, t+t∆t; t, t+2∆t; etc. In the moving master setting, shift are computed by using a different master and slave couple at each comparison using sequential scenes such as t, t+∆t; t+t∆t, t+2∆t; etc. The full workflow of the proposed procedure is outlined in the following (Figure 2).
According to the adopted master/slave relationship, the band histograms of the slave scene are transformed to match the ones of the master scene ( Scenes are then classified to extract binary raster masks distinguishing target sand pixels from the other background land cover classes (e.g. build-up, crops, etc.), that can be excluded from the analysis of displacements. An unsupervised procedure based on the K-Mean clustering coupled with the Maximum Likelihood Classifier is used for the classification (Duda, Canty, 2002). Masks are derived from the classified scenes by means of manual reclassification (Figure 2b). The unsupervised classification was preferred to a supervised procedure for this early study due to the long digitalization time to generate multiple training samples, that would be necessary for classifying a large number of scenes. However, the higher accuracy that is generally provided by supervised classification algorithms (Lu, Weng, 2007) may be critical to improve the land cover displacements analysis performance. Considerations on the above are reported in Section 4.1.
The movement detection considers two masks derived from scenes acquired at time t (master) and time t+∆t (slave) over the same geographical area. The shifts between mask pixels are estimated by means of local cross-correlation analysis (You et al., 2017) applied iteratively to each pair of masks in the multitemporal stack. Shifts are detected for each pixel by using a moving window (i.e. a subset of neighbouring pixels). Two homologous windows of equal shape are defined respectively on the master and the slave masks. At each iteration, the slave scene subset is shifted in all directions around its central pixel and, for each shift, a cross-correlation coefficient (ρ) is computed with respect to the master scene subset (Figure 2c). The features of the shift (i.e. directions and displacement) that maximize ρ are stored into output raster maps -which have equal extension and resolution of the input scenes -in the position of the moving windows central pixel on the master scene. The moving window is then centred on another pixel and the procedure is repeated such as information on the most probable shift is computed for each pixel in the scene.
Outputs of the local cross-correlation analysis mainly consist of three raster maps for each time step (Figure 2d). These include a map with shift directions (degree from the North), a map with displacements (measured in pixels unit), and a map of root mean square cross-correlation errors (from 0 to 1) (Fienup, 1997) that are metrics of the reliability of each detected shift. A temporal series of directions, displacements and errors maps is thus obtained from the procedure applied to the whole input stack.
The cross-correlation analysis is performed through the dedicated functions of the scikit-image Python library coupled with GRASS GIS. This kind of analysis is popular in digital image processing while it has been seldom applied in the literature -to the extent of the authors' knowledge -on land cover displacement using satellite imagery. Cross-correlation of imagery was originally proposed to track movements of the atmosphere, oceans, and glacier fringes (You et al., 2017) that leaves large room for future research on its application to inland land cover displacements. Figure 2. Schematic of the proposed workflow for the desert dunes movement analysis using multispectral data

Results synthesis
The final purpose of the current study is to investigate dunes movement patterns at a global spatial and temporal scale, that means for the whole study area along the considered study period. To that end, a simple visual comparison of maps aided by GIS software may not be viable due to the amount of information (i.e. pixel-wise displacement features) to be considered. In view of the above, the use of summarizing graphs is here suggested to achieve assessments on dunes movement. In particular, the use of windrose diagrams is leveraged in this work to create snapshots of the dominant movement directions and displacements. The windrose diagram allows displaying oriented histograms of cumulative shifts (in percentage) that are read directly from the temporal series of output maps, obtained by the local cross-correlation analysis. Results can be filtered by removing from the diagram all those shifts that present a crosscorrelation error higher than an arbitrary threshold. Therefore, each detected pixel shift on the output maps -showing a crosscorrelation error lower than the selected threshold -is included in the diagram (Figure 3). During the results synthesis, displacements can be multiplied by the scenes pixel resolution to display them in meters. The creation of windrose diagrams was accomplished through the use of the windrose Python library (Roubeyrie, Celles, 2018). The run times required by each computational step are reported in Table 1.

PRELIMINARY APPLICATIONS
Two experiments were performed on a subset of the stacks to qualitative assess both the capability of the procedure in detecting dune movements as well as the relative influence of the parameters set up on the results. Experiments and their settings are summarized in Table 2. Results are included in figures 4 and 5. Results of the experiments are finally compared with dominant wind blow directions registered in the study areas ( Figure  6) which are here considered as the main driving force of the desert sand dunes movement. Results discussion is included in the following.

Results overview and discussion
The proposed procedure requires to set up a number of parameters for running computations and -in turn -implies extensive tuning operations depending on the target land cover features of the considered study area. Preliminary results show a marginal agreement -in terms of directions -between the detected dunes movement and the dominant winds which are here considered as the main driving force of desert sand movements (see figures 4, 5, and 6). The best agreement is found under experiment E1 settings for the Sentinel-2 stacks on both areas. Displacements appear moderately overestimated according to reference data in the literature, see e.g. (Lancaster, 2013) and (Hassoup, 2019) where the estimated net movement rate for dunes in the studies areas varies from dozens of meters (Egypt) to few meters per year (UAE). However, dunes generally show heterogeneous dynamics both in time ad space due to the local influence of wind storms and the terrain morphology (Haijiang et al., 2008). According to the above, Sentinel-2 imagery outperform the Landsat 8, suggesting the pixel resolution being a critical parameter for dunes movement applications. Despite the pixel resolution, the larger moving window of experiment E1 than experiment E2, combined with a longer time step and the use of a fixed master scene, procures more clustered results -in terms movement directions -to the dominant winds blow directions.
In considering the above, the adoption of long time steps (e.g. one year) on a fixed master scene can limit the detection of seasonal or periodical displacements (Hassoup, 2019) which may mask baseline or principal dunes movement patterns. On the other hand, the use of shorter time steps and a moving master should be preferred when the dynamic of periodical displacements is the focus of the analysis. Second, the size selection for the moving window should reflect the spatial configuration of the land cover phenomena under investigations (e.g. average dune extension, maximum expected dunes shift, etc.) in order to improve results (You et al., 2017). To that end, the selection of the error threshold should be also supported by sensitivity analyses. Furthermore, the classification step is likely to have a strong impact on the cross-correlation analysis. In this study, binary masks depicting dune and not-dune patterns were employed to speed up cross-correlation computations on wide scenes. However, the same analysis can be applied to land cover maps with multiple classes as well as on a combination of the original image bands. . It is intuitive that the more accurate is the classification the better is the detection of movements which suggests -in turn -the need for a proper selection of classification procedure providing the best compromise between accuracy and applicability on large images stack.
Additional sources of error can also affect the results. In particular, the satellite imagery georeferencing offered by the data providers has a Root Mean Square Error (RMSE) for the positional accuracy that ranges between 10 m (RMSE Sentinel-2) and 12 m (RMSE Landsat 8). These errors may manifest as uncontrolled shifts among images that are not due to the actual dynamics of the dunes then biasing the movement detection.
In addition, the minimum displacement that can be detected is equal to or larger than the pixel resolution. Therefore, movements of a few meters that likely occur by considering short time steps cannot be sensed by the procedure because of the native resolution of input data, to which the minimum mappable unit is limited. This is particularly true for the Landsat 8 imagery.
Finally, validation of results by means of comparison with ground truth information -intended as field data or remote sensing products with higher resolution and accuracy than the target satellite imagery -is advised. Suitable validation data sources are represented by e.g. Very High-Resolution (VHR) multitemporal imagery, Laser Imaging Detection and Ranging (LIDAR), and derived Digital Surface Models (DSM).
The issues discussed above were not directly tackled in this early stage of work. Nonetheless, these provide an exhaustive list of directions for the future developments of the proposed procedure.

CONCLUSIONS
In the current paper, an experimental and semi-automatic procedure for desert sand dunes movement analysis by means of multispectral and multitemporal imagery was outlined. The computational steps were described in details, according to both data and software requirements. Preliminary application on two study areas allowed obtaining insights into desert sand dunes dynamics. Preliminary results showed marginal agreements in terms of movement directions. Concerning displacements, additional analyses are instead required to produce better quantitative results.
In view of the above, the need for an extensive review of the procedure was remarked. Comprehensive disclosure of leading improvements was provided. These include additional tests on best parameters configuration (i.e. time step, master/slave relationship, moving window size, and error threshold), on the selection of best performing classification algorithm as well as on the refinement of pre-processing operations for the selected satellite imagery, by accounting in particular for the spatial and temporal traits of the phenomenon under investigation. Indeed, the future developments of this work are planned to tackle all edges of the above, coupled with a systematic review of the raw Python code produced for the preliminary tests. To that end, the procedure leverages exclusively FOSS GIS software and open satellite imagery. This is of paramount importance to provide the analysis with a potential to be empowered, replicated and improved.
Although being experimental, the procedure promises to produce valuable results for investigating sand dune movements. Resulting maps and summary graphs can be readily employed into sand mitigation planning by supplying baseline information on dunes dynamics in a timely and cost-effective manner.
The suitability for applications of the modern high-resolution open satellite imagery also in the private sector -such as in transport infrastructure -opens new opportunities to unpin both the economic and scientific potential of this data by adding value to the open data policies while enforcing and motivating investments in satellite Earth Observations.