Digital Image Correlation of Google Earth Images for Earth’s Surface Displacement Estimation

: An increasing number of satellite platforms provide daily images of the Earth’s surface that can be used in quantitative monitoring applications. However, their cost and the need for speciﬁc processing software make such products not often suitable for rapid mapping and deformation tracking. Google Earth images have been used in a number of mapping applications and, due to their free and rapid accessibility, they have contributed to partially overcome this issue. However, their potential in Earth’s surface displacement tracking has not yet been explored. In this paper, that aspect is analyzed providing a speciﬁc procedure and related MATLAB ™ code to derive displacement ﬁeld maps using digital image correlation of successive Google Earth images. The suitability of the procedure and the potential of such images are demonstrated here through their application to two relevant case histories, namely the Slumgullion landslide in Colorado and the Miage debris-covered glacier in Italy. Result validation suggests the e ﬀ ectiveness of the proposed procedure in deriving Earth’s surface displacement data from Google Earth images.


Introduction
In the last decades, an increasing number of satellite platforms are persistently acquiring Earth observation data in the form of optical, multispectral, hyperspectral, and radar products (e.g., GeoEye, WorldView, Sentinel, Jilin, Radarsat, Envisat). Such products form a very important database for Earth surface control and mapping, offering a general medium to high resolution (0.3 to >50 m single-sided pixel dimension) and a revisiting time of a few days to a few weeks. The highest resolution products (<10 m single-sided pixel dimension) are often used for Earth's surface deformation monitoring. Landslides, glaciers, earthquakes, and dune-field displacement tracking are examples of applications of satellite products that, covering small to large areas, have the potential to provide an exhaustive overview of the process [1][2][3][4].
Among the most effective techniques that can be used to monitor Earth's surface displacement, the increasingly used interferometry technique called Differential Interferometric Synthetic Aperture Radar (DInSAR, [5,6]) has the potential to support low-rate deformation tracking applications and mapping [7][8][9][10]. Moreover, due to the characteristics of the sensors, it is able to work in all-weather and day and night conditions. However, this approach is not always suitable due to limitations related to the maximum detectable velocity [11] and potential high cost, especially in the context of high-resolution analyses (e.g., persistent scatterer approaches, [12,13]). This limitation has been partially overcome through the implementation of specific low-resolution approaches [14][15][16]. The use

Slumgullion Landslide
The Slumgullion landslide is located in the San Juan Mountains of southwestern Colorado, USA, and consists in a 3.9 km long, slow-persistently moving, active landslide associated with a larger and dormant landslide deposit (Figure 1a; [30,32]). The active landslide ranges in elevation from 2950 to 3650 m asl and the slope of the ground surface is of approximately 8 • . Total landslide volume has been estimated as~20 × 10 6 m 3 and the average thickness is around 13 m [29]. The landslide has a typical flow-like surface morphology and therefore, although the movement occurs mainly by sliding along discrete basal and lateral slip surfaces, it has been described as an earthflow [29]. The Slumgullion landslide has been actively moving since at least 300 years ago. Landslide velocity ranges between approximately 0.15 and 7 m/year [30]. Higher velocities have been observed along the longitudinal center of the middle sector of the landslide and lower velocities have been observed at the head and the toe. Velocity variations in space materialize distinct kinematic elements bounded by deformational structures [30]. The landslide moves more slowly during winter when most precipitation falls as snow and faster during spring due to snowmelt and rainfall events [30]. The authors of [33] indicate that, over the next century, the landslide movement rate is expected to decrease due to a change in regional Remote Sens. 2020, 12, 3518 3 of 14 moisture balance conditions. This seems to be consistent with estimations completed by the authors of [16,34], indicating a maximum velocity from 4 to 4.4 m/year in the period 2011-2013.
Remote Sens. 2020, 12, x FOR PEER REVIEW 3 of 14 estimations completed by the authors of [16,34], indicating a maximum velocity from 4 to 4.4 m/year in the period 2011-2013.

Miage Debris-Covered Glacier
The Miage Glacier is in the eastern Alps along the southwestern side of Mont Blanc in the Veny valley of northern Italy (Figure 1b). It consists of a 6 km debris-covered glacier with a total surface of approximately 11 km 2 [31]. The debris cover extends from an elevation of approximately 2400 m asl, downslope, toward the frontal lobes forming the lower end of the ablation tongue. Grain size of debris forming the glacier cover ranges from rock boulders to fine pebbles and sand. Debris cover thickness ranges between a few cm of the upper sector to 1.5 m of the lower sector and is consistently

Miage Debris-Covered Glacier
The Miage Glacier is in the eastern Alps along the southwestern side of Mont Blanc in the Veny valley of northern Italy (Figure 1b). It consists of a 6 km debris-covered glacier with a total surface of approximately 11 km 2 [31]. The debris cover extends from an elevation of approximately 2400 m asl, downslope, toward the frontal lobes forming the lower end of the ablation tongue. Grain size of debris forming the glacier cover ranges from rock boulders to fine pebbles and sand. Debris cover thickness ranges between a few cm of the upper sector to 1.5 m of the lower sector and is consistently related to surface slope. Its presence is sustained by continuous rockfalls and rock avalanches [31] and its movement is consistently related to glacier-tongue seasonal displacement. The upper and middle sectors of the tongue move with an estimated average velocity of approximately 60 m/year [31,35]. In the lower sector of the ablation tongue, velocity decreases to an average of 45 m/year with further notable reduction toward the lower end of the glacier tongue. In this zone, velocity decreases to less than 4 m/year and reaches approximately 20 m/year next to the beginning of the southern lobe [31].

Materials and Methods
For estimating the displacement of both the medium and lower sectors of the active part of the Slumgullion landslide and the medium and lower sectors of the Miage debris-covered glacier, color Google Earth images saved as 4800 × 2674 jpeg files were used. Images depicting the active part of the Slumgullion landslide were acquired on 6 June 2013 and 14 October 2015. Images representative of the Miage glacier valley were acquired on 29 August 2009 and 28 May 2011. In both cases, following the master/slave logic of the digital image correlation technique, the first image was used as master and the second one as slave.
The procedure used for image selection, acquisition, and processing is reported in Figure 2. Prior to image selection, the study area (for each case study) was identified into a GIS environment and a rectangular polygon with a dimension of 2400 × 1337 m (H/L ratio~0.55) was created and exported into a shp format. The dimension of the polygon was arbitrarily estimated considering the maximum available resolution of Google Earth images and the need for an acceptable final single-sided pixel dimension (in this case around 0.5 m). Corner coordinates of the polygon were exported in a separated text file for subsequent use in image rectification. Once imported into Google Earth, the polygon was used as reference for image selection. Specifically, after fitting the polygon to the screen by adjusting both image position and scale, acquisition dates of interest were selected, and the corresponding color images were acquired ( Figure 3). In this case, acquisition dates were selected as a function of image sharpness. The images were saved with the reference polygon draped on them for subsequent processing. Specifically, after importing images into a GIS environment, they were digitally rectified considering four pairs of corner coordinates corresponding to the polygon vertex and affine method [36]. This procedure consists in the georeferentiation of the images through digital rectification, allowing the improvement of master and slave image co-registration. The error related to image rectification (i.e., co-registration) multiplied by 2 was considered as the accuracy of detected image feature position and was considered in displacement field map restitution. Once rectified, images were exported as geotiff selecting two different resolutions (i.e., pixel dimension). A first set (i.e., master and slave) was exported selecting full resolution (i.e., 0.5 m single-sided pixel dimension) and a second one (master image only) selecting the prospective resolution of the displacement grid used for final displacement mapping (i.e., 50 m single-sided pixel dimension). The first, high-resolution, set was directly used in tracking analysis. The second set was simply used for subsequent mapping, and its resolution (50 m was selected for this analysis) can be arbitrarily selected as a function of specific needs.
Once generated, images were processed using a specifically developed two-step procedure, coded in MATLAB™ (R2019b, © 1994-2020 The MathWorks, Inc., Natickc, MA, USA) (see Figure 2). Such a procedure, based on an eigen feature identification approach (i.e., corner points; [37]), uses the Kanade-Lucas-Tomasi (KLT) feature-tracking algorithm [27,28] which is applied to greyscale-rectified images to estimate the displacement of corner points in the time interval between successive image acquisitions. The KLT method allows to obtain displacement in the form of tracked feature displacement vectors, displacement components, 2D displacement distribution graphs, and displacement field maps. Algorithm parameterization was consistently completed using a trial and error procedure aimed at maximizing the number of trackable corner points and reducing the displacement estimation error. Minimum accepted quality of corners in the image, in the form of a fraction of the maximum corner metric value, was considered equal to 0.01. A maximum forward-backward bidirectional error threshold of 0.05 pixels was adopted to eliminate points that could not be reliably tracked. Additionally, since the KLT-based point tracker uses image pyramids, where each level is reduced in resolution by a factor of two compared to the previous level, and considering that the number of pyramid levels is Remote Sens. 2020, 12, 3518 5 of 14 related to the magnitude of the displacement to handle, seven pyramid levels were considered for the analysis.
Remote Sens. 2020, 12, x FOR PEER REVIEW 5 of 14 be reliably tracked. Additionally, since the KLT-based point tracker uses image pyramids, where each level is reduced in resolution by a factor of two compared to the previous level, and considering that the number of pyramid levels is related to the magnitude of the displacement to handle, seven pyramid levels were considered for the analysis. In the first step (i.e., Image acquisition and pre-processing and Image processing and deviation analysis of Figure 2), the procedure accounts for the (i) preparation of a reference mapping grid through single ultra-low-resolution image (ULRis; i.e., 50 m single-sided pixel dimension), (ii) estimation of the full-resolution images' pixel dimension, and (iii) definition of the displacement deviation (i.e., local error) through a controlled displacement analysis in known non-moving zones In the first step (i.e., Image acquisition and pre-processing and Image processing and deviation analysis of Figure 2), the procedure accounts for the (i) preparation of a reference mapping grid through single ultra-low-resolution image (ULRis; i.e., 50 m single-sided pixel dimension), (ii) estimation of the full-resolution images' pixel dimension, and (iii) definition of the displacement deviation (i.e., local error) through a controlled displacement analysis in known non-moving zones (see Figure 3). The reference mapping grid (i) is derived from the ULRis by emulating the image grid through vector meshing. Pixel dimension (ii) is consistently estimated from image-pixel coordinates and is stored in a specific variable that forms the basis for subsequent displacement estimation. Non-moving zones (iii) are selected by the user, creating a mask through the roipoly interactive tool implemented in MATLAB™. The code can be arranged with a selectable number of non-moving zones. In this case, since both image sets are characterized by a central moving zone, because of the presence of landslide and glacier, the code is arranged with two subsequent steps of non-moving zone identification. Once identified, the non-moving zone masks are used to support displacement analysis for local deviation estimation. This deviation is related to change in scale due to the potential presence of reliefs and image distortion that has not been successfully corrected by digital image rectification. Deviation estimation is completed using the corner point tracking approach described above and provides both X and Y components of deviation in the form of distribution graphs and raster maps. Such maps, derived by the spatial interpolation of data over the reference mapping grid, form the basis for the deviation compensation processing. Data interpolation is consistently completed using the natural neighbor method [38].
presence of reliefs and image distortion that has not been successfully corrected by digital image rectification. Deviation estimation is completed using the corner point tracking approach described above and provides both X and Y components of deviation in the form of distribution graphs and raster maps. Such maps, derived by the spatial interpolation of data over the reference mapping grid, form the basis for the deviation compensation processing. Data interpolation is consistently completed using the natural neighbor method [38].
In the second step (i.e., Displacement analysis of Figure 2), the displacement analysis is completed without mask constraints, so that is extended to the whole study area depicted in the images. The KLT algorithm tracks identified corner points providing X and Y coordinates in the master and slave images. Such coordinates form the basis for displacement component estimation. Once estimated, displacement components are used to derive the displacement module and its direction through simple geometrical relations. The code provides a distribution graph of displacement. Subsequently, displacement components estimated at specific positions are used to derive X and Y displacement rasters through the spatial interpolation of data over the reference mapping grid. Such rasters are used to derive the displacement field map. Before deriving the displacement field map, X and Y component deviation rasters are subtracted from the X and Y displacement raster to remove the error related to residual image distortion. The final displacement field map shows the displacement distribution across the selected study area following a Lagrangian description scheme. In the second step (i.e., Displacement analysis of Figure 2), the displacement analysis is completed without mask constraints, so that is extended to the whole study area depicted in the images. The KLT algorithm tracks identified corner points providing X and Y coordinates in the master and slave images. Such coordinates form the basis for displacement component estimation. Once estimated, displacement components are used to derive the displacement module and its direction through simple geometrical relations. The code provides a distribution graph of displacement. Subsequently, displacement components estimated at specific positions are used to derive X and Y displacement rasters through the spatial interpolation of data over the reference mapping grid. Such rasters are used to derive the displacement field map. Before deriving the displacement field map, X and Y component deviation rasters are subtracted from the X and Y displacement raster to remove the error related to Remote Sens. 2020, 12, 3518 7 of 14 residual image distortion. The final displacement field map shows the displacement distribution across the selected study area following a Lagrangian description scheme.
Once the displacement field is obtained, this is validated using a comparison with available local velocity data and considering the expected earthflow and glacier spatial kinematics [16,30,[39][40][41].

Slumgullion Landslide
The analysis indicates that the middle and lower sectors of the active part of the Slumgullion landslide consistently moved between 6 June 2013 and 14 October 2015 (860 days). A total of 7069 eigen features were successfully identified and tracked. Their distribution is shown in Figure 4a. As observed, the density of recognized features in the moving zone is much higher than in the non-moving zone and no outliers (displacement vectors related to inconsistently tracked objects) are visually recognizable. Image rectification, completed using the coordinates of image corners, indicates a rectification RMSE of~1 m.  Once the displacement field is obtained, this is validated using a comparison with available local velocity data and considering the expected earthflow and glacier spatial kinematics [16,30,[39][40][41].

Slumgullion Landslide
The analysis indicates that the middle and lower sectors of the active part of the Slumgullion landslide consistently moved between 6 June 2013 and 14 October 2015 (860 days). A total of 7069 eigen features were successfully identified and tracked. Their distribution is shown in Figure 4a. As observed, the density of recognized features in the moving zone is much higher than in the nonmoving zone and no outliers (displacement vectors related to inconsistently tracked objects) are visually recognizable. Image rectification, completed using the coordinates of image corners, indicates a rectification RMSE of ~1 m.  Deviation analysis, conducted considering the NW and SE non-moving zones separated by the moving landslide, indicates for the X and Y displacement components a deviation ranging between −4 and +5 m, and −7 and +1 m, respectively (Figure 5a). The modes of distributions are, in both cases, at −1 m. X and Y deviation spatial distribution is reported in Figure 6a,b. Displacement component distribution as well as resultant 2D displacement are reported in the graph of Figure 5b, and displacement distribution is reported in the displacement field map of Figure 4b. As indicated in the graph, overall, a displacement ranging from −11 to + 5 m was registered along the X-axis of the image, while a displacement ranging from −7 to + 9 m was registered along the Y-axis. 2D resultant displacement ranges between 1 and 12.2 m. Considering the accuracy of positioning estimated as the RMSE multiplied by 2, that in this case is estimated as 2 m, 2D resultant displacement ranges between 0 and 10.2 to 14.2 m (i.e., 12.2 ± 2 m) for a maximum velocity ranging between~4 and~6 m/year. Deviation analysis, conducted considering the NW and SE non-moving zones separated by the moving landslide, indicates for the X and Y displacement components a deviation ranging between −4 and +5 m, and −7 and +1 m, respectively (Figure 5a). The modes of distributions are, in both cases, at −1 m. X and Y deviation spatial distribution is reported in Figure 6a,b. Displacement component distribution as well as resultant 2D displacement are reported in the graph of Figure 5b, and displacement distribution is reported in the displacement field map of Figure 4b. As indicated in the graph, overall, a displacement ranging from −11 to + 5 m was registered along the X-axis of the image, while a displacement ranging from −7 to + 9 m was registered along the Y-axis. 2D resultant displacement ranges between 1 and 12.2 m. Considering the accuracy of positioning estimated as the RMSE multiplied by 2, that in this case is estimated as 2 m, 2D resultant displacement ranges between 0 and 10.2 to 14.2 m (i.e., 12.2 ± 2 m) for a maximum velocity ranging between ~4 and ~6 m/year.

Miage Debris-Covered Glacier
The analysis indicates that the middle and lower sectors of the Miage debris-covered glacier consistently moved between 29 August 2009 and 28 May 2011 (637 days). A total of 19,067 eigen features were successfully identified and tracked. Their distribution is shown in Figure 7a. As observed, the density of recognized feature in the moving zone is much higher than in the nonmoving zone, although it is not fully covered. No outliers (displacement vectors related to inconsistently tracked objects) are visually recognizable. Image rectification, completed using the coordinates of image corners, indicates a rectification RMSE of ~1.5 m.  moving landslide, indicates for the X and Y displacement components a deviation ranging between −4 and +5 m, and −7 and +1 m, respectively (Figure 5a). The modes of distributions are, in both cases, at −1 m. X and Y deviation spatial distribution is reported in Figure 6a,b. Displacement component distribution as well as resultant 2D displacement are reported in the graph of Figure 5b, and displacement distribution is reported in the displacement field map of Figure 4b. As indicated in the graph, overall, a displacement ranging from −11 to + 5 m was registered along the X-axis of the image, while a displacement ranging from −7 to + 9 m was registered along the Y-axis. 2D resultant displacement ranges between 1 and 12.2 m. Considering the accuracy of positioning estimated as the RMSE multiplied by 2, that in this case is estimated as 2 m, 2D resultant displacement ranges between 0 and 10.2 to 14.2 m (i.e., 12.2 ± 2 m) for a maximum velocity ranging between ~4 and ~6 m/year.

Miage Debris-Covered Glacier
The analysis indicates that the middle and lower sectors of the Miage debris-covered glacier consistently moved between 29 August 2009 and 28 May 2011 (637 days). A total of 19,067 eigen features were successfully identified and tracked. Their distribution is shown in Figure 7a. As observed, the density of recognized feature in the moving zone is much higher than in the nonmoving zone, although it is not fully covered. No outliers (displacement vectors related to inconsistently tracked objects) are visually recognizable. Image rectification, completed using the coordinates of image corners, indicates a rectification RMSE of ~1.5 m.

Miage Debris-Covered Glacier
The analysis indicates that the middle and lower sectors of the Miage debris-covered glacier consistently moved between 29 August 2009 and 28 May 2011 (637 days). A total of 19,067 eigen features were successfully identified and tracked. Their distribution is shown in Figure 7a. As observed, the density of recognized feature in the moving zone is much higher than in the non-moving zone, although it is not fully covered. No outliers (displacement vectors related to inconsistently tracked objects) are visually recognizable. Image rectification, completed using the coordinates of image corners, indicates a rectification RMSE of~1.5 m.
Deviation analysis, conducted considering the NW and SE non-moving zones separated by the moving glacier, indicates for X and Y displacement components a deviation ranging between −5 and +8 m, and −5 and +10 m, respectively (Figure 8a). The modes of distributions are at +2 and +4 m, respectively. X and Y deviation spatial distribution is reported in Figure 9a,b. Displacement component distribution as well as resultant 2D displacement are reported in the graph of Figure 8b, and displacement distribution is reported in the displacement field map of Figure 7b. As indicated in the graph, overall, a displacement ranging from −6 to +34 m was registered along the X-axis of the image, while a displacement ranging from −20 to + 43 m was registered along the Y-axis. 2D resultant displacement ranges between 0.3 and 53 m. Considering the accuracy of positioning estimated as the RMSE multiplied by 2, that in this case is estimated as 3 m, 2D resultant displacement ranges between 0 and 59 m for a maximum velocity of~33 m/year. Remote Sens. 2020, 12, x FOR PEER REVIEW 9 of 14 Deviation analysis, conducted considering the NW and SE non-moving zones separated by the moving glacier, indicates for X and Y displacement components a deviation ranging between −5 and +8 m, and −5 and +10 m, respectively (Figure 8a). The modes of distributions are at +2 and +4 m, respectively. X and Y deviation spatial distribution is reported in Figure 9a,b. Displacement component distribution as well as resultant 2D displacement are reported in the graph of Figure 8b, and displacement distribution is reported in the displacement field map of Figure 7b. As indicated in the graph, overall, a displacement ranging from −6 to +34 m was registered along the X-axis of the image, while a displacement ranging from −20 to + 43 m was registered along the Y-axis. 2D resultant displacement ranges between 0.3 and 53 m. Considering the accuracy of positioning estimated as the RMSE multiplied by 2, that in this case is estimated as 3 m, 2D resultant displacement ranges between 0 and 59 m for a maximum velocity of ~33 m/year.

Discussion and Validation
The growing archives of satellite optical remote sensing images and the advent of digital image correlation algorithms offer great potential for the detection and monitoring of landslides and glaciers purpose of this research was to present a new low-cost and easily accessible method to perform image correlation types of studies, based on the use of freely available Google Earth imagery. The two case studies used as test sites were the Slumgullion landslide and the Miage debris-covered glacier. Results of the analysis indicate that both the Slumgullion landslide and the Miage debris-covered glacier actively moved in the reference period. This is consistent with previous analyses indicating that, in both cases, long-term movement has been observed and is expected for the near future [30,31]. In both cases, a typical displacement field for both flow-type landslide and glacier was observed. The Slumgullion landslide moved with a maximum velocity, registered at the upper part of its middle sector, slightly below the landslide neck, around 5 m/year. The minimum velocity was registered at the toe with a rate around 1 m/year. These velocities are consistent with those calculated by the authors of [34] and predicted by the authors of [33]. Conversely, they appear to be much lower than those predicted by the authors of [30] for the period 1998-2002. This is consistent with the longterm prediction of the displacement rate completed by the authors of [33]. The displacement field of Figure 4b depicts a typical condition for this kind of landslide with velocity variations related to the width of the landslide corridor and to the distance from the toe. Specifically, as the width of the landslide corridor decreases, velocity increase and vice versa [40]. In addition, the proposed coded procedure seems to be able to recognize (lateral) spreading occurring at the toe of the slide and velocity reduction at cross-section edges.  Figure 7b). This difference may be related to the ongoing reduction of glacier velocity connected to a deficit in mass balance induced by global warming, commonly observed in the Alpine region of Switzerland, France, and Italy [45][46][47]. A similar difference was also observable next to the beginning of the southern lobe where a velocity of 12 m/year was estimated.
In both cases, local inconsistent displacement vectors are present. Although they are mostly located along image edges, which are the most critical parts of images in object tracking perspective (see Figure 6a,b or Figure 9a,b), several occurrences can be recognized within both study regions. Those located along the edges of the images can be interpreted as local artefacts or false positives connected to the lack of an efficient tracking of non-moving points for generalized distortion compensation. Conversely, those located within the study areas may be representative of local error connected to the interpolation method. In addition, they can be considered as the indication of an erroneous interpretation of those areas as non-moving zones. This underlines a potential limit of the proposed procedure indicating that, in the absence of a detailed knowledge of the study area, results from the analysis have to be carefully interpreted and validated.

Conclusions
The analysis indicates that the proposed procedure and related MATLAB™ code are able to provide quantitative displacement data, in the form of displacement field maps and graphs, from the digital image correlation of Google Earth images. Such images require a specific acquisition scheme and need to be pre-processed to improve image co-registration and reduce internal distortion related to local-scale variation. The application to the Slumgullion landslide in Colorado (USA) and the Miage debris-covered glacier in Italy confirmed the potential of the procedure to derive a consistent estimation of displacement distribution and, at the same time, underlines its potential limit related to

Discussion and Validation
The growing archives of satellite optical remote sensing images and the advent of digital image correlation algorithms offer great potential for the detection and monitoring of landslides and glaciers [18,[42][43][44]. In [18], the authors showed the use of Pléiades panchromatic images for the analysis of a landslide located in the South French Alps. They proposed the use of a multiple pairwise image correlation (MPIC) technique (to obtain horizontal displacement fields) and different multi-temporal indicators for a more accurate detection and quantification of surface deformations. Diversely from the research presented in [18], which aimed to improve the accuracy of image correlation techniques, the authors of [43] implemented a fuzzy Hough transform to analyze the GoLIVE dataset (https://nsidc.org/data/golive) with the goal of measuring seasonal velocity changes in glaciers, enhancing the extraction of the timing and location of ice flow events. A quantitative assessment of digital image correlation methods has been recently presented by the authors of [42], whose work demonstrated that all algorithms are capable of quantifying slope instability displacements, with average errors ranging from 8% to 12% of the observed maximum displacement.
In light of these elements of evidence, it is clear that image correlation techniques represent a powerful tool in the study of landslide and glacier deformation. However, the cost of the availability of satellite imagery still represents one of the main limitations of these techniques. In this context, the purpose of this research was to present a new low-cost and easily accessible method to perform image correlation types of studies, based on the use of freely available Google Earth imagery. The two case studies used as test sites were the Slumgullion landslide and the Miage debris-covered glacier.
Results of the analysis indicate that both the Slumgullion landslide and the Miage debris-covered glacier actively moved in the reference period. This is consistent with previous analyses indicating that, in both cases, long-term movement has been observed and is expected for the near future [30,31]. In both cases, a typical displacement field for both flow-type landslide and glacier was observed.
The Slumgullion landslide moved with a maximum velocity, registered at the upper part of its middle sector, slightly below the landslide neck, around 5 m/year. The minimum velocity was registered at the toe with a rate around 1 m/year. These velocities are consistent with those calculated by the authors of [34] and predicted by the authors of [33]. Conversely, they appear to be much lower than those predicted by the authors of [30] for the period 1998-2002. This is consistent with the long-term prediction of the displacement rate completed by the authors of [33]. The displacement field of Figure 4b depicts a typical condition for this kind of landslide with velocity variations related to the width of the landslide corridor and to the distance from the toe. Specifically, as the width of the landslide corridor decreases, velocity increase and vice versa [40]. In addition, the proposed coded procedure seems to be able to recognize (lateral) spreading occurring at the toe of the slide and velocity reduction at cross-section edges.
The Miage glacier moved with a maximum velocity of 33 m/year between August 2009 and May 2011. This velocity was registered within the upper part of the lower sector of the ablation tongue where an average velocity of 45 m/year was registered by the authors of [31] between 1975 and 1991 ( Figure 7b). This difference may be related to the ongoing reduction of glacier velocity connected to a deficit in mass balance induced by global warming, commonly observed in the Alpine region of Switzerland, France, and Italy [45][46][47]. A similar difference was also observable next to the beginning of the southern lobe where a velocity of 12 m/year was estimated.
In both cases, local inconsistent displacement vectors are present. Although they are mostly located along image edges, which are the most critical parts of images in object tracking perspective (see Figure 6a,b or Figure 9a,b), several occurrences can be recognized within both study regions. Those located along the edges of the images can be interpreted as local artefacts or false positives connected to the lack of an efficient tracking of non-moving points for generalized distortion compensation. Conversely, those located within the study areas may be representative of local error connected to the interpolation method. In addition, they can be considered as the indication of an erroneous interpretation of those areas as non-moving zones. This underlines a potential limit of the proposed procedure indicating that, in the absence of a detailed knowledge of the study area, results from the analysis have to be carefully interpreted and validated.

Conclusions
The analysis indicates that the proposed procedure and related MATLAB™ code are able to provide quantitative displacement data, in the form of displacement field maps and graphs, from the digital image correlation of Google Earth images. Such images require a specific acquisition scheme and need to be pre-processed to improve image co-registration and reduce internal distortion related to local-scale variation. The application to the Slumgullion landslide in Colorado (USA) and the Miage debris-covered glacier in Italy confirmed the potential of the procedure to derive a consistent estimation of displacement distribution and, at the same time, underlines its potential limit related to the knowledge of the study area as a pre-requisite needed to support result interpretation, especially in the absence of further quantitative data. This particular aspect, connected to the need of tracking displacement in a non-moving zone, is related to the need of reducing internal image deformation that may compromise displacement estimation over the whole area of interest.