ORTHORECTIFICATION OF A LARGE DATASET OF HISTORICAL AERIAL IMAGES: PROCEDURE AND PRECISION ASSESSMENT IN AN OPEN SOURCE ENVIRONMENT

: The availability of data time series spanning a long period is crucial for landscape change analysis. A suitable dataset, both in terms of time span and information content, must be available for the use with a GIS. In Italy, one of the most important historical source of land cover analysis is the GAI (Gruppo Aereo Italiano) photogrammetric survey (”Volo GAI”) commissioned in 1954 by the Italian national mapping agency, Istituto Geograﬁco Militare Italiano (IGMI). The survey covers the whole Italy, but so far only some Regions, namely Lombardia and Veneto, have carried out the image rectiﬁcation and the successive analyses to map land cover and use. This work describes the process of image orthorectiﬁcation of the Volo GAI images for the Province of Trento (Provincia Autonoma di Trento). Image orthorectiﬁcation must be performed to transform the images in maps available for analysis. This procedure corrects the geometry according to the terrain surface described by a Digital Terrain Model (DTM) to create an image compatible with the cartographic projection in use. To this end, the orthorectiﬁcation modules available in GRASS GIS have been used, with the advantage of using the same GIS environment which will be used for the landscape analysis. The dataset covering the whole Province contains almost 100 images, this paper presents the preliminary results of the orthorectiﬁcation of a quarter of the images. A reduced dataset has been used to test the results obtained using different settings with respect to: digital image resolution, DTM resolution and number of Ground Control Points (GCPs) used for the external orientation. These preliminary tests show that for the average quality of the Volo GAI images scan resolution beyond 600 DPI and DTM resolution above 10m do not provide signiﬁcant improvements for orthorectiﬁcation images. The minimum number of GCPs to guarantee the requested accuracy can vary from image to image, depending on the image quality and recognizable features position, but it is usually in the 15-20 points range.


INTRODUCTION 1.1 A brief introduction to photogrammetry
Historic data series of aerial photographs represent nowadays one of the most important source of geographical information (Rocchini et al., 2012) to perform multitemporal comparison of land use change and therefore detect landscape changes (Tattoni et al., 2010), (Ciolli et al., 2012), (Tattoni et al., 2017).However, aerial images present the same problem as topographic maps: the process of embedding Earth surface into a plan introduces distortions, because we are forcing a rugged and curved surface to fit on a two dimension flat plan; furthermore, photographs are sensible to the atmospheric conditions and light exposure (Cohen et al., 1996).Geometric errors alter the perceived position and size of some feature: tilt displacement is due to inclination of the focal plan with respect to the ground and relief displacement is due a change in elevation (Paine and Kiser, 2003).Figure 1 shows a clear example these distortions: the road in the non processed image (left) appears to be crooked, while is actually straight (processed image on the right) (OSSIM Development Team, 2015).Figure 1.A comparison between a distorted image (left) and the same, rectified image (right).The road appears to be crooked on the left image, while in reality it is straight (right, processed image), courtesy from OSSIM Wiki.(OSSIM Development Team, 2015) The International Archives of the Photogrammetry, Remote Sensing and Spatial Information Sciences, Volume XLII-4/W8, 2018FOSS4G 2018 -Academic Track, 29-31 August 2018, Dar es Salaam, Tanzania Radiometric errors are distortions in tone or colour of the photographs.They can be caused by vantage points, condition and calibration of the camera or by the post-processing of the image, during the filtering and film emulsion (Cohen et al., 1996) (Jensen, 1996).Photogrammetry is primary concerned in obtaining precise results from aerial photographs: in particular the process of compensating both geometric and radiometric errors is called orthorectification (Morgan et al., 2010).In this paper we will explain briefly elaboration of a large set of photographs, presenting a procedure implemented in GRASS GIS (Geographical Resources Analysis Support System, Geographic Information System) the Free and Open Source Software used for the procedure of orthorectification (GRASS Development Team, 2015).
The frames used are part of the "Volo GAI", a photogrammetric survey commissioned in 1954.In particular we will focus on the survey carried out in Provincia Autonoma di Trento (PAT, a north eastern Italian region of 6.200 km 2 ) that needed rectification before being used in GIS software.

The dataset of 1954 Gruppo Aereo Italiano (GAI) aerial imagery
After World War II the Istituto Geografico Militare Italiano (IGMI) commissioned a series of flights over the whole Italy to survey the territory.The importance of this dataset lies in the fact that this is the first stereoscopic photographs depiction of the Italian region.
The height of the flight was around 10.000 m a.s.l. for the mountainous areas, with a mean scale of 1:45.000(Geoportale Regione Lombardia, 2018).Each photograph is square shaped with a dimension of 230 x 230 mm and covers an area of about 140 km 2 (IGMI, 2016).
Currently IGMI still holds the rights for original data, but every Region of Italy acquired the digital version of these frames and it is possible for everyone interested to access the data (IGMI, 2016).Figure 1   The set of photographs used for the current work covers the entire Provincia Autonoma di Trento.All the photographs present very little or no cloud coverage.Shadowy areas are a major drawback: this is because lot of photographs were taken during the first hours of sunlight, when Sun is lower on the horizon and shadows are longer in mountainous areas.
The physical copy has been digitalized through scanner taking into account the final file size of the scanned image and its final mean ground resolution.The equation 1 was used to calculate ground resolution:

Orthorectification
Orthorectification is the process of modifying the geometry of an image to make it compatible with a carthographic projection.
This procedure is performed in three steps: internal orientation to evaluate the position of the image with respect to the camera; external orientation to evaluate the position of the camera with respect to the external reference system; orthorectification to re-project the image.
The first two steps use sets of equations whose parameters must be evaluated using external information, usually the coordinates of points recognizable on the image.
To perform the internal orientation the identification of the position of 4 or more fiducial markers on the original photograph is required.
A set of Ground Control Points (GCPs) is necessary to perform the external orientation.These GCPs are points of geometrical features on the image whose coordinates are known in a reference system.Novak (Novak, 1992) reviewed the most common used methods for orthorectification.Polynomial rectification consists in transforming the original image into the orthorectified one by means of polynomial equations: where x, y = original image coordinates; x , y = rectified image coordinates; A,B = coefficient matrices of the polynomials.
The International Archives of the Photogrammetry, Remote Sensing and Spatial Information Sciences, Volume XLII-4/W8, 2018FOSS4G 2018 -Academic Track, 29-31 August 2018, Dar es Salaam, Tanzania The maximum order of these polynomials depends on how many GCPs are available.
Differential rectification is based upon a set of equations involving information about geometry of the camera sensor and the Digital Terrain Model (DTM) to correct for relief displacement.This equations, called collinearity equations, rectify the original image by shifting, rotating and scaling each of its pixel (Novak, 1992): where c = focal length of the used camera; x, y = transformed image coordinates; X0, Y0, Z0 = coordinates of projection centre X, Y, Z = object coordinates of the original image; rxx = orientation of the camera (components of the rotation matrix).
Usually the calibration certificate of the camera provides information about coordinates of the projection centre and focal length (Morgan et al., 2010).The algorithm implemented in GRASS creates a grid on output orthorectified image and uses equations 3 to find the corresponding point on the original image.Every pixel of the grid is assigned a grey value (or color), using different interpolating methods (Novak, 1992): 1. starting from the original pixel location it transforms its coordinates and assigns the grey value to the nearest pixel (nearest neighbour); 2. starting from the original pixel location it transforms its coordinates and assigns the grey value interpolating the values of nearest pixels using bilinear methods; 3. starting from the original pixel location it transforms its coordinates and assigns the grey value interpolating the values of nearest pixels using cubic convolutions.

GRASS GIS procedure for orthorectification
The choice of a software based upon Free and Open Source philosophy (such as GRASS) allows the final user to access the source code behind every module (Neteler and Mitasova, 2008).It is worth noting that full access to the program code represents a major advantage with respect to a proprietary solution: the possibility of checking every step of an algorithm is useful to guarantee the appropriate robustness of the output (Rocchini et al., 2012).Futhermore, it is possibile to modify the procedures and the algorithm according to the user requirements, provided that personnel with sufficient programming skills is available (Ciolli et al., 2017).GRASS GIS version 7.4 includes a suite of commands for orthorectification similar to the one already implemented in the previous 6.4 version.
The main body of the suite is divided in modules, which can be invoked singularly or in sequence: i.ortho.targetThis module sets a target location for the image to be rectified.The target location is a subset of maps in GRASS containing the maps used during the procedure (DTM and maps for GCPs) and implicitly defines the output orthoimage reference systems and projection.
i.ortho.elevThis module allows the user to select a DTM for the ortho-rectification process.In the current study different DTMs with different resolutions (1m, 5m, 10m and 25m) have been used for testing purpose.
i.ortho.cameraThis module sets the parameters used for the internal orientation of the image.The focal length of the camera and the coordinates of the fiducial marks on the image must be provided.Since the calibration certificate of the camera is not available for the current dataset, the coordinates of have been determined measuring the distance between the marks on the physical copy of the photographs.These measures revealed a distance of 232 mm between east-west markers and 233 mm between north-south markers.The internal orientation is performed using a 6 parameters affine transformation, therefore the coordinates of at least 3 points must be provided.
i.ortho.initThis module can be used to provide initial values (and their RMS) for the unknowns in the iterative least square adjustment for the evaluation of the external orientation parameters.The need to provide initial values of the 7 unknown parameters, representing 3 translations, 3 rotations and a scale factor, is due to the non linearity of the equations.Initial approximate values are automatically calculated If no value is provided.

g.gui.photo2image
This module provides a graphical interface to indicate the position of the fiducial markers on the digital image.
g.gui.image2targetThis module provides a graphical interface to locate GCPs on the image and provide their coordinates.It is possible to either manually enter the coordinates or locate the same point on a (geo-referenced) map.The height is provided by the DTM.
i.ortho.rectifyThis module performs the actual rectification of the image.

TESTS
A series of test have been carried out to assess the influence of the features of the input data and of the GCPs choice on the geometric proprieties of the output orthophoto.
The main aim of these tests is to determine whether increasing the resolution of the image and/or of the DTM provides a significant advantage in terms of geometric accuracy of the features on the orto-rectified image.Given the large number of images, the possibility of using a reduced resolution reduces considerably the size of the dataset.The feasibility of the use of a medium resolution DTM makes possible to perform the orthorectification procedure even where a high resolution DTM is not available.
An example of the relationship between digital image resolution, file size and ground resolution is given in The analysis of the precision of the external image orientation is based on the root mean square error (RMSE) for each GCP (equation 4).
where RM SEi = Root Mean Square Error of the i-th GCP; µi = residual of the i-th GCP in the east-west direction; νi = residual of the i-th GCP in the north-south direction.
These RMSEs depend on the accuracy of the GCPs selection and on their correspondence the model used in the Least Square Adjustment for the parameters estimation (Rocchini et al., 2012).
An "overall RMSE" is defined as follows: where RM SEi = Root Mean Square Error of the i-th GCP; N = Total number of GCPs.
Tests on the influence of different image and DTM resolutions on the results have been carried out for image 4954, which is representative for the datasets in terms of coverage and height variation.
Using a 1200 DPI scan obviously improves the mean ground resolution of the non rectified image and the overall precision in the rectified image at the cost of larger disk space occupation (table 1), higher memory demand during the orthorectification process and longer processing times.Table 2 shows that with the same number and positions of GCPs the overall RMS is basically constant in terms of pixels, therefore it decreases with the same ratio as the resolution in meters.The DTM resolution effect on the overall RMS has been tested by performing the orthorectification procedure on the same image with the same resolution and the same GCPs, varying only the DTM resolution.Differences between using a 1m resolution DTM and a 10m resolution DTM hardly exceeds 1 pixel, so the optimal choice in terms of memory allocation and processing time for this set of photographs (at a resolution of 600 DPI) is 10m DTM.
These results are only a small example of the overall RMS for the large set of rectified images, but it is representive the vast majority of the results obtained.

RESULTS
A grand total of 90 photographs from volo GAI have been rectified so far, using a 600 DPI resolution and a DTM of 10m resolution, following the indications emerged from the tests described in section 3.
The nearest neighbour method was chosen as interpolation method, because it is less computationally intense compared to the other two available methods (Novak, 1992).The processing has been carried out using the WGS84 reference system with UTM projection, zone 32 N.The quality of the original images varies frame to frame, while some images easily allow the recognition of features such as trees or isolated structure (figure 5), in many images urban areas appear as a white blur, with poor detail (figure 6).Furthermore, another issue arises in areas with high slopes.There the DTM acts like a discontinuous function and the algorithm is not able to fit precisely the photograph to the DTM.This causes a stretch and a step pattern in the rectified image.Figure 7 shows a comparison between the step pattern obtained with different resolutions DTM: using a finer resolution causes a decrease in size of the pattern.
Figure 7. Effects of high slopes in rectified photographs.For the same area on the left a 1 m DTM was used, on the right a 10 m DTM was used.

Precision assessment
Once image and DTM resolution is set, the only variable influencing the precision of orthorectification are number and positions of GCPs.
Their positioning strongly depends upon the availability of recognizable features on the territory depicted in photographs.Ideally, GCPs should be distributed homogeneously through the photograph to assure homogeneous results with good results for every part of the image.On this basis, Bernstein proposed that a sufficient number of GCPs should be 16 (Bernstein and Colwell, 1983).However, it is not always possible to place 16 GCPs homogeneously.For example, some photographs included in Volo GAI depict the glacier zone of "Parco dell'Adamello-Brenta", which were completely occupied by glaciers in 1954 and have changed dramatically.Furthermore, the mean altitude of these areas is above 2000 m a.s.l.where there are few buildings, so features that can be used as GCPs are very rare and not homogeneously distributed.
In such cases orthorectification is still possible: GRASS GIS algorithm can orthorectify an image with 4 non-aligned points.The result needs to be carefully inspected: with g.gui.mapswipetool in GRASS GIS it is possible to compare features (such as roads, lake shores and mountain ridges) between the processed photograph and a reference rectified image or map to check of differences.Although it is possible that they have changed shape through the years when using maps with different dates, an error in the rectification process is often recognizable, as figure 8 shows.et al., 2010), (Maimeri, 2018)).In both cases a land use classification was carried out to extract the forested and pasture areas; these were then resampled with a pixel resolution of 10m.
In those studies the required precision was in the same order of magnitude as the RMSE of orthophotos.For other kinds of studies this precision could not be enough: before using 1954 Volo GAI rectified imagery it is important to manually check both precision and quality to avoid bias in the results.
5.2 Discrepancies between GRASS 7.4 and 6.4 GRASS 6.4 and 7.4 both implement apparently the same version of the orthorectification algorithm.However, using the same in-puts (image, GCPs, fiducial markers and DEM) on the same computer can lead to two different results.Figure 9 is obtained using the same photograph (3862) and same inputs in GRASS 7.4 and 6.4.Using r.mapcalc the pixel by pixel grey level difference have been evaluated for the two images.The absolute differences between the two output orthophoto is shown in figure 9, in a grey scale representative with lower values in black and higher values in white.A correlation of the higher values with the slope values can be observed, but it is not clear why these discrepancies occur: further analysis are required on versions 7.4 and 6.4 of GRASS algorithm implemented in the i.ortho.photomodule.

Conclusions
As stated in (Rocchini et al., 2012) differential rectification in a region characterized by a very diverse and ragged territory has proven its potential in georeferencing a large set of photographs, obtaining precision on RMSE around 10m for each photograph.In mountainous areas using just a polynomial rectification could lead to read tilt and displacement errors which will affect analysis based on remote sensing and image classification (Aronoff, 2005) (Morgan et al., 2010).
Results of section 3 show that an increase in image and DTM resolutions does not provide significant advantages for the orthorectification process with the 1954 Volo GAI dataset.Other limitations due to the original image quality play a more important role in the information that can be obtained from the orthorectified images, as seen in section 4.
GRASS GIS has proven to be a proper solution for rectification: 1. it is a free and open software, for many researchers it allows collaboration and research reproducibility; The International Archives of the Photogrammetry, Remote Sensing and Spatial Information Sciences, Volume XLII-4/W8, 2018FOSS4G 2018 -Academic Track, 29-31 August 2018, Dar es Salaam, Tanzania 2. it represents a way to save costs and still work with a powerful GIS tool capable of orthorectifiyng; 3. its code could be edited to give better performances in orthorectification; 4. it is user-friendly and troubleshooting is made easier by an active community of developers.
Future studies will investigate the differences between GRASS 6.4 and 7.4 results examining the source code to understand the different results.Further analysis will be conducted on the effects of the different available interpolating methods (nearest neighbor, bilinear and bicubic interpolation) during the orthorectification process.
However, GRASS GIS is not the only open source solution, as other FOSS systems, such as OSSIM (OSSIM Development Team, 2015), can also perform orthorectification.
shows an example of how frames are arranged.The fiducial markers of the frame, used for the internal image orientation, are highlighted in green.Parameters of the camera and of the flight era highlighted in blue: name of the camera, focal length of the camera lens, flight height, date and hour.The progressive number of the image in the survey is highlighted in red.

Figure 2 .
Figure 2. A non rectified photograph from 1954 GAI flight.For the meaning of highlighted areas we refer to paragraph 2.1.
Comparison of overall RMS for image 4954 with different scan and DTM resolutions.

Figure 3
Figure3shows the tile index for the orthorectified images.They cover the vast majority of Trentino territory, while images for some other areas, not covered by the tiles in figure3, namely Monte Baldo, Paneveggio, Tesino and Val di Fassa have been already rectified in previous works(Rocchini et al., 2012) (Tattoni et al., 2010) (Cattani, 2015) (Maimeri, 2018).

Figure 3 .
Figure 3. Tile coverage of the orthorectified photographs.

Figure 4
Figure 4 shows one rectified photograph (3853): the displacement effect of the altitude variations is clearly visibile in the south (bottom) and nort east (top right) part of the image.The image is also rotated counter-clockwise due to the difference of direction between the flight path and the datum axes.

Figure 5 .
Figure 5.A pasture near Campitello (Val di Fassa): buildings, trees and open areas are clearly recognizable.

Figure 6 .
Figure 6.Depiction of the village of Campitello (Val di Fassa).It is not clear how buildings and streets are arranged.

Figure 8 .
Figure 8.A comparison between the same lake in a 1954 rectified photograph with few GCPs (right) with an already rectified 2006 image (left).Note the displacement of the edges of this particular shape.

Figure 9 .
Figure 9. Discrepancies between the same photograph (3862) rectified in both 7.4 and 6.4 version of GRASS.In white are represented pixels with different grey values, in black pixels with the same grey value.

Table 1 .
File size and ground resolution for 600 DPI and 1200 DPI resolutions 23x23cm images.

Table 3 .
Variation of the precision in photograph 4954 for different choices of DTM and scanning of 600 DPI.
Table 3 shows the resulting overall RMSEs using different DTMs on image 4954 with a 600 DPI resolution.