From Archived Historical Aerial Imagery to Informative Orthophotos: A Framework for Retrieving the Past in Long-Term Socioecological Research

Aerial photographs have been systematically collected from as early as the 1930s, providing a unique resource to describe changes in vegetation and land cover over extended periods of time. However, their use is often limited by technical constraints, such as the lack of ground control information and precise camera parameters, which hamper an accurate orthorectification of the raw imagery. Here, we describe the historical aerial photographs orthorectification (HAPO) workflow, based on a conventional photogrammetric procedure (the direct linear transformation (DLT) Method), integrated as a geographic information systems (GIS) procedure, in order to perform the image orientation and orthorectification, thereby converting historical aerial imagery into high-definition historical orthoimages. HAPO implementation is illustrated with an application to a rugged landscape in Portugal, where we aimed to produce land-cover maps using an aerial photograph coverage from 1947, as part of a study on long-term socioecological dynamics. We show that HAPO produces highly accurate orthoimages and discuss the wider usefulness of our framework in long-term socioecological research.


Introduction
Historical aerial photographs (HAP) are the oldest form of remotely sensed data, yielding a spectrum of valuable historical information on vegetation and land cover [1][2][3]. Such information has enormous potential to provide crucial baseline data for quantitative mapping and monitoring of landscape composition and heterogeneity over time [3][4][5][6][7][8][9], and to relate land-use history to changes in biodiversity and ecosystem services [2,[10][11][12]. Aerial photographs started to be produced from as early as the 1930s [1,9], whereas satellite imagery only became available in the 1970s and, at the beginning, had coarse spatial resolution [9,13], which limits their usefulness to analyze long-term landscape change. Despite their value, however, the potential of HAP remains clearly underexplored due in part to technical constraints [14], thus hindering a thorough understanding of the spatial and temporal patterns of land use and land cover, and their resulting environmental impacts, particularly during the pre-satellite era [12].
The use of HAP is challenging due to the need for accurately georeferencing the photos to a ground coordinate system, so they can be compared with other geographic datasets in geographic information systems (GIS) [2,[15][16][17][18]. In general, GIS users do not have easy access to photogrammetry software tools. Additionally, HAP orthorectification is a complex and difficult process since historical ground control information, calibration certificates, and precise or even basic knowledge of camera parameters are often not available [16][17][18][19]. Furthermore, processing HAP presents additional challenges in landscapes undergoing regime shift processes [2,[16][17][18][19], due to the dramatic land-cover changes that can occur over time [20], making it difficult to find adequate ground control features that can be georeferenced from information in modern orthoimages. Probably because of these problems, to our knowledge, there is still no fully automatic pipeline for generating orthoimages from HAP [14].
To address these challenges, some studies have exploited a relatively fast and semiautomatic process by using the implicit depth information contained in overlapping imagery, and by subsequently introducing the extracted digital surface model (DSM) information into an object-based image analysis (OBIA) process (e.g., [18,21]). However, automatic extraction of DSMs with image-matching algorithms can deliver undesirable results (spiky surface models), resulting from unmodeled distortions, errors in the interior orientation parameters, and/or low radiometric quality [19]. Moreover, historical aerial stereo pairs are often not available, due to deterioration by age or other circumstances, as relocations or accidents in the archive facilities (e.g., [15]), which can damage part of the collections. Further, this approach cannot be used in the case of isolated images or incomplete data sets.
To overcome these problems, we developed a new workflow providing accurate historical aerial photographs orthorectification (HAPO) and georeferencing, which is implemented through a program we called "Projective 3D" (PRJ3D). GIS software, such as ArcGIS or QGIS, implements the projective transformation as a 2D image georeferencing method, which is only suitable for aerial photographs of flat terrain. Our method brings the third dimension into the process, combining direct linear transformation (DLT) for an individual image orientation, and conventional photogrammetric orthorectification procedures with a digital elevation model (DEM) for resampling individual historical aerial imagery into high-definition historical orthoimages. The methodology brings the old DLT method, normally not used in photogrammetric software tools, into an image-georeferencing process, overcoming the limitations found by GIS users in dealing with HAP. In this study, we describe the HAPO workflow, detailing all the steps required to process isolated or marginally overlapping sets of HAP. Workflow implementation is illustrated with a case study in Portugal, where we aimed to produce land-cover maps using aerial photograph coverage from 1947. The study shows the potential of HAPO to solve problems commonly encountered when orthorectifying historical aerial photographs and underlines its potential wide applicability in similar studies addressing long-term landscape dynamics.

Background
Aerial photographs have been systematically collected over decades for producing military and civilian maps and are archived in many institutions worldwide [15,22]. This resource has been used in a number of scientific fields, and its importance has increased in recent years to address questions involving long-term changes in socioecological systems. This is shown, for instance, by a simple literature search performed in ISI Web of Science (in December 2018), covering the period 1990-2018, where we used the search string TS = ("historical aerial photographs*" OR "Orthorectification* historical aerial photographs*" OR "historical imagery*" OR "orthorectification* historical imagery*" OR "historical photography*" OR "orthorectification* historical photography*"). The 260 papers retrieved by this search were spread over different research topics, such as geosciences (78 papers), geography (61), remote sensing (24), and environmental sciences/ecology (97), and there was a steady increase in the number of papers published since the 1990s (19) until the 2010s (241). Most studies using HAP in environmental/ecology applications were carried out in the USA (50%), while European studies represented a smaller proportion (29%).
These previous studies have used several different approaches to extract information from HAP. Some studies used a simple approach involving only HAP relative image orientation rather than georeferencing and orthorectification [15,19,[23][24][25][26], which provides a quick solution when a coarse spatial resolution is sufficient to detect changes [27,28]. These studies generally addressed stable areas over time, usually urban or man-constructed features [2], with analyses of change usually carried out over small spatial extents or at low resolutions. Other studies retrieved from the literature search used previously orthorectified (e.g., [9,29]) or only georeferenced (e.g., [30][31][32]) images but provided little information on image processing. In other cases, the camera information (metadata) was complete and available (e.g., [33][34][35]), but often there was no reference to the orthorectification process itself (e.g., [36][37][38][39]). More detail is usually given in technical studies dealing specifically with remote sensing, with several authors describing feature-based registration methodologies for historical aerial image orientation and orthorectification, such as the use of TIL (time-invariant lines) (e.g., [17]) or OBIA processes (e.g., [21]). However, the application of a standard bundle adjustment [16,18,19,36] is not possible when only even or odd HAP along strips are available, because the marginal overlap is around 20% and not the standard 60% of stereo pairs.
Overall, there are still significant technical constraints to retrieve meaningful information from HAP. Such constraints are particularly challenging when HAP are to be used in long-term socioecological studies, for which coarse thematic (and spatial) resolution classification systems may be insufficient to describe and understand the patterns and processes of landscape change [40]. For instance, these studies often require the identification of small patches or region-specific habitats of particular ecological importance [41], and to map their distribution at the landscape scale is not possible when using simple approaches, such as relative HAP orientation. Therefore, methods are still needed to deal successfully with the challenges of accurate HAP orthorectification, and to retrieve fine-scale information essential for high-resolution, spatially explicit long-term socioecological reconstructions [3]. In this paper, we address the main challenges of HAP georeferencing and orthorectification, describing a new workflow (HAPO) that can be used in the absence of ground control information or camera calibration and can be applied to isolated HAP or when stereo pairs are unavailable.

Orthorectification: The HAPO Workflow
Although in a future step, the HAPO workflow may be integrated in a GIS software (an ArcGIS toolbox or a QGIS plugin), for the time being, it was implemented as C programming language program (called "Projective 3D"-PRJ3D) and compiled with the Microsoft Visual Studio 2017 C compiler. The executable runs in the MS Windows command line, version 10 and older, taking several input files, and performs an independent orientation of individual photos using a direct linear transformation (DLT) [42] approach. The DLT method does not require camera information [16,18] but requires a point-based registration method to determine the exterior orientation by the collinearity model [17].
The DLT model is expressed by the following equations [42]: where (C,R) are the column and row coordinates of a point in the image matrix, in pixel units; (X,Y,Z) are the terrain coordinates (easting, northing and height) of the point in some cartographic reference system; and (a 1 , a 2 , . . . , a 11 ) are the coefficients characteristic of the image. These coefficients mix the exterior orientation parameters, as well as camera parameters, and can be determined if appropriate ground control points (GCPS) are available. In fact, in order to extract high-quality data from archival aerial photography, where there may be only a small number of images available and the redundancy may be low, orthorectification should be based on a reliable high-quality network of ground control information [15][16][17] well distributed in the block, especially when camera calibration information is unavailable [43]. The PRJ3D tool deals with the geometric aspects of the image orientation, to read image and terrain coordinates collected for a given photo, and to apply the least squares method to solve for the 11 unknown coefficients. For a given GCP, all coordinates will be available, so we can write 2 equations with the 11 unknown coefficients. A minimum of 6 noncollinear points would be needed to determine all the coefficients, by solving a system of equations. However, to obtain a more accurate solution, while controlling for errors, it is recommended to use more points and least-squares adjustment in parameter estimation. By providing residual errors, it is also possible to have a quality control of the process. The least-squares adjustment produces residuals in image coordinates for each point, (r C , r R ), for which the norm is calculated: The program calculates the maximum of these values and the root mean square error, RMSE, in the form: where n is the number of points used. The RMSE and the maximum norm are used to assess the quality of the process. A residual of a few pixels can be accepted, considering the accuracy of GCP coordinates as well as the acceptable accuracy, according to some predefined map scale. The HAPO workflow for HAP orthorectification with the new "Projective 3D" (PRJ3D) program is outlined in Figure 1 and comprises three major steps: (A) acquisition of 2D and 3D coordinates; (B) preparation of ancillary files; and (C) orthorectification processing. In step A, modern orthorectified imagery (MOI) is used to help with orienting raw HAP with our HAPO workflow, namely by establishing the terrain coordinates (X,Y) for GCP's acquisition (Equation (1)) and, later on, to check the planimetric accuracy of HAP orthomosaic. In a first approach, aerial photos are approximately georeferenced in ArcGIS, using a projective transformation, without any orthorectification, which is generally deemed acceptable in flat areas. In areas with significant relief, this procedure yields poor accuracy, often with discrepancies of 100 m, when compared to other precise geographic information. Nonetheless, this approximation facilitates the identification of common ground points when overlaid on modern orthoimages. After loading each HAP into ArcGis and fitted to display with MOI, at least 10 common points (GCPs) between HAP and MOI should be identified. The 2D coordinates (column C, row R) are provided by historical imagery, and 3D coordinates are given by MOI (easting X, northing Y) and a digital elevation model (DEM; height Z) (Output 1, Output 2, and Output 3). Image and terrain coordinates information (C, R, X, Y, Z) is then rearranged into the Coordinates.txt file, which corresponds to Output 4 of Step A.
From this point on (Steps B and C), the PRJ3D program and the obtained output data are run by the open source geospatial software OSGeo4W Shell (http://osgeo4w.osgeo.org) and the Geospatial Data Abstraction Library (GDAL) (http://www.gdal.org/) [44]. GDAL is a computer software library for reading and writing raster and vector geospatial data formats, which is released under the permissive X/MIT style free software license by the Open Source Geospatial Foundation. As a library, it presents a single abstract data model to the calling application for all supported formats. It may also be built with a variety of useful command line interface utilities for data translation and processing [44].   yellow arrows indicate outputs generated by tasks; and white arrows indicate data input required for each task.
PRJ3D also reads a DEM of the area and projects points from ground coordinates, with elevation obtained from the DEM, onto the image space, using Equation (1). Using GDAL command lines as shown in Figure 2, the DEM and each HAP are prepared to comply with the formats required by our PRJ3D software (generating Output 5 and Output 6 of Step B), and Output 4 is run by PRJ3D in order to create residual files (Output 7 and Output 8 of Step B) for quality control. In this way, all the necessary ancillary files are obtained to proceed to Step C. At this point, PRJ3D is run and Outputs 4-8 with GDAL command lines to perform HAP orthorectification. For the obtained (C,R) image coordinates, it interpolates a greyscale value from the input image, using bilinear interpolation. In this way, an orthoimage is produced, with pixel data being brought to their correct geographical position. The orthoimage is created with a resolution specified by the user, normally larger than the DEM resolution. Nevertheless, the ideal resolution for the orthoimage is the mean ground resolution of the original photo. Finally, the obtained orthorectified HAP is converted into a tiff raster file.

Test Area and Remote Sensing Data
The test case was conducted in a mountainous region with a Mediterranean-type climate in Northeast Portugal, covering an area of 118,636 ha and ranging in elevation from 92 m to 1199 m a.s.l. (Figure 2). The area is included in the European Union's Natura 2000 network of nature conservation sites [45], and in the Baixo Sabor Long Term Ecological Research Site (LTER_EU_ PT_002). This area provides a valuable test case for our methodology, because it is a rugged landscape where an accurate orthorectification process is crucial to obtain detailed information on historical land cover. Furthermore, the landscape studied has suffered major land-cover changes over the years, mainly due to large scale farmland abandonment since the 1960s, thus presenting additional challenges [2,[16][17][18][19]25,26], especially regarding GCPs acquisition. Finally, the availability of detailed information on land-cover dynamics was required in the scope of long-term socioecological studies carried out at this LTER site, aiming to describe and understand long-term trends in biodiversity and ecosystem services [46], which could then be used in adaptive land management strategies [8,47,48]. For our test case, we used 74 HAP of the Portuguese Army Geographic Institute (IGeoE) 1947 RAF collection covering the study area. This HAP collection derives from a systematic survey of Portugal made by the British Royal Air Force (RAF) in 1947 and is one of the largest and oldest collections of vertical aerial photos available in the archives of IGeoE in Lisbon [15]. Despite their relevance [15], they have rarely been used in ecological studies and remain relatively unexplored overall. In fact, the set of photographs used in this study had never been acquired from IGeoE in more than 65 years, nor have they ever been orthorectified.
The photos were digitalized by IGeoE in a photogrammetric scanner with a resolution of 0.021 mm [18] and presented a scale of ca. 1:29,500. There was neither a calibration certificate nor any information about the used camera, calibrated photo coordinates for fiducial marks and principal

Test Area and Remote Sensing Data
The test case was conducted in a mountainous region with a Mediterranean-type climate in Northeast Portugal, covering an area of 118,636 ha and ranging in elevation from 92 m to 1199 m a.s.l. (Figure 2). The area is included in the European Union's Natura 2000 network of nature conservation sites [45], and in the Baixo Sabor Long Term Ecological Research Site (LTER_EU_ PT_002). This area provides a valuable test case for our methodology, because it is a rugged landscape where an accurate orthorectification process is crucial to obtain detailed information on historical land cover. Furthermore, the landscape studied has suffered major land-cover changes over the years, mainly due to large scale farmland abandonment since the 1960s, thus presenting additional challenges [2,[16][17][18][19]25,26], especially regarding GCPs acquisition. Finally, the availability of detailed information on land-cover dynamics was required in the scope of long-term socioecological studies carried out at this LTER site, aiming to describe and understand long-term trends in biodiversity and ecosystem services [46], which could then be used in adaptive land management strategies [8,47,48].
For our test case, we used 74 HAP of the Portuguese Army Geographic Institute (IGeoE) 1947 RAF collection covering the study area. This HAP collection derives from a systematic survey of Portugal made by the British Royal Air Force (RAF) in 1947 and is one of the largest and oldest collections of vertical aerial photos available in the archives of IGeoE in Lisbon [15]. Despite their relevance [15], they have rarely been used in ecological studies and remain relatively unexplored overall. In fact, the set of photographs used in this study had never been acquired from IGeoE in more than 65 years, nor have they ever been orthorectified.
The photos were digitalized by IGeoE in a photogrammetric scanner with a resolution of 0.021 mm [18] and presented a scale of ca. 1:29,500. There was neither a calibration certificate nor any information about the used camera, calibrated photo coordinates for fiducial marks and principal point, nor any lens distortion information. Other several classical HAP issues [26] were identified, namely: (i) a poor radiometric quality with non-uniform luminosity or image darkness; (ii) high radiometric variability within images from the same flight; (iii) a nonconstant overlap of the images; and (iv) the lack of a flight map or historical ground control points (GCPs). As the images were not stereo pairs, they presented a low and nonconstant spatial overlap and therefore were treated separately.
After loading each 1947 HAP into ArcGIS and fitting it to display with 2007 modern orthoimages provided by Portuguese Forestry Resources Institute (DGRF), we identified 10 common points per photo between 1947 images and 2007 orthophotos. To help to establish a reliable and high-quality collection of historical GCPs given the considerable land-use changes that occurred over time, we used complementary historical cartographic products at the scale of 1:500,000, including Portuguese Roads Map (1951), Diagram of the Portuguese Railways (1929) and Itinerary Map of Portugal (1937). Such historical road and railway maps often have rather unique and regionally specific characteristics, which can be useful for georeferencing HAP in regions such as Europe [3]. For each photograph, GCPs were acquired in landscape features that remained relatively unchanged and recognizable over time, such as tree fences, stone walls or roads.
The data to adjust Equation (1) included the 2D coordinates (C,R) of the GCPs identified in the 1947 historical imagery, and 3D coordinates (X,Y,Z) for the same GCPs collected from 2007 orthophotos and a digital elevation model (DEM). The DEM was generated from modern altimetry data (contours and spot heights) at a scale of 1:25,000, produced by IGeoE, with a 6 m × 6 m resolution. We then applied the HAPO workflow as described in Section 3.1. For each photograph, the root mean square error (RMSE) was calculated, as well as the maximum residual value. Their values, initially in pixel units, were converted to meters. An accuracy of 5 m (0.2 mm in a 1:25,000 reference scale) was taken as an acceptable RMSE.
In the final step, images were reformatted for mosaic building, discarding the peripheral areas, more prone to distortion problems, and focusing on the central portion of the photograph (i.e., the portion with minimum relief displacement). Since the small overlaps between adjacent photos were not regular, the area to discard was not the same for all the photos. Therefore, each photo was treated separately to obtain its maximum usable surface.

Photointerpretation and Land-Cover Mapping
After successfully being orthorectified and mosaicked, historical aerial images were subjected to visual photointerpretation to produce land-cover maps for test areas. Specifically, we produced a regular tessellation of the study area with 2 km 2 hexagon units, which was considered a sampling unit of sufficient size to capture spatial variations in landscape composition and heterogeneity [49,50]. A total of 56 test areas (≈10% of the study area) were then selected through a stratified random sampling approach [51], with strata corresponding to municipalities.
The procedure was based on visual interpretation, because although time-consuming, it provides reliable results [19], with trained aerial interpreters being able to successfully identify a range of different land-cover types [2,3,8,52]. Automatic or semiautomatic techniques were not used, because the radiometric quality of our historical images was not satisfactory for such approaches [19,21]. Furthermore, when a single panchromatic image is available (such as a gray-scale image obtained by HAP digitizing), an automatic image classification, in a similar way to what is done with multispectral images, is simply not possible [53].
Visual photointerpretation was made with the support of ancillary geospatial data, such as DEM, topographic maps, and other orthophotos of the same area [5,54]. We also used complementary historical data [3], such as Census data [55,56], as well as official agrarian statistical and technical reports [57][58][59][60], to improve the interpretation of landscape features. These records contained demographic as well as natural resource and agricultural information, which helped to evaluate the dominant land cover/land uses in the landscape at the time. Additionally, a site survey was carried out to check small arboreal patches that could be identified in 1947 photos. The rules of photointerpretation followed as closely as possible those used on the latest official land-cover map of Portugal (COS2007) [61], to make the data comparable across time periods.
Land-cover patches and classes were interpreted at a 1:4000 scale and hand-digitized using ArcGIS 10.3 [62]. The minimum mapping unit (MMU) was 1 hectare according to the criteria used in COS2007. Using Conefor Sensinode 2.2 [63], polygons (of one same class) under 1 hectare within 60 m from each other were grouped into multipolygon features until they complied to the MMU, assuring that (i) fine-scale information was successfully acquired and maintained in the analysis and (ii) comparability of multitemporal data sets was ensured.
The accuracy of the land-cover photointerpretation using 1947 photos could not be established quantitatively, because there were no historic reference data, and subsequent land-cover changes rendered field validation impossible [2]. Although field visits allowed the detection of some arboreal patches identified in 1947 photos, others had disappeared or changed considerably and so could not be used to validate historical land-cover mapping. To minimize these problems and thus enhance mapping accuracy, particular care was taken to cross-check photointerpretation with other sources of information, such as census data [55,56] and the official agrarian statistical and technical reports [57][58][59][60].

Historical Orthophotographs and Orthomosaic
Application of the HAPO methodology successfully processed all images (Figure 3), producing an orthomosaic of 74 orthorectified HAP for our test area (Figure 4). Visual inspection of the orthorectified images overlaid on modern 2007 orthomaps and road networks indicated a good matching between old and recent geographic information ( Figure 5). This was confirmed by the low average value and standard deviation of RMSE (1.97 m and 1.00 m, respectively), which ranged between aerial photographs from 0.90 m to 5.88 m. Only one out of 74 images processed had an RMSE above 5 m (Figure 6), which was taken as the reference value for acceptable accuracy. The maximum errors were slightly higher, still with only 2 values above 10 m. Overall, images covering areas with steeper slopes and with more pronounced landscape changes over the years had the highest RMSE values. Nevertheless, the accuracy obtained was well within that acceptable for the reference map scale of 1:25,000.

Land-Cover Mapping
Our historical aerial photographs (Figure 6a) allowed the identification of eight different land-cover classes in the landscapes of 1947: annual crops (covering 61.24% of the total photointerpreted area), perennial crops (14.21%), shrublands (11.77%), natural forests (7.88%), seminatural forests (1.33%), planted forests (3%), water (0.28%), and urban (0.28%) (Figure 6b). It was possible to confidently identify and delineate small to medium-sized and large patches, ranging from 0.1 to 200 ha, allowing a detailed socioecological description of the 1947 landscape in our test area, despite the marked landscape change that occurred over time.

Discussion
Our study described and provided a successful application of the new HAPO workflow to orthorectify historical aerial photographs (HAP). The local case study developed in Portugal contained many of the problems commonly faced by researchers wanting to use HAP and showed that HAPO can produce spatial information with the accuracy required to describe and understand long-term socioecological patterns and processes. We showed that HAPO can deal with non-

Discussion
Our study described and provided a successful application of the new HAPO workflow to orthorectify historical aerial photographs (HAP). The local case study developed in Portugal contained many of the problems commonly faced by researchers wanting to use HAP and showed that HAPO can produce spatial information with the accuracy required to describe and understand long-term socioecological patterns and processes. We showed that HAPO can deal with non-overlapping or marginally overlapping HAP sets, lack of a flight map, difficulties to identify ground control points (GCPs) in changing landscapes, poor radiometric quality, and high radiometric variability between images, among other factors. Despite these problems, the final orthoimages had high geometric accuracy and could be used to produce land-cover maps that will later feed into studies addressing questions related to long term socioecological dynamics. Therefore, we believe that HAPO can be used by a wide range of researchers dealing with similar problems elsewhere, representing a significant advance to harness the power of historical information available in HAP.
Through the application of HAPO to our particular test case, it was clear that achieving high accuracy in HAP orthorectification and georeferencing requires that care should be taken with a number of critical aspects. First, a good resolution digital elevation model is required to achieve the best results, particularly in rugged landscapes such as ours, where significant changes in elevation over small areas represent a challenge for the correct application of our procedure. Second, it is important to have access to modern orthorectified and georeferenced imagery, or other adequate cartographic information, without which it may be nearly impossible to obtain accurate GCPs. Third, GCPs should be well distributed along the historical image and preferentially not placed near the borders to avoid highest correspondence errors between HAP and modern orthoimages. Fourth, when covering highly rugged areas, steep slopes and landscapes suffering pronounced changes over the years, it is recommended to acquire a larger number of GCPs to help to reduce positional errors. Finally, we highlight the importance of using complementary cartographic products such as historical road and railway maps, as well as historical datasets such as Census data and unpublished reports contemporary to the HAP, which can be very helpful to identify GCPs in the orthorectification process and to support the visual photointerpretation to produce land-cover maps.
Overall, the new HAPO workflow presented and tested in this study represents a valuable addition to the toolkit of researchers interested in landscape ecology, socioecological patterns and dynamics, long-term ecology, and global change, by offering the possibility to capture the unique and extremely valuable historic memory of our landscapes, and permitting a detailed and objective look-back in time [21]. Reconstructing and recovering historical patterns and processes is key to understanding how present conditions came about, with the analysis of landscape change drivers providing an essential basis [64] not only for monitoring landscape diversity [57,58] and for investigating changes in vegetation and landscape structure [35][36][37], but also to understand human-environment interactions and how they have shaped contemporary landscapes [64,65]. This unique information can be used by managers and researchers for a great variety of environmental studies and applications [1,2] and can definitely contribute to more informed policy, management, and restoration decisions [48,[64][65][66][67][68][69].
Future work may consider the adaptation of the program to the Python programming language in order to include it as an additional GIS command, for example, in the form of a QGIS plugin. This will make it more friendly for nonspecialist GIS users.
Author Contributions: A.T.P. conceived and designed the work, collected data, and wrote the paper; J.A.G. developed the PR3J program; A.T.P., P.B., and J.H. analyzed the data; all authors supported paper writing and edits.

Conflicts of Interest:
The authors declare no conflict of interest.