## 2 Materials and Methods
### Accuracy Assessment of Bundle Block Adjustment, Orthomosaic and Digital Surface Model  
The accuracies of the BBA_traditional and BBA_PPK, namely of BBA (X, Y and Z), orthomosaic (X and Y) and DSM (Z) were assessed with the help of the VPs. The coordinates of VPs, as computed by BBA_traditional and BBA_PPK, were retrieved through Pix4 Dmapper, namely by using Checkpoints. They were marked in the same manner as GCPs, i.e., by measuring the center of the target (in this study, we marked each VP at 30 images). Then, the difference between the initial and computed position of the checkpoints (for X, Y and Z) was calculated and shown in the quality report [48]. The calculated differences were further used in comparative analyses. Coordinates of the VPs of the orthomosaics (X and Y) and DSMs (Z) were retrieved by using GIS software, namely ArcGIS Pro 2.7.3 [55]. The VP targets were vectorized based on the orthomosaic. To accurately determine the center of the target, the procedure of vectorizing was carried out in three steps: (i) In the first step, a polygon with the actual size of target was made, orientated in the same direction as the target, and centralized according to the center of the target, such that it was covered by a polygon. (ii) In the second step, two polylines were drawn from each polygon's angle, such that they crossed in the middle of the target. Lastly, (iii) in the third step, the point that represents the center of the target was placed at the crossing of the polylines. The location of the point was used to calculate the X and Y location from the orthomosaic. The whole procedure was performed three times, separately [37], and, in the end, an average value of three coordinates was used as the final one, which was used in further analyses. The Z coordinate was extracted from the DSM layer, based on the planimetric (XY) locations of VPs as measured in the field, namely by using the function Extract Values to Points. The coordinate differences, as mentioned above, were then used to compute several statistics. These concern the minimum, maximum, median, mean, standard deviation, and root mean square error (RMSE) of the differences between the coordinates measured on photogrammetric-derived products (X\({}_{\text{M}}\), Y\({}_{\text{M}}\) and Z\({}_{\text{M}}\)) and given (X\({}_{\text{VP}}\), Y\({}_{\text{VP}}\) and Z\({}_{\text{VP}}\)) coordinates of the VPs (Equations (1)-(3)) [63]:  
diff\[{}_{\text{X}}=\,X_{\text{VP}}-\,X_{\text{M}}\] (1)  
diff\[{}_{\text{Y}}=\,Y_{\text{VP}}-\,Y_{\text{M}}\] (2)  
diff\[{}_{\text{Z}}=\,Z_{\text{VP}}-\,Z_{\text{M}}\] (3)  
The RMSE, given by Equation (4) for the differences concerning each coordinate (C) separately (n is a total number of coordinates), was computed as an accuracy value in planimetry (Equation (5), [63]), and a total accuracy value (Equation (6), [63]) was also computed as an accuracy value in planimetry (Equations (4)-(6); [63]):  
\[\text{RMSE}_{\text{C}}=\sqrt{\sum_{1}^{\text{n}}\frac{\text{diff}_{\text{C}} ^{2}}{\text{n}}}\text{ with C}=X,\,Y,\,Z \tag{4}\]  
\[\text{RMSE}_{\text{xy}}=\sqrt{\text{RMSE}_{\text{x}}{}^{2}+\text{RMSE}_{ \text{y}}{}^{2}} \tag{5}\]  
\[\text{RMSE}_{\text{xyz}}=\sqrt{\text{RMSE}_{\text{x}}{}^{2}+\text{RMSE}_{ \text{y}}{}^{2}+\text{RMSE}_{\text{z}}{}^{2}} \tag{6}\]The computed differences (diff\({}_{\rm{X}}\), diff\({}_{\rm{Y}}\) and diff\({}_{\rm{Z}}\)) are shown in the form of histograms for each photogrammetric product and georeferencing method separately (Section 3.1, Section 3.2, Section 3.3). Additionally, we also plotted the differences spatially to study any correlation with the terrain topography; the differences were plotted for each photogrammetric product and georeferencing method (Section 3.4). To compare the accuracy of both the AT, using the BBA_traditional and BBA_PPK (X, Y and Z), and the produced orthomosaics (X and Y) and DSM (Z) for each of the three study areas, two statistical tests were carried out.