1. Introduction

Figure 1. Types of Digital Elevation Models (DEM) [7] (p. 47).

Figure 1. Types of Digital Elevation Models (DEM) [7] (p. 47).

•

DTM (digital terrain model): A DEM that records the boundary between the lithosphere and the atmosphere, without biosphere and anthroposphere, also called a “bare-earth” DEM. The treatment of hydrosphere, cryosphere, and voids (e.g., excluded buildings, water and trees) must be specified and clearly localized, e.g., by respective masks [3] (p. 7) (Figure 2).

•

DSM (digital surface model): A DEM that records the lower boundary of the atmosphere [5] (p. 7) (and either the lithosphere, hydrosphere, cryosphere, biosphere or anthroposphere) (Figure 2).

2. Materials and Methods

2.2. The case study: the plateau of Aptera

Aptera is an ancient city situated in Greece, near Souda Bay, Crete. It is located on a plateau that spans approximately 250 hectares.

Figure 3. Google maps – the Souda bay.

Figure 3. Google maps – the Souda bay.

Figure 4. Google maps – the Aptera plateau.

Figure 4. Google maps – the Aptera plateau.

Now managed by the Greek Ministry of Culture, Ephorate of Antiquities, over the centuries, the site has experienced various occupations and influences:

• Minoan civilization (2600 BC - 1100 BC): the oldest occupation of Aptera. Remains of houses, tombs, and everyday objects have been found.
• Second Minoan civilization (2600 BC - 1100 BC): the oldest occupation of Aptera. Remains of houses, tombs, and everyday objects have been found.
• Classical period (480 BC - 323 BC): Aptera became an independent Greek city-state and the theater was built. Roman period (323 BC - 395 AD): during this period, the baths, the thermal baths, and the Roman forum were built.
• Byzantine period (395 AD - 1204 AD): new walls and fortifications were built to protect the city from Arab invasion, and it is believed that an earthquake in the 7th century caused the disappearance of the ancient city.
• Venetian period (1204 AD - 1669 AD): during this period, Aptera was occupied by the Venetians, who built a castle and fortifications.
• Ottoman period (1669 AD - 1898 AD): during the Ottoman occupation, the castle was used as a prison and a defensive place.

Figure 5. Topographical plan of Aptera [20] (p. 17)

Figure 5. Topographical plan of Aptera [20] (p. 17)

2.3. Aerial photogrammetric survey

The surface area to be mapped is approximately 250 hectares. To ensure a comprehensive and accurate survey of the area, it was decided to proceed via aerial photogrammetry using a low-cost drone, with automated flight through planning and support - during the photogrammetric reconstruction phase - of GPS coordinates acquired through a GNSS antenna.

2.3.1. Unmanned aerial system (UAS)

For the aerial photogrammetry work via UAV, the DJI Mini 2 drone was used. Our Unmanned Aircraft System (UAS) included two DJI Mini 2 drones, two DJI-RC-N1 controllers, 9 batteries, and two power banks for the continuous recharging of depleted batteries. As we will see in the flight planning phase, equipping ourselves with multiple drones, controllers, and batteries proved necessary to face the vast extents of the area and the windy conditions (from 8 to 30 km/h) prevailing in the project area. Wind indeed represents a critical factor that impacts the battery life and the efficiency of the survey mission. It is also crucial to pay close attention to ensure that the wind or any gusts do not exceed the dynamic capabilities of the drone itself which, in the case of the Mini 2, guarantees a maximum wind resistance of up to 36 km/h.

2.3.2. Initial flight strategy

Given the extensive nature of the case study and the necessity to adhere to tight survey timelines for the workshop, an automated survey flight through planning was chosen rather than manual. Flight planning allowed for the optimization of survey times, the reduction of operator errors, and ensured the acquisition of predefined results, thanks to the a priori definition of parameters such as: this:

•

the percentage of vertical and horizontal overlap between successive shots,

•

flight altitude,

•

confinement of the flight area,

•

the pattern of the flight path,

•

UAV speed,

•

flight duration,

•

the number of batteries to be used,

•

the desired resolution of the orthophotogrammetric output.

However, it’s important to emphasize that office flight programming is not sufficient to ensure the mission’s success. Indeed, it’s necessary to inspect the area before initiating the flight to identify potential obstacles or dangers that could not be detected beforehand, such as poles, pylons or electrical wires, and potentially modify the flight plan accordingly, as described in Chapter 2.3.3.

For this case study, we opted for the Dronelink software solution, which offers a web app for PC flight mission planning and a smartphone app for mission planning, drone connection, and control during flight. Dronelink offers various features depending on the chosen plan, divided between professional plans (periodic payment with three solutions: starter, growth, and enterprise) and hobbystic plans (one-time payment with three solutions: basic, premium, and elite). For this case study, we adopted the premium plan, which proved ideal for planning areas and micro-areas with nadir flight at a constant height.

The total area was divided into 5 zones, each of which was further divided into smaller areas, varying in number from 3 to 6. The size of the individual areas varies between 8 and 12 hectares, although there are exceptions determined by the morphology of the area to be surveyed.

Figure 6. Initial flight strategy.

Figure 6. Initial flight strategy.

The flight altitude was determined based on the desired spatial resolution, expressed by the Ground Sampling Distance (GSD). The GSD is the distance, measured on the ground, between two adjacent pixels in the image. For the present study, the desired GSD is about 2 cm/pixel, therefore the flight altitude was identified at an altitude of approximately 40-50m.

However, the Premium plan does not allow the flight altitude to be set as a variable function relative to the ground but as a constant function relative to the “home” point, which makes the definition of flight altitude dependent on the altitude of the starting point and the resulting GDS variable according to the morphology of the surveyed area. Moreover, it should be noted that, as the workshop’s focus was on surveying the upper and predominantly flat area of the plateau, the altitudes and the GDS refer to the average elevation of this area, to the detriment of the quality of sloping areas or at the foot of the plateau that can descend up to 150 meters in elevation, but which do not constitute the primary objective of the work.

The vertical overlap between shots and the lateral overlap between strips was set at 75% to ensure a larger safety margin. Although a reduction in overlap would have resulted in time and battery savings, it would have compromised the safety of the area’s survey coverage, which was not considered acceptable given the time constraints required by the workshop.

The different areas overlap by about 40 meters, equivalent to about 2 stripes, to ensure continuity in the photo processing process and prevent any gaps in the image.

The UAV speed was set at the default value of 16.1 km/h to avoid flying too quickly and obtaining blurry photographs.

The average extent of each area, around 10 hectares, was deemed optimal considering the limitation of not exceeding a maximum distance of 1.5 km between the operator and the drone for safety reasons.

Below is a schema that summarizes the defined input parameters and the output results obtained from Dronelink’s “Mission Preview” function, in terms of completion time, batteries used, and resulting GSD with respect to the home point, based on the entered data.

Table 1. Initial flight strategy.

Table 1. Initial flight strategy.

Zone	Area	hectares	Altitude(m)	Flight duration	N° batterires	GSD (cm/pix)
1	1	12.1	40	00:36:40	9	1.42
	2	9.7	40	01:06:47		1.42
	3	9.0	40	01:33:20		1.42
	4	9.3	40	02:01:46		1.42
2	1	11.1	40	00:37:40	9	1.42
	2	12.3	40	01:15:20		1.42
	3	12.8	40	01:52:56		1.42
	4	4.3	40	02:09:54		1.42
3	1	10.4	40	00:33:16	7	1.42
	2	9.0	40	01:00:04		1.42
	3	10.1	40	01:32:25		1.42
4	1	9.9	40	00:32:07	8	1.42
	2	12.9	40	01:12:48		1.42
	3	11.0	40	01:47:04		1.42
5	1	7.0	40	00:21:35	7	1.42
	2	9.1	40	00:43:21		1.42
	3	11.4	40	01:12:38		1.42
	4	11.8	40	01:44:32		1.42

Please note that, in the absence of a precise indication of wind strength, the software’s estimated duration of battery life is 15 minutes.

2.3.3. Executed Flight plan

The flight plan prepared a priori proved to be fairly consistent with the situation found on site, despite the need for some modifications and the exclusion of some areas for reasons of force majeure. In particular, it was confirmed that the planned altitude settings could be maintained, as there were no poles or pylons in the vicinity.

However, the main problem was the impossibility of entering certain areas due to private property or prohibitive morphological conditions. Therefore, it was necessary to reduce the area covered by the survey and exclude some areas due to the too high distance that would be created between the operator and the UAV, considering possible signal losses and the obstructive action of the wind with respect to the control of the UAV. In general, distances greater than 1.5 km from the operator were avoided, as at these distances the signal was fluctuating and the energy consumption of the battery became too high compared to the promised effectiveness.

Below is a diagram representing the areas that were actually surveyed.

Figure 7. Executed flight strategy.

Figure 7. Executed flight strategy.

Below, the changes made to flight parameters and area extensions are being discussed.

Table 2. Executed flight strategy.

Table 2. Executed flight strategy.

Zone	Area	hectares	Altitude(m)	Flight duration	N° batterires	GSD (cm/pix)
1	1	12.1	40	00:36:40	9	1.42
	2	9.7	40	01:06:47		1.42
	3	9.0	40	01:33:20		1.42
	4	9.3	40	02:01:46		1.42
2	1	11.1	40	00:35:34	9	1.42
	2	12.3	40	01:14:03		1.42
	3	12.8	40	01:54:42		1.42
	4	4.3	40	02:09:06		1.42
3	1	10.9	50	00:27:54	5	1.78
	2	7.7	50	00:47:32		1.78
	3	9.0	50	01:09:55		1.78
4	1	9.6	50	00:24:26	4	1.78
	2	8.3	55	00:53:23		1.96
	3	-	-	-	-	-
5	1	6.4	40	00:22:05	4	1.42
	2	8.3	40	00:47:20		1.42
	3	-	-	-	-	-
	4	-	-	-	-	-

The extent of the areas and zones defined in the first approach phase was consistent with the objectives of field survey efficiency and digital reconstruction effectiveness in the office, except for some aspects of overlapping areas that could be optimized. The overlap between adjacent areas set at 40m in the first flight hypothesis was indeed slightly redundant in the post-first campaign processing phase: an overlap of 20 meters, equivalent to about the area swept by a stripe, would have been sufficient, even in terms of safety coverage. For this reason, and for the usual reasons of time and battery savings, the adjacent areas belonging to zones 4 and 5 were surveyed by overlapping a single stripe between their paths (e.g. areas 3.1-3.2.3.3 and 4.1-4.2).

To avoid being too far from the drone, the UAV’s starting (and return) point was moved with each battery change, trying to center as much as possible in the area to be surveyed. It was noted that the maximum distance of 1.5 km represented the limit beyond which the signal could regularly be lost, albeit momentarily.

After the first processing of the images, the quality achieved was found to be satisfactory, so it was decided to confirm the flight altitudes and possibly raise them for the benefit of the acquisition times (e.g. zones 3), provided that any additional altitude variations were made in the case of significant average altitude variations within the area compared to the starting/return point (e.g. zone 4).

The vertical overlap between shots and lateral between strips was deemed satisfactory based on the first digital reconstructions and for this reason it was maintained as previously indicated.

The UAV speed of 16.1 km/h proved to be adequate for the purposes of the survey, although less favorable when returning home (RTH) due to the large distance from the home point at low speed, which involved a significant battery consumption. Despite this, it was decided to maintain this speed and, if necessary, to manually manage the return home without incurring, in this way, the speed limits of the automatic flight.

Finally, we end the report on the actual estimates of time and battery consumption: in general, the estimates were respected, with an autonomy of about 15-20 minutes per battery, except in the acquisition phases where the wind was blowing stronger, reducing battery life and consequently increasing the actual survey time.

2.5. Photogrammetric data processing

At the end of each day of field data collection, once we returned to the office, the focus shifted to the data processing phase. This essential step allowed to verify the quality of the information collected in the field and to provide the working teams involved in the workshop with fundamental material, such as point clouds, meshes, DEMs and orthophotos, in a very short time, for their daily commitment.

For the photogrammetric reconstruction activity, Agisoft Metashape software was used.

In this chapter, it will be presented the workflow used to generate the typical outputs for a digital reconstruction project of a territorial-scale object.

In the next chapter, instead, we will focus on the analysis of the quality and accuracy of the results obtained: as already mentioned, for this second purpose we will examine a smaller area within the same plateau, for which the same workflow that we will expose in the next paragraphs was followed.

Below, it is presented a brief summary of the phases of the reconstruction process that we will soon discuss in more detail:

•

Photo alignment

•

GCP and QCP definition

•

Alignment optimization

•

Dense cloud

•

Point classification

•

2D-DEM (DTM)

•

Contour lines

•

DSM-TIN and DTM-TIN (Mesh with texture)

•

Orthomosaic

2.5.1. Photo alignment

The first step in our process involved importing the 8532 photos acquired, organized into folders corresponding to individual missions. Subsequently, we estimated the camera positions using the “Align Photos” command.

The parameters and values used are as follows:

•

Accuracy: High

•

Generic preselection: yes

•

Referece preselection: yes, Source

•

Reset curent alignment: -

•

Key point limit: 40.000

•

Tie point limit: 4.000

•

Apply mask to: none

•

Exclude stationary tie points: no

•

Guided image matching: no

•

Adaptive camera model fitting: yes

Higher accuracy settings help to obtain more accurate camera position estimates. Lower accuracy settings can be used to get the rough camera positions in a shorter period of time. In the Generic preselection mode, the overlapping pairs of photos are selected by matching photos using lower accuracy setting first, benefiting the processing time. In the Source preselection mode, the overlapping pairs of photos are selected based on the measured camera locations (in this case, the GPS position present in the metadata of the drone photos). Key point and Tie point limit were indicated as recommended by the Metashape guide. Adaptive camera model fitting helps to prevent the divergence of some parameters for datasets with weak camera geometry, like a typical aerial dataset.

At this stage, it is possible to explore the Tie Points and observe, as indicated in blue in the following image, the positioning of the cameras.

Figure 8. Tie points and camera locations.

Figure 8. Tie points and camera locations.

Subsequently, the 10 detected points from the GNSS antenna were identified by placing their respective markers in all the photographs where those points were present.

Figure 9. Tie points and markers identification.

Figure 9. Tie points and markers identification.

The next step involved georeferencing the model by importing the point coordinates into Metashape.

Firstly, the reference system was modified from the default WGS84 (EPSG::4326) to WGS84 (EPSG::32635) which corresponded to the system of the point coordinates obtained through GNSS. Subsequently, the camera information was automatically converted by selecting the appropriate option.

Table 1. GNSS survey – points coordinates.

Table 1. GNSS survey – points coordinates.

Marker	Easting (m) - survey	Northing (m) - survey	Altitude (m) - survey
1	512663,286	3924086,250	217,233
2	512680,138	3924117,120	213,147
3	512665,874	3924204,077	205,814
4	512699,736	3924220,669	203,286
5	512747,623	3924166,389	206,227
6	512748,647	3923848,650	211,528
7	512506,672	3923891,421	218,460
8	512390,339	3924087,504	219,545
9	512654,232	3924270,943	208,973
10	513234,929	3924361,905	192,365

Subsequently, the corresponding .csv file was loaded using the “Import Reference” command, and the Markers were automatically updated with the coordinates from the file.

At this point, we moved on to the alignment optimization phase.

The GPS information associated with the drone-derived cameras was deactivated as it was not desired to use them as a reference for alignment optimization due to their lower reliability compared to the newly inserted GPS coordinates.

Markers to be used for optimization were then selected, distinguishing between GCPs (Ground Control Points) and QCPs (Quality Control Points). For the entire plateau model, due to the lack of GCPs, all Markers were used as GCPs. However, for the detailed model used in the accuracy analysis chapter, as will be seen, various combinations of GCPs and QCPs were tested to study different results.

In the reference settings panel, the assumed accuracy of the measured values was specified (0.005m). The “Optimize Camera Alignment” command was used to further refine parameters related to the photographic data.

Once completed, the new camera alignment was defined, and new tie points were created.

It is now possible to proceed to the next phase: generating the dense point cloud.

2.5.2. Dense point cloud generation and management

In order to achieve a good result in a timely manner, the dense point cloud was generated with the following settings:

Dense cloud generation:

•

Quality: Medium

•

Depth filtering: Mild

•

Calculate point colors: yes

•

Calculate point confidence: yes

This resulted in a point cloud consisting of 256,572,086 points.

Figure 10. Dense point cloud generated.

Figure 10. Dense point cloud generated.

Subsequently, the “Filter by Confidence” command was used to assess the point cloud created in terms of confidence. This function allows filtering points within a confidence range defined by the user, which can range from 0 to 255. It was decided to remove points with values in the range from 0 to 2.

Figure 11. Dense point cloud – confidence level classification.

Figure 11. Dense point cloud – confidence level classification.

Figure 12. Dense point cloud – 0-2 range confidence.

Figure 12. Dense point cloud – 0-2 range confidence.

As a result, the number of points in the point cloud decreased from 256,571,864 to 213,416,628, which represents approximately 80% of the original number.

To obtain a Digital Terrain Model (DTM), it is now necessary to classify the point cloud and remove points that do not represent the terrain, such as vegetation, man-made objects, water bodies, etc. The Metashape software provides options for manual or semi-automatic point classification, allowing for iterative refinement of the process by manipulating key parameters.

• Classify ground points.
  A semi-automatic method that allows the recognition of only ground points through the definition of values for 4 parameters is as follows:

  a.

  Max angle (deg): sets limitation for an angle between terrain model and the line to connect the point in question with a point from a ground class.

  b.

  Max distance (m): sets limitation for a distance between the point in question and terrain model.

  c.

  Cell size (m): Determines the size of the cells for point cloud to be divided into as a preparatory step in ground points classification procedure.

  d.

  Erosion radius (m): determines the indentation (in meters) from unclassified points, ground points within that area will be unclassified.

• Classify points.
  Another semi-automatic method that allows for classifying points into different categories (including an uncategorized category) by defining a confidence value ranging from 0 to 1.
  The categories include:

  a.

  Ground

  b.

  High vegetation

  c.

  Building

  d.

  Road surface

  e.

  Car

  f.

  Man-made object

• Assign class:
  A manual method that allows for selecting points of interest and modifying their class. This can be a valuable support to semi-automatic classification to quickly correct the classification of specific areas that were not properly classified by the software.

The classification workflow consisted of multiple phases:

• Ground point classification.
  Through an iterative method, values were determined that allowed for an initial satisfactory classification of the point cloud.

  a.

  Max angle (deg): 15

  b.

  Max distance (m): 0.05

  c.

  Cell size (m): 20

  d.

  Erosion radius (m): 0

  The point cloud was then classified into the following categories:

  a.

  Unclassified points (grey)

  b.

  Ground points (brown)

  c.

  Low noise points (pink)

2.

Assign points.

Some points that clearly belonged to the “ground points” category were not included. Therefore, these points were selected and their classification was manually set. Among these points, for example, there were some rocks that the software may have mistaken for trees, bushes, man-made objects, or background noise.

Figure 14. Manual points classification.

Figure 14. Manual points classification.

3.

Classify points.

Before proceeding to the second phase of classification, it was necessary to reset the existing classification to avoid the “wrong point count” error. This operation was performed after removing points classified and confirmed as non-ground, preserving and refining the previous work. As a result, the number of points in the point cloud decreased from 213,416,628 to 101,715,919.

Subsequently, the “classify points” operation was executed. After a series of iterative tests, it was determined that a confidence value of 0.7 was appropriate. The classification was limited to the “ground” and “high vegetation points” classes, with the goal of refining the classification of the terrain point cloud by removing further vegetation components. Points belonging to the “high vegetation” category were then eliminated, resulting in a reduction of the point cloud from 101,715,919 to 93,193,523 points.

Figure 15. Point cloud classification.

Figure 15. Point cloud classification.

4.

Assign points.

At this stage, all the remaining points were manually reclassified as “ground points.” Finally, some less significant outlier points were removed from the point cloud. As a result, the number of points in the point cloud decreased from 93,193,523 to 87,865,204. Therefore, compared to the initial point cloud, the final point cloud contains 34% of the starting points.

Figure 16. Definitive ground points classification.

Figure 16. Definitive ground points classification.

2.5.3. DEM-2D (DTM) generation

At this point, we proceeded with the creation of the two-dimensional Digital Terrain Model (DTM). To achieve this, we selected the newly created terrain point cloud model and executed the “Build DEM” command.

The following are the parameters and values used:

•

Projection type: Geographic, GGRS87

•

Source data: Dense cloud

•

Interpolation: Enabled (default)

•

Point classes: Ground

•

Resolution: 0.111829 (not editable)

•

Total size (pix): 27936x21837

Figure 17. 2D-DTM model.

Figure 17. 2D-DTM model.

2.5.4. Contour lines generation

Using the produced DEM, contour lines of the model can now be generated with the “Generate contours” command.

The following parameters were set for this purpose:

•

Source data: DEM

•

Min altitude (m): -4.52

•

Max altitude (m): 231.70

•

Interval (m): 10 or 1

•

Prevent intersections: yes

Figure 18. 2D-DTM model with contour lines (10m interval).

Figure 18. 2D-DTM model with contour lines (10m interval).

2.5.5. DEM-TIN (DSM and DTM) generation

The DEM-TIN models for both the DSM and DTM were created using the “Build mesh” command to construct the geometry and the “Build texture” command to apply color data.

For the DEM-TIN model of the DSM, the following parameters and values were set:

•

Source data: Dense cloud (the optimized point cloud model based on confidence)

•

Surface type: Arbitrary (3D)

•

Quality: -

•

Face count: Medium (23.000.000)

•

Interpolation: Enabled (default)

•

Depth filtering: -

•

Point classes: All

•

Calculate vertex colors: No

•

Use strict volumetric masks: -

•

Reuse depth maps: -

A mesh consisting of 26,452,815 faces was then generated. For the DEM-TIN model of the DSM, the following parameters and values were set:

•

Source data: Dense cloud (the final DTM point cloud model)

•

Surface type: Height field (2.5D)

•

Quality: -

•

Face count: Medium (5.857.860)

•

Interpolation: Enabled (default)

•

Depth filtering: -

•

Point classes: All

•

Calculate vertex colors: No

•

Use strict volumetric masks: -

•

Reuse depth maps: -

A mesh with 5,769,130 faces was generated. The mesh was then cleaned in the outer areas, resulting in 5,700,490 faces.

Figure 19. Comparison: DSM model detail.

Figure 19. Comparison: DSM model detail.

Figure 20. Comparison: DTM model detail.

Figure 20. Comparison: DTM model detail.

Lastly, the texture was generated, using identical parameters for both the DTM and DSM models:

•

Texture type: Diffuse map or Occlusion map

•

Source data: Images

•

Mapping mode: Generic

•

Blending mode: Mosaic (default)

•

Texture size/count: 8000x8 o 16000x1

•

Enable hole filling: yes

•

Enable ghosting filter: no

Figure 21. DTM (TIN) model with diffuse map texture.

Figure 21. DTM (TIN) model with diffuse map texture.

Figure 22. DSM (TIN) model with occlusion map texture.

Figure 22. DSM (TIN) model with occlusion map texture.

2.5.6. Orthomosaic generation

Finally, an orthomosaic of the plateau was created using the “Build orthomosaic” command.

The following settings and parameters were applied:

•

Projection Type: Geographic – GRS87 (EPSG::4121)

•

Surface: Mesh DTM (2.5D)

•

Blending mode: Mosaic

•

Refine seamlines: no

•

Enable hole filling: no

•

Enable ghosting filter: no

•

Enable back-face culling: no

•

Pixel size (m): 3 cm/pix

At the “surface” option, you can choose between Mesh or DEM. However, the results demonstrate that the 2.5D Mesh provides better image projection onto the surface compared to DEM.

Figure 23. Orthomosaic.

Figure 23. Orthomosaic.

Figure 24. Orthomosaic detail 1.

Figure 24. Orthomosaic detail 1.

Figure 25. Orthomosaic detail 2.

Figure 25. Orthomosaic detail 2.

3. Results

3.1. Subsection

The work carried out has led to the generation of various types of output that can be used by workshop participants in various ways.

In particular, through exporting from Metashape, the following outputs were delivered:

•

Dense point clouds: exported in *e57 open format. These were primarily used as a record of the surveying work and for creating evocative images of the archaeological site. Views were captured through the appropriate command within Autodesk Recap, a software for managing the point cloud that is intuitive and useful for studying the cloud using bounding box, measurement, point visibility and color management, screenshot, and annotation functions. In order to not have to proceed more than once with the lengthy import of the *e57 file into Recap, a *rcs point cloud model was also created (Autodesk proprietary).

•

Textured DTM/DSM meshes: exported in *obj format for geometry, *jpg for texture image, *mtl for geometry-texture correlation. The mesh models represent the most manageable and exploitable output for three-dimensional project analysis and details. The DTM can be used in particular to extract environmental morphological sections, while the DSM carries three-dimensional information concerning also the vegetation and buildings present in the archaeological site. Two types of textures were generated: diffuse map and occlusion map. The “diffuse map” type presents color data from the mosaicking of image projections while the “occlusion map” type presents a grayscale color data from ambient occlusion. The diffuse map gives a suggestive realistic effect to the model, while the occlusion map highlights its geometries. To accommodate various solutions, these types of textures were created in two forms: a higher definition one (8000x8000pix with 8 textures) and a lower one (16000x16000pix with 1 texture). The obj format is manageable through different software, within the workshop the most used were Rhinoceros3D and Sketchup.

•

Contour lines: exported in *dxf format. The contour lines, which are vector in nature, were used to represent the terrain’s slope in plan in the two-dimensional CAD drawings.

•

Orthomosaic: exported in *tiff format. The orthomosaic represents the main two-dimensional raster output of the work. During export from metashape, it is possible to define the pixel size of the image (cm/pix, 3 in this case, as per orthomosaic generation), the total dimensions (36947x29437 pix, in this case, data derived based on the size of the model and the pixel size defined previously), the possibility of including the alpha channel, the background color, the creation of a Word file with associated metadata, the level of compression, etc... A great possibility of the *tiff is that it can be imported onto CAD at a 1:1 scale. The orthomosaic is particularly useful as a background for CAD drawings, to trace specific elements in nadir view or as a project masterplan or various detail documents.

Figure 26. Evocative view of the 3D point cloud.

Figure 26. Evocative view of the 3D point cloud.

3.2. Accuracy analysis

In this paragraph, the focus will be on the accuracy analysis of the digital terrain model.

As mentioned in section 2.4.2, the topographic survey was conducted with a limited number of GPS points, mainly concentrated on the main archaeological site for a specific project. Therefore, it has been decided to analyze the digital model within the larger contextual area, which measures 525x490m. To achieve this, the 1898 cameras related to the area were selected, and a new Metashape project was created using these cameras.

Figure 27. Tie points and camera locations.

Figure 27. Tie points and camera locations.

As observed previously, the cameras were aligned, tie points were generated, and the GPS points were accurately identified through their respective markers (9 out of the total 10 detected). Their accurate identification was achieved in all the images where they are present.

Figure 28. Tie points and markers identification.

Figure 28. Tie points and markers identification.

At this stage, a different approach was taken compared to what was done in chapter 2.5, where all the GPS points detected were used as GCPs. The tie points were duplicated eight times, in order to create multiple models to analyze different combinations of GCPs and Quality Control Points (QCPs).

Figure 29. GCP (red) and QCP (yellow) combinations.

Figure 29. GCP (red) and QCP (yellow) combinations.

For each of these models, the alignment of the cameras was optimized based on the respective GCPs. The accuracy of the corresponding QCPs was then analyzed using method 2 outlined in chapter 2.6.

Method 1, on the other hand, was only used in cases where all nine markers were used as GCPs. As mentioned before, it is more appropriate to discuss the precision of the model in this scenario rather than accuracy. The interpretation of the analysis result differs from that of method 2 cases, and this will be explained in more detail in chapter 4.

Below is a table displaying the initial data and the accuracy values obtained for each combination. The initial data was exported using the appropriate “Export reference” command.

Table 1. Accuracy assessment strategy: survey and digital model coordinates.

Table 1. Accuracy assessment strategy: survey and digital model coordinates.

Accuracy assessment strategy				Survey coordinates			Digital model coordinates
Combination	Marker	GCP	QCP	Easting xvi (m)	Northing yvi (m)	Altitude zvi (m)	Easting xci (m)	Northing yci (m)	Altitude zci (m)
1	1	1	1	512663,286	3924086,25	217,233	512663,24	3924086,259	217,221
	2	1	1	512680,138	3924117,12	213,147	512680,124	3924117,127	213,172
	3	1	1	512665,874	3924204,077	205,814	512665,927	3924204,142	205,78
	4	1	1	512699,736	3924220,669	203,286	512699,707	3924220,696	203,321
	5	1	1	512747,623	3924166,389	206,227	512747,457	3924166,393	206,263
	6	1	1	512748,647	3923848,65	211,528	512748,555	3923848,587	211,52
	7	1	1	512506,672	3923891,421	218,46	512506,748	3923891,274	218,454
	8	1	1	512390,339	3924087,504	219,545	512390,463	3924087,524	219,558
	9	1	1	512654,232	3924270,943	208,973	512654,327	3924271,02	208,923
2	1	1	0	512663,286	3924086,25	217,233	512663,261	3924086,281	217,234
	2	0	1	512680,138	3924117,12	213,147	512680,14	3924117,159	213,245
	3	0	1	512665,874	3924204,077	205,814	512665,944	3924204,169	205,747
	4	0	1	512699,736	3924220,669	203,286	512699,708	3924220,713	203,386
	5	1	0	512747,623	3924166,389	206,227	512747,504	3924166,427	206,237
	6	1	0	512748,647	3923848,65	211,528	512748,578	3923848,607	211,529
	7	1	0	512506,672	3923891,421	218,46	512506,72	3923891,312	218,459
	8	1	0	512390,339	3924087,504	219,545	512390,425	3924087,517	219,545
	9	1	0	512654,232	3924270,943	208,973	512654,311	3924271,014	208,963
3	1	1	0	512663,286	3924086,25	217,233	512663,263	3924086,252	217,256
	2	0	1	512680,138	3924117,12	213,147	512680,172	3924117,118	213,443
	3	0	1	512665,874	3924204,077	205,814	512665,95	3924204,108	205,94
	4	0	1	512699,736	3924220,669	203,286	512699,772	3924220,662	203,784
	5	0	1	512747,623	3924166,389	206,227	512747,636	3924166,394	207,043
	6	1	0	512748,647	3923848,65	211,528	512748,636	3923848,65	211,529
	7	1	0	512506,672	3923891,421	218,46	512506,682	3923891,399	218,444
	8	1	0	512390,339	3924087,504	219,545	512390,357	3924087,509	219,555
	9	1	0	512654,232	3924270,943	208,973	512654,239	3924270,958	208,956
4	1	1	0	512663,286	3924086,25	217,233	512663,257	3924086,348	217,256
	2	0	1	512680,138	3924117,120	213,147	512680,17	3924117,221	213,151
	3	0	1	512665,874	3924204,077	205,814	512666,066	3924204,242	204,491
	4	0	1	512699,736	3924220,669	203,286	512699,873	3924220,75	202,099
	5	1	0	512747,623	3924166,389	206,227	512747,597	3924166,419	206,212
	6	1	0	512748,647	3923848,65	211,528	512748,568	3923848,597	211,531
	7	1	0	512506,672	3923891,421	218,46	512506,735	3923891,301	218,451
	8	1	0	512390,339	3924087,504	219,545	512390,412	3924087,549	219,543
	9	0	1	512654,232	3924270,943	208,973	512654,46	3924271,256	206,385
5	1	0	1	512663,286	3924086,25	217,233	512663,247	3924086,301	217,275
	2	0	1	512680,138	3924117,12	213,147	512680,134	3924117,173	213,274
	3	0	1	512665,874	3924204,077	205,814	512665,936	3924204,177	205,753
	4	0	1	512699,736	3924220,669	203,286	512699,699	3924220,72	203,384
	5	1	0	512747,623	3924166,389	206,227	512747,493	3924166,435	206,239
	6	1	0	512748,647	3923848,65	211,528	512748,57	3923848,614	211,529
	7	1	0	512506,672	3923891,421	218,46	512506,716	3923891,317	218,459
	8	1	0	512390,339	3924087,504	219,545	512390,422	3924087,521	219,545
	9	1	0	512654,232	3924270,943	208,973	512654,307	3924271,021	208,962
6	1	1	0	512663,286	3924086,25	217,233	512662,866	3924085,648	216,789
	2	0	1	512680,138	3924117,12	213,147	512679,79	3924116,444	212,857
	3	0	1	512665,874	3924204,077	205,814	512665,654	3924203,407	205,241
	4	0	1	512699,736	3924220,669	203,286	512699,472	3924219,934	203,146
	5	0	1	512747,623	3924166,389	206,227	512747,259	3924165,658	206,611
	6	1	0	512748,647	3923848,65	211,528	512748,045	3923848,072	211,087
	7	0	1	512506,672	3923891,421	218,46	512506,263	3923890,813	216,808
	8	1	0	512390,339	3924087,504	219,545	512390,097	3924087,059	219,024
	9	1	0	512654,232	3924270,943	208,973	512653,992	3924270,221	208,1
7	1	0	1	512663,286	3924086,25	217,233	512663,269	3924086,313	220,17
	2	0	1	512680,138	3924117,12	213,147	512680,188	3924117,164	216,114
	3	0	1	512665,874	3924204,077	205,814	512665,96	3924204,055	207,371
	4	0	1	512699,736	3924220,669	203,286	512699,762	3924220,566	204,909
	5	0	1	512747,623	3924166,389	206,227	512747,624	3924166,372	208,963
	6	1	0	512748,647	3923848,65	211,528	512748,633	3923848,662	211,534
	7	1	0	512506,672	3923891,421	218,46	512506,675	3923891,39	218,45
	8	1	0	512390,339	3924087,504	219,545	512390,35	3924087,515	219,552
	9	1	0	512654,232	3924270,943	208,973	512654,232	3924270,951	208,97
8	1	0	1	512663,286	3924086,25	217,233	512663,224	3924086,25	219,973
	2	0	1	512680,138	3924117,12	213,147	512680,147	3924117,115	215,98
	3	0	1	512665,874	3924204,077	205,814	512665,949	3924204,033	207,303
	4	0	1	512699,736	3924220,669	203,286	512699,752	3924220,554	204,926
	5	0	1	512747,623	3924166,389	206,227	512747,586	3924166,346	209,043
	6	1	0	512748,647	3923848,65	211,528	512748,645	3923848,65	211,528
	7	0	1	512506,672	3923891,421	218,46	512506,666	3923891,124	217,576
	8	1	0	512390,339	3924087,504	219,545	512390,34	3924087,501	219,545
	9	1	0	512654,232	3924270,943	208,973	512654,233	3924270,946	208,973

Table 1. Accuracy assessment results.

Table 1. Accuracy assessment results.

Combination	n	RMSEx (m)	RMSEy (m)	RMSEr (m)	RMSEz (m)	Accuracyr (m)	Accuracyz [NVA] (m)
1	9	0,278	0,216	0,352	0,156	0,627	0,307
2	3	0,186	0,242	0,305	0,297	0,543	0,582
3	4	0,199	0,107	0,226	0,659	0,403	1,291
4	4	0,384	0,406	0,559	1,129	0,995	2,214
5	4	0,188	0,253	0,310	0,286	0,561	0,561
6	5	0,567	0,827	1,002	0,780	1,785	1,528
7	5	0,189	0,223	0,293	1,538	0,521	3,014
8	6	0,185	0,290	0,344	1,438	0,612	2,818

4. Discussion

References

1. Mora, O.E.; Suleiman, A.; Chen, J.; Pluta, D.; Okubo, M.H.; Josenhans, R. Comparing SUAS Photogrammetrically-Derived Point Clouds with GNSS Measurements and Terrestrial Laser Scanning for Topographic Mapping. Drones 2019, 3, 64. [Google Scholar] [CrossRef]
2. Jiménez-Jiménez, S.I.; Ojeda-Bustamante, W.; Marcial-Pablo, M. de J.; Enciso, J. Digital Terrain Models Generated with Low-Cost UAV Photogrammetry: Methodology and Accuracy. ISPRS International Journal of Geo-Information 2021, 10, 285. [Google Scholar] [CrossRef]
3. Remondino, F.; Barazzetti, L.; Nex, F.; Scaioni, M.; Sarazzi, D. UAV PHOTOGRAMMETRY FOR MAPPING AND 3D MODELING – CURRENT STATUS AND FUTURE PERSPECTIVES. The International Archives of the Photogrammetry, Remote Sensing and Spatial Information Sciences 2012, XXXVIII-1-C22, 25–31. [CrossRef]
4. Westoby, M.J.; Brasington, J.; Glasser, N.F.; Hambrey, M.J.; Reynolds, J.M. ‘Structure-from-Motion’ Photogrammetry: A Low-Cost, Effective Tool for Geoscience Applications. Geomorphology 2012, 179, 300–314. [Google Scholar] [CrossRef]
5. Guth, P.L.; Van Niekerk, A.; Grohmann, C.H.; Muller, J.-P.; Hawker, L.; Florinsky, I.V.; Gesch, D.; Reuter, H.I.; Herrera-Cruz, V.; Riazanoff, S.; et al. Digital Elevation Models: Terminology and Definitions. Remote Sensing 2021, 13, 3581. [Google Scholar] [CrossRef]
6. Pratibha P. Shingare; Sumit S. Kale; Review on Digital Elevation Model. International Journal of Modern Engineering Research (IJMER) 2013, 3, 2412–2418.
7. Digital Terrain Mapping and Analysis - Ppt Download. Available online: https://slideplayer.com/slide/5153395/ (accessed on 8 June 2023).
8. Polat, N.; Uysal, M. DTM GENERATION WITH UAV BASED PHOTOGRAMMETRIC POINT CLOUD. The International Archives of the Photogrammetry, Remote Sensing and Spatial Information Sciences 2017, XLII-4-W6, 77–79. [CrossRef]
9. Abdullah, A.F.; Vojinovic, Z.; Price, R.K.; Aziz, N.A.A. Improved Methodology for Processing Raw LiDAR Data to Support Urban Flood Modelling – Accounting for Elevated Roads and Bridges. Journal of Hydroinformatics 2011, 14, 253–269. [Google Scholar] [CrossRef]
10. Barazzetti, L.; Brovelli, M.A.; Valentini, L. LiDAR Digital Building Models for True Orthophoto Generation. Appl Geomat 2010, 2, 187–196. [Google Scholar] [CrossRef]
11. Trepekli, K.; Balstrøm, T.; Friborg, T.; Fog, B.; Allotey, A.N.; Kofie, R.Y.; Møller-Jensen, L. UAV-Borne, LiDAR-Based Elevation Modelling: A Method for Improving Local-Scale Urban Flood Risk Assessment. Nat Hazards 2022, 113, 423–451. [Google Scholar] [CrossRef]
12. Mazzoleni, M.; Paron, P.; Reali, A.; Juizo, D.; Manane, J.; Brandimarte, L. Testing UAV-Derived Topography for Hydraulic Modelling in a Tropical Environment. Nat Hazards 2020, 103, 139–163. [Google Scholar] [CrossRef]
13. Gil-Docampo, M.L.; Arza-García, M.; Ortiz-Sanz, J.; Martínez-Rodríguez, S.; Marcos-Robles, J.L.; Sánchez-Sastre, L.F. Above-Ground Biomass Estimation of Arable Crops Using UAV-Based SfM Photogrammetry. Geocarto International 2020, 35, 687–699. [Google Scholar] [CrossRef]
14. Dzierżek, J.; Kaczorowski, J.; Szymanek, M. Advantages of High-Resolution Digital Terrain Model (DTM) Analysis in Geological Cartography 2016.
15. Hugenholtz, C.H.; Whitehead, K.; Brown, O.W.; Barchyn, T.E.; Moorman, B.J.; LeClair, A.; Riddell, K.; Hamilton, T. Geomorphological Mapping with a Small Unmanned Aircraft System (SUAS): Feature Detection and Accuracy Assessment of a Photogrammetrically-Derived Digital Terrain Model. Geomorphology 2013, 194, 16–24. [Google Scholar] [CrossRef]
16. Dubbini, M.; Curzio, L.I.; Campedelli, A. Digital Elevation Models from Unmanned Aerial Vehicle Surveys for Archaeological Interpretation of Terrain Anomalies: Case Study of the Roman Castrum of Burnum (Croatia). Journal of Archaeological Science: Reports 2016, 8, 121–134. [Google Scholar] [CrossRef]
17. Rusnák, M.; Sládek, J.; Buša, J.; Greif, V. Suitability of Digital Elevation Models Generated by Uav Photogrammetry for Slope Stability Assessment (Case Study of Landslide in Svätý Anton, Slovakia). Acta Scientiarum Polonorum Formatio Circumiectus 2016, 15, 439–450. [Google Scholar] [CrossRef]
18. Ariza López, F.J.; Atkinson Gordo, A.D. Analysis of Some Positional Accuracy Assessment Methodologies. Journal of Surveying Engineering 2008, 134, 45–54. [Google Scholar] [CrossRef]
19. Sanz-Ablanedo, E.; Chandler, J.H.; Rodríguez-Pérez, J.R.; Ordóñez, C. Accuracy of Unmanned Aerial Vehicle (UAV) and SfM Photogrammetry Survey as a Function of the Number and Location of Ground Control Points Used. Remote Sensing 2018, 10, 1606. [Google Scholar] [CrossRef]
20. Tzanakaki, K. Άπτερα Το Ταφικό Μνημείο ΙΙ: Aνασκαφικά Δεδομένα Και Ερμηνευτική Προσέγγιση. H Ελεύθερνα, η Κρήτη και ο Έξω Κόσμος/ Eleutherna, Crete and the Outside World. Πρακτικά Διεθνούς Aρχαιολογικού Συνεδρίου (Ρέθυμνο 2018) Ν. Χρ. Σταμπολίδης, Μ. Γιαννοπούλου (επιμ.) 2020.