Photogrammetry and outputs

Why this matters

Modern mapping software produces attractive results at the press of a button. But when the map splits into two blocks, a stockpile has holes, or the user asks for “ground height” under trees, someone who does not understand how images become maps cannot fix the problem, and may deliver files with the right names but the wrong meaning. This module explains as much of the mechanism as you need to make those decisions.

Depth comes from parallax

One photograph gives the direction from the camera to an object but not its distance. When the camera moves, the same point appears at a different position in the image. This is parallax. Hold a finger in front of your face and close each eye in turn: the nearby finger jumps more than distant objects. Photogrammetry uses the same principle.

For an image pair taken at the same height a distance $B$ apart (the base), the photogrammetry textbook by Wolf and colleagues gives $Z=Bf/p$, where $p$ is the parallax on the image plane.

BASE_M, FOCAL_MM = 16.0, 8.8
PIXEL_MM = 13.2 / 5472

def depth_m(parallax_px):
    return BASE_M * FOCAL_MM / (parallax_px * PIXEL_MM)

ground_px, roof_px = 729.6, 810.7
z_ground, z_roof = depth_m(ground_px), depth_m(roof_px)
print(f"ground {z_ground:.1f} m below camera, roof {z_roof:.1f} m, roof height {z_ground - z_roof:.1f} m")
print(f"1 px parallax error at 80 m changes depth by about {80 ** 2 / (BASE_M * FOCAL_MM) * PIXEL_MM:.2f} m")

ground 80.0 m below camera, roof 72.0 m, roof height 8.0 m
1 px parallax error at 80 m changes depth by about 0.11 m

Structure from Motion

A drone survey uses hundreds of images, not one pair. Structure from Motion (SfM) estimates the position and orientation of every camera together with the shape of the scene. The method grew out of computer vision work such as Snavely and colleagues (2006) and has been widely adopted in the geosciences, as reviewed by Westoby and colleagues (2012).

1. Images and camera model: read the image size, focal length and lens distortion
2. Detect and match features: for example, corners of marks on the ground. Points tracked across several images are tie points; matches that conflict with the geometry are rejected
3. SfM and bundle adjustment: OpenSfM starts from an initial image pair, adds images one at a time, computes 3D points and periodically refines everything together with bundle adjustment, so that points projected back into the images land as close as possible to where they are seen
4. Dense matching: estimate depth for many pixels, giving a dense point cloud
5. Georeference and outputs: use camera positions or GCPs, then create the orthomosaic, DSM or mesh

Reprojection error measures how many pixels a point projected back into an image lies from where it is seen. It is internal consistency, not error on the ground. A model with low reprojection error can still be shifted as a whole if the georeferencing is wrong.

Five kinds of output

OpenDroneMap exports the orthophoto and DEMs as GeoTIFF, the point cloud as LAZ and the mesh as OBJ; DSM and DTM must be switched on as options.

An RGB camera only sees what reflected light reaches, so under dense vegetation there may be no image of the ground at all. LiDAR measures distance from the travel time of laser light and may obtain ground points through gaps in the canopy, but it needs good calibration and navigation. The USGS sets requirements for LiDAR data in its Lidar Base Specification. Point cloud files from the two methods may look alike, but their evidence and limits differ.

Stockpile volume

Volume is computed from a height grid: each cell has area $A_i$ and a height equal to the surface minus the base. The result depends heavily on the chosen base level, and cells with no data must be reported, not counted as zero.

CELL_AREA = 2.0 * 2.0          # m² per cell
dsm = [  # surface height (m), None = no data
    [100.0, 101.0, 102.0, 100.0],
    [101.0, 103.0, 104.0, 101.0],
    [100.0, 102.0, None, 100.0],
]

def volume(base_level):
    cells = [z for row in dsm for z in row if z is not None]
    return sum(CELL_AREA * max(z - base_level, 0.0) for z in cells)

missing = sum(z is None for row in dsm for z in row)
for base in (100.0, 100.5):
    print(f"base {base} m: volume {volume(base):.0f} m3")
print(f"cells without data: {missing} of {sum(len(r) for r in dsm)} -> report, do not treat as zero")

base 100.0 m: volume 56 m3
base 100.5 m: volume 42 m3
cells without data: 1 of 12 -> report, do not treat as zero

Module lab

Common mistakes

Summary

• Depth comes from parallax when the camera moves; heights are sensitive to parallax error
• SfM estimates cameras and shape together; bundle adjustment minimises reprojection differences, and reprojection error is internal consistency
• Choose outputs to match the question: the DSM is the top surface, the DTM is the ground
• Volume depends on the base level, and areas without data must be disclosed

Check your understanding

1. $B$ = 20 m, $f$ = 8.8 mm and the parallax is 2.2 mm. What is the depth?
2. A reprojection error of 0.8 pixels means the map is within 0.8 pixels on the ground. True or false?
3. A user wants to design drainage in an area with trees. Should they use the DSM or the DTM, and what should they watch for?
4. A grid of 5 cells, each 4 m², is 1, 1, 2, 2 and 3 m above the base. What is the volume?
5. Why must no-data cells not be counted as zero?