Using image-based modelling (SfM–MVS) to produce a 1935 ortho-mosaic of the Ethiopian highlands

Abstract

Approximately 34,000 aerial photographs covering large parts of Ethiopia and dating back to 1935–1941 have been recently recovered. These allow investigating environmental dynamics for a past period that until now is only accessible from terrestrial photographs or narratives. As the archive consists of both oblique and vertical aerial photographs that cover rather small areas, methods of image-based modelling were used to orthorectify the images. In this study, 9 vertical and 18 low oblique aerial photographs were processed as an ortho-mosaic, covering an area of 25 km², west of Wukro town in northern Ethiopia. Using 15 control points (derived from Google Earth), a Root Means Square Error of 28.5 m in X 35.4 m in Y were achieved. These values can be viewed as optimal, given the relatively low resolution and poor quality of the imagery, the lack of metadata, the geometric quality of the Google Earth imagery and the recording characteristics. Land use remained largely similar since 1936, with large parts of the land being used as cropland or extensive grazing areas. Most remarkable changes are the strong expansion of the settlements as well as land management improvements. In a larger effort, ortho-mosaics covering large parts of Ethiopia in 1935–1941 will be produced.

Keywords:

• land use
• oblique aerial photograph
• ortho-mosaic
• PhotoScan
• Structure from Motion–MultiView Stereo

1. Introduction

The aerial photographs of Ethiopia taken by the Italian Military Geographical Institute (IGMI) between 1935 and 1941 are the oldest remote sensing data of the country. They cover large parts of the central and northern highlands and form a unique dataset that allows the study of the status of the environment, in a period that until now is only accessible from terrestrial photographs (Nyssen et al. 2009; Frankl et al. 2011) or narratives (Pankhurst 1997). In total, ca. 34,000 individual aerial photographs were registered. Until recently, they were stored in old ammunition boxes of the Italian army which, together with age, strongly reduced their quality. Using multi-cameras, one vertical photograph, two low oblique photographs (28°30′ to the vertical) and one high oblique photograph (57°00′ to the vertical) were exposed simultaneously across the flight line (Figure 1; Nyssen and Petrie 2013).

Figure 1. Flight plan of the nine sets of aerial photographs used to produce the ortho-mosaic. The high-oblique photographs (shown as dashed lines) could not be incorporated in the production of the ortho-mosaic.

Stereographic analysis of the aerial photographs and subsequent mapping are originally made for the use of manual and analytical photogrammetric instruments that are no longer in use. The Santoni Stereocartograph and Stereosimplex instruments allowed corrections for the angular values of the oblique cameras when producing the stereo models (Nyssen and Petrie 2013). Producing orthophotographs from this material using conventional photogrammetric software did not deliver satisfactory results, as (1) stereo pairs of the vertical photographs cover small areas (ca. 4 km²) and (2) it was not possible to produce stereo models from overlapping oblique photographs. Moreover, fiducial marks were not present on the aerial photographs and no information was available on the camera calibration parameters.

Image-based modelling techniques allow textured 3D models to be created from large sets of (both terrestrial and airborne) photographs of the same surface taken from different positions. The methodology is based on semi-automated Structure from Motion–MultiView Stereo (SfM–MVS) workflows (Seitz et al. 2006; Verhoeven et al. 2013), which are integrated in software like PhotoScan (Agisoft). In recent years, detailed and accurate 3D photo models were created for small cultural heritage objects (Stal et al. 2012), buildings (Stal et al. 2011) and landforms (James and Robson 2012; Frankl et al. 2014). Using 20 vertical aerial photographs from 1945 over the Kreuttal region of Austria, Verhoeven et al. (2013) produced an historical ortho-mosaic. In a similar way, Gomez (2012) produced an ortho-mosaic of the Iwaki volcano (Japan) using vertical aerial photographs of 1965–1994. This reveals that multiple oblique and vertical aerial photographs taken over the same surface – as the 1935–1941 dataset of Ethiopia – are potentially suitable to be processed using image-based modelling.

The objective of this paper is, therefore, to use image-based modelling as a methodology to produce ortho-mosaics from multiple overlapping vertical and oblique aerial photographs taken from Ethiopia in 1935–1941. These ortho-mosaics will allow spatially explicit studies on past land use and cover or geomorphic processes.

2. Materials and methods

The scale of the vertical photographs collected by IGMI in 1935 is approximately 1:14,000, while the low and the high oblique photographs have scales of 1:15,000 and 1:26,000 (valid at their principal point). Sequential exposures have a longitudinal overlap of ca. 60%. The ground coverage of one assemblage covers a strip of approximately 11.5 km across the flight line, of which 1 km is recorded by the central vertical photograph (Figure 1).

In this study, from one flight line, 9 vertical and 18 adjacent low oblique aerial photographs were used to create an ortho-mosaic (Figure 1). The high-oblique aerial photographs could not be processed by image-based modeling due to (1) insufficient overlap as successive high-oblique aerial photographs were taken alternatively from left to right, (2) the very large image deformations induced by the large tilt of the camera, (3) very large relief displacements in the image and (4) the occurrence of large ‘blind’ areas hidden behind imposing mountain ranges. In total, the set of 27 aerial photographs covers an area of 25 km². The area is located 10 km west of the town Wukro (North Ethiopia), a semi-arid limestone plateau dissected by the Suluh River valley. Elevations in the photographs scenery range between 1880 and 2140 m a.s.l.

First, the 27 aerial photographs (10 × 15 cm large) were scanned at 600 dpi with a desktop scanner (Nyssen and Petrie 2013). Then, the spatial extent of the aerial photographs was mapped on high-resolution Google Earth images (2006 Digital Globe, spatial resolution of ca. 0.5 m). Second, the individual vertical and low-oblique aerial photographs were pre-processed by optimising contrast and brightness in Adobe Photoshop. Overlapping vertical and low-oblique photographs were loaded in PhotoScan Pro, black borders where masked (in order to avoid erroneous matching) and the integrated semi-automatic SfM–MVS algorithm was run (Verhoeven et al. 2013). The SfM algorithm allows the creation of a 3D structure from multiple 2D images (i.e. aerial photographs). This included a feature point detection and description, image matching (cfr., defining tie-points), image triangulation and bundle adjustment (resulting in a sparse 3D point cloud of image matching points based on the projection of image pixels in 3D). As fiducial marks were not present on the aerial photographs and no information was available on the camera calibration parameters, a conventional photogrammetric inner orientation could not be performed. Hence, in PhotoScan, a separate camera model was estimated for each frame, neglecting the interior orientation but corresponding with the best-fitting linear geometry.

The absolute orientation was performed after the image alignment by assigning Ground Control Points (GCPs) to the images. These GCPs were defined by the unambiguous detection of characteristic points on the ancient images and Google Earth imagery (Frankl et al. 2013). The georeferencing of the model into WGS1984 was done by using points that occur on at least on two overlapping aerial photographs. Identifying GCPs on both Google Earth imagery and the aerial photographs was time-consuming, and depending on the area and image contrast, the number of points that could accurately be identified varied. In total, 76 GCPs could be identified that were well distributed in X, Y and Z. In order to determine the number of GCPs that needed to be selected to optimise the accuracy, models were produced with from 3 up to 60 GCPs, corresponding to a density of 0.12–2.4 points per km² (Figure 2A). The Root Means Square Error (RMSE) was computed for the GCPs, and the Mean Error (ME) and Standard Deviation of the Error (SDE) was computed for the remaining checkpoints (not used in the georeferencing). This follows the guidelines of Smith, Chandler, and Rose (2009), which define the ME as a good measure to detect systematic errors and the SDE as a measure for surface variability or precision.

Figure 2. Accuracy assessment of the image-based models in function of an increase in number of GCPs. (A) Density of GCPs in function of number of GCPs, (B) RMSE (m) in function of number of GCPs, (C) ME (m) in function of checkpoints and (D) ME (m) in function of checkpoints.

Finally, the MVS algorithm in PhotoScan Pro produced a dense estimate of the surface geometry in the form of a Triangulated Irregular Network (TIN), which enabled a georeferenced ortho-mosaic to be generated. To get acquainted with the reality on the ground, a field visit was organised on 27 November 2013.

3. Results and discussion

Selecting the optimal number of GCPs for the absolute orientation was achieved by identifying the number of GCPs for which an additional increase in GCPs did not significantly affect the values of RMSE, ME and SDE (Figure 2B–2D). For all accuracy indicators (RMSE, ME and SDE), increasing the number of GCPs from 3 to approximately 15 resulted in more varied errors. Considering the RMSE in X, Y and Z, the values stabilised at an error of ca. 30 m in X and Y and ca. 50 m in Z when the number of GCPs was 15 or more (Figure 2B). Using 15 GCPs or more also caused the ME in X and Y approximately zero, meaning that no systematic errors occur in the planimetry (Figure 2C). For Z, a systematic 30-m overestimation of the elevation persists. This could be related to the poor resolution (90 m) of the Shuttle Radar Topographic Mission elevation data compared to the 3D geometry reconstructed by the SfM method. As the primary objective of this study was to produce an ortho-mosaic (being a true representation of the planimetry), a more detailed error analysis in Z was not done. We hereby stress that the models are not suitable for extracting elevation data or to study changes therein. The SDE of the checkpoints is similar to the RSME of the GCPs, although in X and Y increasing the number of GCPs did not seem to have much effect on improving the model accuracy (Figure 2D). From this, we retained 15 GCPs to produce the TIN and ortho-mosaic, corresponding to a GCP density of 0.6 GCPs km^–2 (Figure 2A). Table 1 summarises the resulting RMSE, ME and SDE. Considering that the SfM–MVS method yields accuracies of 1/1000 of the viewing distance, the planimetric accuracy of our ortho-mosaic (approx. 30 m in X and Y) is rather low; with the flight height being 2400 m above the ground surface, the maximal accuracy would be 2.4 m. However, this study uses aerial photographs of >70 years old that were stored in suboptimal conditions. Given (1) the deformation of the aerial photographs, (2) the very small intersecting angles of the corresponding rays, (3) the poor quality of the images and the lack of metadata (EXIF files) and (4) the geometric quality of the Google Earth imagery, the obtained error statistics can be viewed as optimal. From this dataset, improving the accuracy might be obtained when using GCPs that were recorded in the field with accurate Global Navigation Satellite System (GNSS) receivers or when using higher-resolution Digital Elevation Models. However, as indicated by James and Varley (2012), using planimetric information from digital globes can yield accurate results. In addition, using GNSS receivers would limit the methodology to local studies, making the production of a country-wide 1935–1941 otho-mosaic unfeasible. Using original negatives (that include fiducial marks) may allow the calculations of an inner orientation of the aerial photographs and also improve the accuracy. These could, however, not be accessed for this study.

Table 1. Summarising statistics describing the accuracy of the retained photo-model with 15 GCPs.

	X (m)	Y (m)	Z (m)
RMSE	28.5	35.4	50.7
ME	−4.7	−3.1	30.3
STE	34.5	30.7	53.5

In order to determine if the accuracy of the models is uniform for the scenery covered by the vertical and low-oblique aerial photographs, RMSEs were computed using GCPs identified on vertical and low-oblique aerial photographs, respectively. As indicated by a paired t-test (α = 0.05), the average errors in XYZ were not significantly different between the vertical and the low-oblique aerial photographs (p = 0.92).

Figure 3 shows the ortho-mosaic of the Ethiopian highlands of 1936, based on 15 CPs and the corresponding Google Earth imagery of 2008. The ortho-mosaic can also be viewed on http://geoweb.ugent.be/physical-geography/research/environment-ethiopia. The white stripes on the ortho-mosaic are areas for which a 3D surface could not be produced by PhotoScan and correspond to narrow gaps in ground coverage between the vertical and the low-oblique aerial photographs of one set, most probably due to manual cropping of the original aerial photographs when preparing the assemblages. Areas for which image matching was sub-optimal resulted in local distortions in the ortho-mosaic. The large variability in brightness between the different aerial photographs is also visible in the ortho-mosaic. The central part, which mainly corresponds to the vertical photographs, is generally brighter than the areas of north and south covered by the low-oblique aerial photographs.

Figure 3. (A) 1936 ortho-mosaic (arrow indicates areas with insufficient overlap for ortho-mosaic creation) and (B) 2006 Google Earth imagery for the same area (with the Suluh River mapped in blue). (a)–(d) are explained in the text and in Figures 4 and 5. The ortho-mosaic can also be viewed dynamically at http://geoweb.ugent.be/physical-geography/research/environment-ethiopia.

When comparing land use and management between the two periods, we can observe that the land use remained largely similar between 1936 and 2008 (Figure 3). Large parts of the land were used as cropland and extensive grazing areas of sparse shrubland existed on steep slopes or where soils were too shallow to be cultivated. The areal expansion of the village of Zahina (a) is among the most remarkable changes between 1936 and 2008 (Figures 3 and 4). Moreover, the village area includes many more trees nowadays than ca. 80 years ago. The steep slopes that drain towards the Suluh River were incised by gullies in 1936 (b), and were covered by sparse shrubland and cropland. By 2008, land management on these steep slopes strongly improved with the introduction of stone bunds along contour lines (c) indicates a remnant forest which had the same extent in 1936 as in 2008. Preservation of the forest is due to the existence of the ‘May Cholot’ spring in the forest, which is considered as holy water by local communities. The vegetation cover in the Suluh River valley (d) did not change dramatically, although eucalypt plantations were introduced on the floodplain. The 2013 situation of (a), (b), (c) and (d) can also be viewed on terrestrial photographs in Figure 5.

Figure 4. Detailed view of Zahina village in 1936 (left) and 2006 (right, black and white Google Earth imagery) indicating the large increase in housing and trees. Although the situation has changed drastically, some features remained similar, as the cross-point of footpaths in both periods (black arrows). The red angle shown in the 2006 situation corresponds to the viewpoint of Figure 5(a).

Figure 5. Views of the mapped area in 2013. (A) Zahina village, (B) steep slope with stone bunds draining to the Suluh River, (C) remnant forest and (D) the Suluh River valley.

In order to develop a methodology for producing an ortho-mosaic, only 18 aerial photographs were used in this study. However, depending on the hardware used, ortho-mosaics can be produced for much larger areas, using more photographs. The maximum number of photographs that can be used mainly depends on the hardware. Studies are known that use more than 1000 photographs to produce image-based models (Castillo et al. 2013), indicating the potential for producing ortho-maps from the ca. 25,000 low-oblique and vertical aerial photographs of Ethiopia collected in the second half of the 1930s.

4. Software

PhotoScan recommends a hardware configuration of: CPU: Intel Sandy Bridge 2600 or 2700 RAM: DDR3-1600, 4 × 4 GB (16 GB total) or 4 × 8 GB (32 GB total) GPU: NVidia GeForce GTX 580. More information can be found on http://www.agisoft.ru/support. System requirements for Google Earth can be found on https://support.google.com/earth/?hl=en. For OS, Windows XP or later is required and for Mac, Mac OS X 10.6 or later.

5. Conclusions

In this study, multiple vertical and oblique aerial photographs taken from Ethiopia between 1935 and 1941 were used to produce a historical ortho-mosaic. This was done from image-based modelling methodologies integrated in PhotoScan. As a result, the Ethiopian landscape can be analysed from orthorectified large-scale aerial photographs, taking into account an planar error of ca. 30 m. Characteristics of historical land use and cover, geomorphological processes, etc., which until now could only be studied from archival terrestrial photographs or narratives, can now be studied from ortho-mosaics much in the same way as it is done for more recent periods. In a larger effort, multiple ortho-mosaics will be merged in order to produce an ortho-map of the early twentieth-century Ethiopian landscape.

Acknowledgements

The 1935–1941 aerial photographs were retrieved and scanned through a Memorandum of Understanding between the Ethiopian Mapping Agency (EMA), Mekelle University (Ethiopia) and Ghent University (Belgium); we particularly extend our thanks to Mr. Sultan Mohamed, Mr. Gezahegne Gebremeskel, Ms. Yetemwork Dejene and Ms. Shewaye (EMA), to Dr. Kiros Meles and Prof. Mitiku Haile (Mekelle University) and to Mr. Martijn Debever (Ghent University). Special thanks also to the people of Zahina village for their hospitality and to Seifu Gebreselassie and Dawit Lsanwerk for their support in the field.