Vector-based approach for combining ascending and descending persistent scatterers interferometric point measurements

Abstract The combination of ascending and descending persistent scatterers interferometric (PSI) data by means of resampling and/or spatial interpolation, separately for each synthetic aperture radar geometry, is a commonly followed procedure, limited though by the reduced spatial coverage and the introduced uncertainties from multiple rasterization steps. Herein, an alternative approach is proposed for combining different PSI line-of-sight observables in the vector domain, based on the geographic proximity of PS point targets. In the proposed nearest neighbour vector approach all necessary analysis steps are performed by means of attributes transfer and calculations between features geodatabases, prior to any rasterization. By increasing the number of input point vectors during subsequent interpolation, the overall error budget coming from spatial interpolation is being reduced. The advantages of the proposed vector-based approach compared to the commonly used grid-based procedure are being demonstrated using real data.


Introduction
With the introduction of the multi-baseline synthetic aperture radar (SAR) interferometry concept by Ferretti et al. (2001), various approaches arose for the generation of SAR displacement time series following different processing schemes (Werner et al. 2003;Perissin & Wang 2012). Since then, persistent scatterers interferometric (PSI) techniques are routinely used for research purposes, while their extensive validations in different environments (Agudo et al. 2006; Terrafirma Project) led to more operational frameworks for many geohazard monitoring applications (European Space Agency (ESA) 2013). In the course of their development, novel approaches are presented involving advanced solutions for complex non-linear motion scenarios (Berardino et al. 2002;Mora et al. 2003;Hooper et al. 2004) as well as wide area products (Adam et al. 2013;Duro et al. 2015). Further progress is expected by exploiting the capabilities of missions like the Copernicus Sentinel-1, which offer systematic acquisitions with shorter repeat cycles and large swath coverage (https://sentinel.esa.int). The open and free data policy is an additional asset of the Sentinel missions, promoting scientific research and stimulating innovation.
As imposed by the acquisition geometry of SAR sensors on-board of satellite platforms, having a slanted view (tilted with respect to the vertical direction) referred to as line-of-sight (LOS), only part of the three-dimensional motion filed is captured. To overcome this limitation, the combination of independent 1D SAR Interferometric (InSAR) measurements from ascending and descending satellite orbits is a common adopted practice (Fialko et al. 2001;Gudmundsson & Sigmundsson 2002;Wright et al. 2004;Foumelis et al. 2013), permitting the extraction of the actual Vertical (Up) and E-W motion components from the SAR LOS observables (Hanssen 2001).
PSI results correspond to measurements over point-like scattering targets, typically exported in a common American Standard Code for Information Interchange (ASCII) format, and visualised as vector layers in Geographic Information Systems (GIS) for presentation purposes and further analysis. The combination of poiSnt measurements of different spatial distribution and density in a post-processing frame unavoidably involves resampling to a common regular grid. Resampling (also referred to as gridding) could be limited to cells where original data are available (Tamburini et al. 2010;Tofani et al. 2013), or include spatial interpolation to estimate values over regions of missing data (Hung et al. 2011;Raspini et al. 2012;Dehghani et al. 2013).
In practice, when ascending and descending PSI results are combined following simple resampling, issues arise from the reduced coverage in regions where data from both SAR geometries exist. Moreover, spatial averaging of values from multiple scatters available within each grid cell is also causing smoothing of the displacement field. When interpolation is considered, independent continuous LOS measurement surfaces are obtained, which are then combined using simple raster math operations to derive the Vertical and Easting components. In the referred to case, the lack of PSI data over SAR layover and shadow zones, especially in mountainous regions, as well as due to coherence loss for rapidly decorrelating areas (e.g. vegetated lands), lead to random spatial distributions and data gaps. This adds difficulties to any interpolation algorithm. Additional issues are caused by itself; the inherited disadvantages of the interpolation procedure implemented. Thus, reducing the effects of resampling and interpolation, given that compensating for the SAR acquisition geometry factors is not straightforward, could lead to improved products and possibly widen their applicability.
An alternative procedure is presented herein, using PS data in their initial vector format and by taking advantage of vector manipulation capabilities as well as the geostatistical modules already embedded in many GIS software packages. Although such operations are often adopted in the geospatial science domain, dealing with manipulation and analysis of vector data (Yianilos 1993), they are not sufficiently introduced to medium-resolution geodetic InSAR applications. Several algorithms are though proposed for the geometrical fusion of multi-track three-dimensional (3D) point clouds delivered by PSI or SAR Tomography (TomoSAR) techniques using very high resolution SAR data over urban environments (Gernhardt et al. 2012;Schack & Soergel 2014;Wang & Zhu 2014). Our goal is to reduce the drawbacks of rasterization and interpolation when combining multi-track PSI results, retaining both spatial coverage and statistical properties of the inputs, while working with PS results in their original vector format. Investigation of the performance of various spatial interpolation techniques is not part of the study, as comparable issues in terms of data handling applies.
Nisyros volcano was selected as a pilot site for checking the performance of the proposed approach ( Figure 1). Selection criteria included the limited spatial extent to reduce the computational load, as well as the radial shape of the volcanic island providing the required variations in terms of slope gradients and aspects. These are significant contributing factors for the availability of PS targets due to the inclined viewing geometry of the spaceborne SAR systems.
In the following sections, pre-processing of SAR data to arrive at the required input PSI data-sets is briefly presented, followed by a detailed description of the proposed methodology for combining the ascending and descending PSI measurements. Intermediate investigations within the frame of the aimed analysis on the effect of resampling and spatial interpolation are also presented. This aims to clearly demonstrate the influence of the various post-processing aspects when dealing with such an application. The obtained results are demonstrated in a separated section, whereas subsequent discussion is made on the advantages and limitations of the approach and the planning of steps forward.
Interferometric analysis followed already described processing scheme (Foumelis et al. 2013), for SAR co-registration at sub-pixel accuracy, spectral filtering and interferogram generation using the GAMMA software packages (Wegmüller et al. 1998). Single master, as well as consecutive multi-reference interferograms in full resolution, were formed and used as inputs for the Persistent Scatterers (PS) and the Small Baselines (SB) analysis of the Stanford Method for Persistent Scatterers (StaMPS) software (Hooper 2008). The detection of the persistent scatterers was based on the amplitude dispersion technique (Ferretti et al. 2001). In order to maintain consistency, similar processing parameters and procedure as well as reference area were considered for both data stacks, identifying a total number of 14,123 and 14,890 point targets for ascending and descending geometry, respectively.
The detected PS targets were then geocoded using the Digital Elevation Model (DEM) aided back-geocoding approach implemented in StaMPS. For both co-registration and geocoding to map geometry a 20-m resolution DEM derived from digitised 1/50,000 scale topographic maps was utilised (Hellenic Military Geographical Service (HMGS) 1972). Usage of common heights and geocoding procedure ensures proper spatial overlap, with the geolocation accuracy being at a fraction of the pixel (less than 10 m in horizontal map projection) for both the ascending and descending tracks. Further tuning of PS positioning could be achieved by compensating for potential range and azimuth timing errors as well as for uncertainties of reference points heights per track (Ketelaar et al. 2007). However, the obtained geolocation accuracy for such medium resolution SAR system is considered sufficient for the purpose of the current application demonstration.
Following steps involved PS-SB processing for estimation of displacement histories, presentation of which is out of the scope of the current work. During processing a number of points are rejected based on various quality criteria. The reduced set of PS was considered throughout the subsequent analysis as a more realistic input to the intended investigation. The final reduced PS stacks consist of 13,151 and 12,212 points for the ascending and descending geometry, respectively (Figure 2(a)). As expected, the spatial distribution of the ascending and descending targets, apart from the effect of temporal decorrelation (changes in land cover/use), is highly related to the orientation of slopes, and thus, correlates well with the terrain aspect ( Figure 2(b)). Average densities of 318 and 295 points per km 2 were obtained for the ascending and descending PSI results, respectively (Figure 3(a)-(b)).

Methodology
Starting from a set of two PS point data-sets (ascending and descending), first step involved data import into GIS-specific format ensuring both the correct definition of the projection system and ingestion of the actual PS measurement within the attribute geodatabase.
The approach presented herein is based mainly on tools and functionalities available in most GIS software packages. In our case, ESRI's ArcMap and ArcView software packages (ESRI 2011) were used for the majority of the analysis steps, while external scripts, publicly available (Jenness 2004(Jenness , 2012Tchoukanski 2006) or developed within the frame of the study, were integrated in some parts of the processing.
The analysis followed both commonly applied grid-based approach (from herein referred to as 'GRID') as well as the proposed nearest neighbour vector (NNV) approach. The processing workflow for both approaches is presented graphically in Figure (4).  Processing was repeated for multiple grid spacing, starting from a fine grid of 20-m pixel to a coarse one of 260 m with a step of 20 m, in order to understand and quantify the effect of pixel size selection on the outputs. Details on the interpolation parameters are provided in the following sections. Due to the considerable number of steps involved for final outputs generation and calculation of statistical properties, the 'Model Builder' (ESRI 2012), a visual programming language embedded in ArcMap for geoprocessing in batch mode, was used to build and automate parts of the work flows. Subsequent analysis concerned the identification for each point in one geometry of its closest, in the second geometry by means of NNV analysis (Clark and Evans 1954). At this point, the calculation involved only one nearest neighbour point between the two geometries and not an average of several closest points, providing a more straightforward procedure and avoiding spatial smoothing of the initial PSI measurements. Since PS points are not regularly distributed in space, the distances from their closest neighbour in the opposite geometry is varying, in our case from was less than one metre to 332 m with an average at 52 m. Representation of these distances is shown graphically in Figure (5), enabling proper inspection of the effect of slope aspect. Indeed, larger distances are required to connect points located over slopes viewed only by single SAR geometry Figure (6). For gently sloping regions, closest point distances remain short due to availability of PS targets from both ascending and descending passes.
To avoid the connection of points separated by long distances since they might not be sharing common behaviour in terms of ground displacement, a maximum search radius was introduced as a constraint, limiting the analysis within distances shorter than 150 m. To this end, the distance constraint was arbitrary defined as being three times higher than the average distance between PS points. Points failed to meet the above requirement were rejected from further analysis. Similar constraint was applied when defining search window parameter for GRID approach. Although constraining distance does not affect actually the performed analysis, it ensures reduced extrapolation of values and thus, more reliable displacement maps.
Having defined the links between ascending and descending point targets, was allowed to join their attributes in the geodatabase domain, while the independent point data-sets were combined in a unique vector layer. The increase in density obtained by merging the PSI measurements in a vector domain is shown in Figure 3(c). The calculation of the parameters of interest from the PSI LOS measurements available for each point, namely the vertical and E-W motion components, was then achieved using simple mathematic operations between the corresponding fields in the built geodatabase.
Apart from implementing the proposed NNV approach, the conventional GRID approach was applied to permit the critical assessment of the obtained results. In short, GRID post-processing, contrary to NNV, involved the interpolation of each vector data independently, followed by the combination of the output interpolated surfaces by applying simple calculations between the raster layers. In that case, the final error was estimated by also combining the interpolation errors that describe the independent raster layers as proposed by Agudo et al. (2006).
Of importance is the fact that for both GRID and NNV approaches, processing parameters were equally defined (e.g. maximum interpolation distance in GRID vs. maximum search radius in NNV) to facilitate robust comparisons between them and to investigate potential effects on their outputs.

Resampling effects
To demonstrate the effect of resampling on point measurements, independently of the interpolation technique, the initial analysis involved the evaluation of the real-valued pixels (ones containing PS points), contrary to the 'interpolated' ones, for which no PS are available; hence, spatial interpolation would be then required to obtain a PSI estimate. At this stage, pixel values are not examined, but the availability or not of PS points within the resolution cell for each grid spacing.
Since the NNV approach is addressed to the combination of ascending and descending PS, only the pixels with available point targets from both geometries within the resolution cell were assigned to 'data' class. Otherwise, if only one component was present the pixel was allocated to the 'interpolated' class. However, to permit a proper comparison, the NNV-combined PS data-set was further constrained according to the analysis pixel size considered each time. Points with nearest neighbour distances between opposite geometries exceeding the largest dimension of the specific grid cell corners (√2x 2 ), for each analysis resolution, were identified, and not further considered in this part of the investigation. In principal, this is comparable to the maximum point search distance of spatial interpolation algorithms. Analogous classification into 'data' and 'interpolated' pixels was also followed for the GRID approach. In that case, the classification was performed during the resampling of the PSI results independently for each SAR geometry and prior to their combination into single raster product (see Figure 4).
The increase in the initial 'data' pixels by NNV can be easily observed in line with the considerable increase of inputs point vectors as indicated by the calculated density maps (Figure 3). PS densities were calculated based on hexagon-shaped grids, as proposed by Agudo et al. (2006) for the generation of PSI density maps.
The spatial distribution of data pixels for both NNV and GRID approaches, additional data pixels for NNV, as well as the ones for which interpolation is required following any combination strategy are shown in Figure 7. Investigation is performed on various resolutions for completeness. It is obvious Figure 6. Map representations of lines connecting ascending and descending PS point targets used to transfer attributes between geodatabases in NNV approach (magnified region in Figure 5) (up) and overlay of PS points on DEM heights showing the effect of terrain aspect vs. SAR acquisition geometry (layover and shadow) on the availability of PS points and NNV distances (down). For simplicity, when mutual exchange of attributes occurs between ascending and descending points, only one connecting line is shown (overlapping symbols).

Figure 7.
Distribution of 'data' and 'interpolated' pixels (see manuscript for definitions) for NNV and GRID approaches and additional 'data' pixels gained by NNV for the various spatial resolutions. The percentage of contribution for each pixel class is presented in Figure 8.
that with the increase of the pixel size (decrease resolution) the 'interpolated' pixels are reduced, while in all cases additional 'data' pixels are gained by the NNV approach. As expected, significant number of 'interpolated' pixels are still dominant in finer resolutions.
Examining in detail the trade between 'data' and 'interpolated' pixels for both NNV and GRID approaches in various resolutions (Figure 8), it could be stated that equal number of 'data' and 'interpolated' pixels is reached by decreasing the grid resolution rather earlier for NNV than for GRID approach. This observation directly implies that for the combination of different LOS geometries, more robust results f, in terms of the number of real-valued data pixel could be obtained in higher resolution by NNV (~60-80 m) than for GRID (~120 m). Characteristic remains the fact that at 40-50 m grid resolution for NNV there seems to be a balance between data and interpolated pixel (50% each), whereas for GRID real data pixels cover less than 20% of the area of interest (Figure 8).

Spatial interpolation effects
The positive effect of the proposed NNV approach in increasing the spatial coverage of actual data pixels has been demonstrated in a more theoretical way in the previous section. However, resampling of irregular distributed vector data by simple averaging of points' values at each resolution cell) is not applied in practice. Most often continuous spatial coverage is of interest and thus, interpolatiion over no data regions is applied.
Most commonly applied interpolation technique for PSI results is based on Kriging (Oliver & Webster 1990), as proposed during the validation of the PSI technique in various environments (Agudo et al. 2006). Although not a deterministic interpolator, among the advantages of kriging is the estimation of the variance of the predicted surface. Kriging variance is independent of the measurements, however, can be used as an estimation of the accuracy of the interpolated surface. In practice, it describes the probability that the actual value correspond to the predicted one plus or minus times the square root of the variance, assuming errors are normally distributed.
Following the recommendations of Agudo et al. (2006), ordinary kriging using an exponential semivariogram model was applied to the PSI vector data. Instead of considering a maximum number of points for the interpolator, resulting in a smoother surface, a predefined search radius as per previous analysis was considered. To ensure spatial continuity of the prediction surface in cases were no points were identified within the defined radius, an additional parameter for minimum number of points (set to eight points) was considered.
The interpolation errors obtained for the NNV and GRID approaches are shown in (Figure 9). The impact of the reduced number of PS points in terms of error budget when interpolation is performed on individual PS results (LOS geometries) is rather evident. Furthermore, standard error increases significantly for GRID when combining the two interpolated surfaces compared to the single NNV output. The above observation is quite reasonable since the accuracy of the predicted surface is reduced with increasing distances from the initial inputs points. In addition, the final error of the GRID is estimated by combining in a similar procedure the ascending and descending errors which results in higher overall error values (Figure 10(c)). It is important to underline that the described improvement gained by NNV refers exclusively to errors introduced by the spatial interpolation procedure which add up to the measurements' uncertainties.

Effects on displacement estimates
Since the analysis was performed based on actual processing results, displacement values calculated for each approach were also examined as a final step. No further interpretation of the results in terms of actual geophysical signal was considered for the purpose of this paper.
The correspondence between GRID and NNV is presented for the vertical motion component at a resolution of 60 m (Figure 10(a)). A relatively high correlation based on ordinary linear regression (R 2 = 0.96) is observed. In fact, NNV-GRID correlation subtly increases with the increase of the size of grid spacing. This is expected as the number of 'interpolated' pixels (see Section 4) are considerably reduced for grids with large pixel sizes. The consistency between the two approaches is basically anticipated, otherwise, significant deviations would indicate a non-proper utilization of the PSI results.
Actually, when examining the spatial distribution of the deviations between the two approaches ( Figure 11), it can be seen that differences exist and appear localized over mountainous regions. This is more visible for finer grid resolutions. The calculated deviations vary from approximately ± 10 mm yr −1 to ± 7 mm yr −1 for finer to coarser grid resolution, respectively. Indeed, the spread observed in the scatterplot of Figure 10a corresponds to these regions. In fact, only limited part of the region of interest exhibits rough terrain, explaining the small number of deviating pixels (low dispersion in Figure 10 and their corresponding spatial distribution in Figure 11). Regions with more dominant terrain would potentially present large extent and magnitude of deviations. Figure 9. Interpolation errors of PS vector points for the ascending and descending geometries (up) and the combined PSI results following NNV and GRID approaches (down) for 80 m grid spacing. Different colour scaling was selected for visualisation purposes.

Results
For the purposes of this study, point targets obtained from the PSI analysis of SAR data over Nisyros volcano were utilized to demonstrate the advantages of the NNV approach. The reduced number of PS points over gently sloping regions should be attributed to temporal decorrelation effects and the lack of point-like scatterers over natural terrain. On the contrary, missing PS over specific slope aspects (W-NW for ascending and E-SE for descending) is directly related to the slant acquisition geometry of the SAR. Modelling the relationship between slope aspect and SAR LOS geometry on SAR measurements and specifically SAR-derived DEMs is further discussed in Foumelis et al. (2015). Given the comparable PS densities for ascending and descending orbits, a thorough comparison was performed, while data gaps presented an opportunity to investigate the actual performance of the proposed approach.
The analysis included the separate examination of each aspect affecting common practices for combining PSI data. When interpolation is not considered, the error sources are dominated by the resampling options, while if continuous surfaces are targeted, then interpolation effects are contributing as well.
In practice, for the conventional GRID approach, by considering the ascending and descending PSI data-sets separately during resampling, there is a substantial reduction of the spatial coverage if interpolation is not applied (limited number of real data pixels). However, if interpolation is considered finally, an increase in uncertainties from multiple interpolations should be expected. For the former case, to overcome coverage issues large cell sizes could be selected at the expense of resolution, which in turns results in a significant smoothing of the initial PSI measurements. For the latter, unavoidably extensive interpolation over data gaps (layover and shadow areas) to obtain a spatially continuous results, introduce errors that cumulatively propagate to the final result increasing the overall measurement uncertainties.
In the NNV approach, the consideration of single nearest neighbour in NNV practically eliminates on the one hand the effect of averaging, on the other hand it allows interpolating PSI results on a finer resolution grid, an important aspect for many geohazard applications. The differences between the approaches in terms of error budget are two folded and clearly caused by the number of interpolations involved and the density of points utilized during the interpolation procedure. For GRID, interpolation is performed twice based on the independent PSI data-sets, while for NNV a single interpolation is considered using combined PS point stacks, increasing thus the density of the inputs.
For what concerns the actual displacement values, comparable results to the conventional approaches were obtained by NNV for flat and gently sloping regions, except from local deviations identified over high terrain slopes (approximately larger than 23-36°). The above fact implies that no significant improvement is expected over gently sloping regions or areas of low relief. Nonetheless, over high relief regions differences of up to 8-9 mm yr −1 are observed, partly augmented by the fact that no spatial smoothing is performed in NNV. As expected, differences over regions of actual data are minor, but become significant in areas without original inputs (pixels with interpolated values). In fact, a difficulty in evaluating the performance arises from the lack of external source of information (e.g. GPS or levelling data) for validation purposes. However, the high correlation between conventional GRID and NNV over relatively flat area signifies the correctness of the approach.
Finally, regarding the combined PSI LOS uncertainties (√((σ a 2 )+(σ d 2 )), where σ a and σ d the uncertainties for ascending and descending orbits, respectively), alike to the displacement values, they should be comparable for both NNV and GRID approaches. Indeed, high correlation is observed, yet having larger spread from the optimum theoretical correlation (Figure 10(b)), compared to the displacement estimates (Figure 10(a)). Deviations (up to 1 mm yr −1 ) are likewise limited over regions of high terrain where larger distances between ascending and descending points occur. It is interesting to note that no relationship between the sign of deviations and the applied approach, that could suggest a systematic over-or under-estimation of the displacements, was recognized.  Figure 11. Deviation in vertical displacement rates between NNV and GRID at 80 m grid spacing for all pixel classes (see Figure 7) (a), NNV and GRID data pixels (b), NNV data and GRID interpolated pixels (c) and NNV and GRID interpolated pixels (d).

Conclusions and discussion
The advantages of the proposed vector-based NNV approach compared to the commonly used gridbased procedure for combining PSI ascending and descending LOS observables is being demonstrated using real data as a case study. The NNV approach combines PSI data in the vector domain, performing necessary attributes transfer and calculations between feature geodatabases prior to any rasterisation. By increasing the number of input vector points during subsequent interpolation the error budget coming from spatial interpolation is being reduced. Elaborating the current findings, it appears that simple resampling is applicable only when sufficiently well-distributed PS targets are available. In that case, smaller pixel sizes could be selected with NNV, ensuring in the meanwhile better spatial coverage. The option to proceed without further interpolation of the combined PSI data-set is an alternative option for NNV (Figure 4). When filling of no data regions is of interest, to obtained a continues surfaces, NNV reduces the interpolationrelated errors simply by applying it only once. Even the use of the combined PS stack in NNV during interpolation reduces potential inaccuracies, especially over mountainous regions.
From the existing spatial interpolation techniques, Kriging method was chosen to better meet the application domain (Agudo et al. 2006). The triangular irregular network (TIN) methods were not given specific focus in this work. This is mainly due to their limited applicability to point measurements of high spatial variability and random density, such as the PSI data. Furthermore, despite being precise over input point locations, the lack of error estimate limits their choice as an interpolator. Triangulation could replace the regular-grid based interpolation steps in both the herein described GRID and NNV approaches. Independently though of its accuracy, since several smoothing options could be applied, when triangulation is considered within GRID post-processing, it suffers from analogous drawbacks of multiple interpolation steps. In fact, there are additional issues when considering triangulation as no combination is possible between surfaces having different triangular facets geometry. Thus, an additional rasterisation step is required to convert TIN surfaces into regular grids, which finally, adds to the overall errors and smoothing of the results.
The increase in output surface resolution together with the reduction of interpolation error budget for the combined PSI results in NNV is of significant importance during modelling of ground motion or integration of PSI measurements with GNSS observations for various geohazard applications. The proposed approach are applicable to any PSI data-set independently of the processing scheme. Yet, a main advantage that needs to be underlined, remains the potential to extract the Vertical and Easting motion components at a vector level without the need to interpolate or even resample. This is more evident when dense and well-distributed PS targets are available, which seems to be the case when using novel algorithms able to deduce displacement time series from targets exhibiting distributed scattering mechanism (Ferretti et al. 2011).
The functionalities of GIS towards such analysis are of valuable contribution, providing an easy way to manipulate vector data both in terms of features editing, as well as attribute operations within geodatabases. The analysis was performed taking into account a substantial range of output grid spacing (from 20 to 260 m), not applicable in practice, with the aim to concurrently investigate both resampling and interpolation effects. Still, the high demands of running numerous GIS operators multiple times for the various outputs resolutions, was addressed by utilizing the model building functionalities of the GIS for subsequent automated batch processing. The benefit of working with vector data, instead of raster layers, especially in saving storage space should be more pronounced when large stacks of PSI or wide area processing results are involved.
Future development should include the optimization of the NNV approach taking into account terrain aspect in relation to the SAR acquisition geometry for constraining the selection of nearest neighbours, avoiding 'linking' points across regions of SAR layover and shadow zones. The current implementation only considers geographic proximity of points from different SAR geometries disregarding variations in LOS motion patterns across ridges and valleys due to changes in terrain orientation. Such procedure would directly benefit cases where significant horizontal motion is expected, such as landslide studies. Further considerations may include a distance weighting procedure for multiple nearest neighbours selection controlling the smoothness of the derived surfaces. Apart from potential future algorithmic improvements, it was demonstrated in the current work that the use of available GIS functionalities can contribute towards more robust combination of ascending and descending PSI data. Finally, the implementation of the NNV post-processing in a more automated way within existing GIS packages should also be a major contribution, especially to non-expert GIS users, allowing a more straightforward exploitation of SAR interferometric results, where actual motion component is required.