## 4 Methods
### DSM Extraction  
As has been already stated, we decided to use exactly the same subareas, reference LiDAR data and methodology as previously used by [PERSON] et al. [19] in order to allow the readers to compare the results obtained with different sensors, matching algorithms, software packages and land covers.  
In the current work, three stereo pairs (two from Deimos-2 and one from WV2) were used independently to extract the DSM in the study area. The software used was again OrthoEngine (PCI Geomatica v. 2018). OrthoEngine implements a hierarchical SGM approach based on the widely known algorithm proposed by [PERSON] [17] to generate the disparity maps after applying an epipolar rectification process to the original stereo images.  
After the sensor orientations for Deimos-2, following the best operational approach provided by the results described in Section 5.1, and after the WV2 stereo pairs were carried out, three grid spacing format DSMs for each subarea were extracted (two DSMs corresponding to the two Deimos-2 stereo pairs and one DSM corresponding to the WV2 stereo pair). All the DSMs were extracted in orthometric elevations using the EGM2008 geoid without filling blanks (no interpolation) and setting the grid spacing to 1 m [19; 21]. In short, 27 unfilled DSMs were extracted corresponding to the nine subareas and the three stereo pairs.  
According to several works [18; 19; 21; 20], the quality of the extracted DSMs was assessed by computing their completeness and vertical accuracy.  
The completeness of every DSM was computed for the nine studied subareas as the ratio between the number of correctly matched points for the stereo-photogrammetrically extracted DSM and the maximum possible number of points corresponding to the selected DSM grid spacing.  
The vertical accuracy assessment of the DSMs derived from the Deimos-2 and WV2 stereo pairs was carried out by using the 3D points from the benchmark LiDAR DSM (see Section 3.4) as ICPs, computing vertical residual (z-residual) at each corresponding point by subtracting the LiDAR height from the photogrammetrically derived height. Note that each ICP will produce a z-residual only if the photogrammetrically derived DSM presents height information in the area around the planimetric position of this ICP.  
The exploratory data analysis included standard accuracy statistical measures such as mean, standard deviation (SD), RMSE\({}_{\text{Z}}\) and 95 th percentile linear error (LE95). These statistics were computed to assess the final vertical accuracy after removing outliers from the z-residuals populations by applying the three-sigma rule [41]. In addition, and due to the likely presence of outliers in the measurements taken from the original dataset, an additional robust statistic such as normalized median absolute deviation (NMAD) [44] was computed over the residuals. NMAD is a measure of scale or variability that can be considered a consistent estimator for the estimation of standard deviation, also offering the advantage of being very insensitive to the presence of outliers. Note that the three-sigma rule was not applied in the case of the NMAD calculation.  
For studying the statistical influence of the two studied factors (i.e., sensor and land cover) on DSM quality, an experimental design based on a factorial model with three samples was used. The Kruskal-Wallis H test [45], a well-known rank-based non-parametric test, was applied to analyse statistically significant differences (\(p<0.05\)) between two or more groups of an independent factor in relation to a quantitative dependent variable (mean, SD, RMSE\({}_{\text{Z}}\), LE95, NMAD and completeness). It is important to note that, the residual populations (z-residuals at ICPs) did not always fit a normal distribution, and therefore the ANOVA test could not be used.