Using aerial photography to estimate wood suitable for charcoal in managed oak forests

Mexican oak forests (genus Quercus) are frequently used for traditional charcoal production. Appropriate management programs are needed to ensure their long-term use, while conserving the biodiversity and ecosystem services, and associated benefits. A key variable needed to design these programs is the spatial distribution of standing woody biomass. A state-of-the-art methodology using small format aerial photographs was developed to estimate the total aboveground biomass (AGB) and aboveground woody biomass suitable for charcoal making (WSC) in intensively managed oak forests. We used tree crown area (CAap) measurements from very high-resolution (30 cm) orthorectified small format digital aerial photographs as the predictive variable. The CAap accuracy was validated using field measurements of the crown area (CAf). Allometric relationships between: (a) CAap versus AGB, and (b) CAap versus WSC had a high significance level (R 2 > 0.91, p < 0.0001). This approach shows that it is possible to obtain sound biomass estimates as a function of the crown area derived from digital small format aerial photographs.


Introduction
Across the global south, charcoal is frequently produced from clear cutting or selective logging of woody plants that resprout (Chidumayo and Gumbo 2013). The sustainable management of resprouters for charcoal production requires data on their structure and growth capacity; which in turn are affected by the species' sprouting ability (strong versus weak) (Bond and Midgley 2001), forest management practices (clear cutting versus thinning; short versus long rotations) (Aguilar et al 2012), and local biophysical conditions (Leonardsson and Gotmark 2015). In other words, good quality data on the resprouters' structure and growth are site-specific and should be generated in situ if adequate sustainable forest management plans are to be designed and implemented.
In order to design sustainable charcoal production programs in data-poor landscapes, cost-effective methods are needed for measuring the aboveground woody biomass productivity of resprouters. One way of attaining this is to estimate the biomass stocks between multiple time intervals, so as to model regrowth rates. However, while past field measurements of biomass stocks are missing for most charcoal production areas of the world, historical aerial photographs could be available.
The most accurate methods for estimating aboveground woody biomass are those involving direct measurement of individual trees by destructive sampling (Brown et al 1989, Brown 1997, Pacala et al 1996, Návar 2009). Common field-measured predictive variables for the development of allometric models are diameter at breast height (DBH) and total tree height (Picard et al 2012). However, these laborintensive methods alone do not provide information about the continuous spatial distribution of woody biomass (Lu 2007). This can be overcome by integrating field-based allometric equations with remote sensing and spatial analysis (Castillo-Santiago et al 2013, Dubayah and Drake 2000, Dong et al 2003, Lefsky et al 2005, Lu 2007, Zolkos et al 2013, a technique that has proved useful even for regional and global wall-to-wall aboveground biomass (AGB) estimates (Saatchi et al 2011, Cartus et al 2014, Baccini et al 2012, Avitabile et al 2011.
Methods less dependent on field data rely on the availability of high spatial resolution products (< 1 m per pixel) such as high-resolution digital aerial photographs, multi-and hyperspectral sensor imagery as well as light detection and ranging (LiDAR) systems coupled with modern computation and data management capabilities. These allow us to infer tree descriptive parameters such as crown area, biomass volume, tree density, tree height and species composition (Drake et al 2003, Leckie et al 2003, Holmgren and Persson 2004, Koukoulas and Blackburn 2005, Naesset and Gobakken 2008, Gougeon and Leckie 2006, Massada et al 2006, Popescu 2007, Breidenbach et al 2010. However, these tools are often not feasible to implement due to budget constraints. But most important, the newer the sensor, the less the availability of past images.
When dealing with high-resolution biomass estimations (such as those needed for local woodland management plans), allometric models have been calibrated by relating field-based DBH measurements with AGB determination in order to develop empiric equations to estimate biomass (Brown et al 1989, Brown 1997, Baccini et al 2004, Aguilar et al 2012, Segura and Kanninen 2005. The estimation of tree biomass from aerial imagery has been around for at least two decades. Earlier approaches (e.g. Brandtberg andWalter 1998, Korpela 2004) were based on either automated or manual analysis of digitized analog aerial photographs and photogrammetrical software. These techniques gave a good overall agreement between image-based, ground-based crown diameter and total stem volumes, but rather limited accuracy regarding individual tree height in dense forest (Korpela 2004, Naesset 2002. But the level of agreement depends on adequate single-tree crown delineation, which is affected by forest density, structure, species composition as well as tree position within the photograph. Korpela (2004) reports that the best delineation occurred in trees close to the photograph nadir. More recent approaches include digital aerial photographs, multispectral sensors and LiDAR sensors.
None of these approaches deal specifically with young sprouts typically exploited for charcoal production, which have quite a different architecture to non-managed stands (Aguilar et al 2012, Chidumayo 2014. Sousa et al (2015), for example, report a canopy to the biomass model that could be useful for bioenergy purposes, but oaks in their study area were never cleared cut, but '…pruned with a periodicity of approximately 5 years in order to enhance fruit production and the stands are predominantly mature'. Furthermore, none of the previous studies dealt specifically with wood suitable for charcoal (WSC), as the portion of AGB actually useful for charcoal making. Oak forests and rangelands (Quercus spp.) in many developing countries across the tropics represent a major woodfuel source, and are intensively coppiced for charcoal on rotations of varying duration, usually by clear cutting in order to save time and effort carrying trunks and thick branches to the kilns (e.g. Aguilar et al 2012, Jomaa et al 2003, Aus-Der-Beek et al 2006, Herrera and Chave 2006, Nygren 2005, Wakeel et al 2005, Shrestha 2003, Barrow andHicham 2000, Somanathan 1991). Oaklands traditionally managed for charcoal are no exception regarding the lack of legal permits for exploitation and commercialization (Masera et al 2015). It has been reported that many charcoal producing areas across developing countries have been traditionally managed without any kind of permits (Masera et al 2015, Sander et al 2013, Zulu and Richardson 2013. Barring a few exceptions (e.g. Chidumayo 2014, Castillo-Santiago et al 2013, Chidumayo 1991, 2013, Chidumayo 1988, charcoal producing areas have been neglected by foresters and ecologists due to their low economic potential as timber sources, and high level of anthropization, respectively.
Similar to other data-poor landscapes, there is a need for high-resolution estimates of WSC.
The present study was designed and performed to calibrate allometric models relating the crown area of young oak sprouts measured in aerial photographs with AGB and WSC at tree and stand level. The underlying objective was to show that high-resolution estimates of AGB and WSC, aimed as a source of information for sustainable charcoal production, can be produced in a simple and relatively cost-effective way with little previous information on the characteristics of the managed woodland. Therefore, this framework could be implemented where applicable in similar settings.

Methods
Photogrammetric measurements of the tree crown area were obtained from very high-resolution (30 cm) small format digital aerial photographs to derive indirect allometric equations and to estimate the basal area (BA), AGB and WSC. To account for the existing variation, the estimates of BA, AGB and WSC were obtained at two different levels of analysis: (1) individual trees, and (2) stand level (aggregation of individual trees). The stand level will be referred to as 'grouped trees' from now on.

Study area
The study area comprises the eastern portion of the Cuitzeo lake basin located in the Trans  The climate is temperate subhumid with dry winters, and from north to south it exhibits a gradient of increasing precipitation and decreasing temperature. The long-term  reports of the Servicio Meteorologico Nacional (National Mexican Weather Service) indicate mean annual precipitation ranges from 695 mm in Cuitzeo (station 16027) in the north to 1116 mm in Acuitzio del Canje at the southmost portion of the basin (station 16001, SMN 2010a, 2010b). While mean annual temperature ranges from 18.2 • C in Cuitzeo town (station 16027) to 16.9 • C in Acuitzio del Canje (station 16001) (SMN 2010a(SMN , 2010b. This corresponds to a transitional zone between the dry and humid temperate climates (López et al 2007). The Cuitzeo basin is a mosaic of land covers, where only 20% of its total extension is covered by mixed pine-oak forests. The dominant land uses are agriculture (rainfed and irrigated) and grasslands covering almost 50% of the total surface of the basin (Mendoza et al 2011). Approximately one million people inhabit the study area (density = 254 pop km −2 ). An important aspect of the population dynamics within the Cuitzeo basin is its emigration to Morelia, the nearest major urban settlement, which accounts for almost 79% of the urban population (López et al 2007). One of the main economic activities of the rural communities is charcoal production using traditional techniques to supply the urban settlements (Castillo-Santiago et al 2013, Serrano-Medrano et al 2014. Quercus castanea Née and Quercus laeta are the most frequently used species for charcoal production, and both are endemic to Mexico (Aguilar et al 2012). These oak species tend to resprout after cutting and therefore may form a tree-trunk structure called a stool in which several boles correspond to one individual. Interviews with land owners indicate that the oak forests in this region have been harvested intensively to produce charcoal for at least 60 years. Clear cutting is the main harvesting practice in the study area. This is because earth mound kilns are built by piling the wood at a location, and so selective logging will significantly increase the effort of carrying the wood into the kiln, which is always done 'by hand'. Furthermore, land owners and cattle grazers are seldom charcoalers, meaning that there is no real interest from key stakeholders for sound oakland management. Typical rotation cycles vary between 5-7 years to 15 or even 20 years, depending on site quality and agreements among land owners and charcoalers within an informal and unregulated production system (Camou-Guerrero et al 2014).

Field survey
No previous information concerning allometric relationships between the crown area and BA for managed oak trees was available for the study area. Therefore, ground inventories were necessary. Two field campaigns were conducted in extensively managed oak forest sites located throughout the eastern region of the Cuitzeo basin (figure 1). The first campaign was conducted during the spring of 2014 and the second during the spring of 2015. Only Quercus spp. individuals were sampled during the field campaigns, and the dominant species in terms of distribution and relative abundance was Q. castanea. Field measurements included: (1) the DBH, typically measured at a height of 1.30 m above the ground and a trunk circumference >10 cm, and (2) the crown diameter (tree crown spread) as the average length of two perpendicular lines across the crown area.
The field BA was obtained from the DBH measurements for each coppice shoot. For tree individuals with  several coppice shoots, we sum up the BA of each shoot to obtain the total BA for each tree (equations 1 and 2).
Since treetops from sprouting oak trees usually exhibit asymmetric architecture, the projected field crown area (CA ) was estimated approximating the crown shape to an ellipse (equation 3, table 1).
where CD fNS and CD fEW represent the field crown diameter in the north-south and east-west direction, respectively. Twelve sites of 9-14 ha were sampled, each one hosting between 35 and 151 trees. Within the sites, the sampled trees were unevenly distributed in 12-25 natural stands of ca. 0.15 ha. Each individual tree was georeferenced using differential correction for improving the accuracy of GPS positioning.

Aerial photography acquisition and processing
High-resolution near vertical small format digital aerial photographs were acquired from an aerial platform (Cessna 182) during March 2014 and April 2015 using a Nikon AF Nikkor 24 mm objective and small format digital cameras Nikon D70 SLR and D300S (figure 2).
Both cameras were previously calibrated to assess lens distortion and swath distance to flight altitude relation. Central image position coordinates were registered with a WAASS-enabled GPS unit (GPSMAP 60CS, Garmin, Inc.), while roll, pitch and yaw angles were recorded during flight by means of an Arduino UNO board coupled to an Arduino compatible logger shield (Adafruit, Inc.) and a 9 degrees of freedom inertial measurement unit (9DOF Razor IMU, SparkFun Electronics, Inc.). Furthermore, aerial photographs

Crown area estimation using aerial photographs
Individual tree crowns were delineated using a mirror stereoscope and high-contrast color prints of the aerial photographs. The tree crown limits were on-screen digitized using GRASS GIS v6.4.3 (GRASS Development Team 2012). Individual tree measurements were aggregated to obtain grouped tree areas (figure 3). Field verification and delineation of interpreted features occurred within the first month after the acquisition of the image.

Biomass estimation using photo-derived crown measurements
To estimate the total AGB and WSC using the crown area measured in aerial photographs (CA ap ) as the predictor variable instead of DBH, first we evaluated the existence of an allometric relationship between the field-measured DBH (expressed as BA) of the oaks' coppice-shoots and CA ap . The data gathered during the two field surveys were divided into two sets as follows: (1) field measurements of the first field survey (2014) were used to evaluate the BA versus CA ap , and (2) field measurements for the second field campaign were used to estimate the biomass and generate the allometric equations between CA ap versus AGB and CA ap versus WSC.
We calculated the BA from the DBH field measurements for Q. castanea and Q. spp. coppice shoots. Since typical sprouting oak trees or stools managed for charcoal making exhibit several coppice shoots, as opposed to unharvested oaks trees, it was necessary to sum up the BA of the coppice shoots belonging to the same tree (equation 2). In this way it was possible to establish a 1:1 relationship between the total individual tree BA and individual tree photo-measured crown area (CA ap ). We followed the same approach to evaluate the allometric relationship between the BA for grouped trees (Q. spp) and the CA ap corresponding to those grouped trees, i.e. we sum up the total stand BA and total stand CA ap .
Regressions models were generated to evaluate the allometric relationship between BA versus CA ap and we identified a strong linear relationship between BA versus CA ap ( figure 4).
Then, we proceeded to calculate the total AGB and WSC. We used as input values to estimate the biomass, the field measurements from the DBH for single coppice shoots (n = 1604 coppice shoots, 366 individual trees). The allometric equations that allow us to estimate both the total AGB and WSC from the field-measured DBH were generated in a previous study conducted by Aguilar et al (2012) in the same study area (table 2).
The AGB and WSC were obtained for each coppice shoot. For tree individuals with several coppice shoots, we sum up the WSC and AGB of each shoot to obtain the total biomass component for each individual tree.

Statistical analysis
Once the biomass components were obtained, linear regression models were generated for the two levels of analysis (individual and grouped trees) and for both sample periods (2014 and 2015) to evaluate the allometric relationships between: (1) CA f versus CA ap , (2) CA ap versus BA, (3) CA ap versus AGB, and (4) CA ap versus WSC. Field-and photo-measured crown area were used as independent variables to generate the allometric models. All variables were log-transformed for the adjustment of linear regressions. The statistical analysis was conducted using the R v3.2.0 software (R Development Team 2015). Notes: Equations are in the form y = a(DBH) b , 'y' in kg per coppice-shoot and 'DBH' in centimeters per coppice shoot. Parameters a (scaling coefficient) and b (scaling exponent) were solved by the Gauss-Newton non-linear least squares algorithm. All scaling exponents differ significantly from 1 (p < 0.0001).

Figure 4.
Linear models between field-(CA f ) and photo-measured (CA pa ) crown area versus tree BA for individual Q. castanea trees (top) and grouped Q. spp. trees (middle) and linear models between field-measured crown area (CA f ) and photo-measured crown area (CA pa ) for individual Q. castanea trees (top) and grouped Q. spp. trees (bottom).

Results
In both field campaigns, we found significant relationships between field-measured CA and photo-measured CA for both individual trees (R 2 = 0.88, p < 0.001) and grouped trees (R 2 = 0.96, p < 0.001) (figures 4 and 5). Thus, the estimation of the BA was expressed as a function of CA ap in the form Ln(AB) = Ln (AC ap ) + b. For individual trees sampled during the spring of 2014, the R 2 values were equal to 0.73 for individual Q. castanea trees and 0.90 for grouped trees (Q. spp). Allometric relationships were highly significant (p < 0.001) in all cases. The adjusted R 2 values were higher when the crown areas of the individual trees were aggregated. As we increased the sample efforts during spring 2015, the relationship between CA ap and BA for the individual and grouped trees became more accurate than that observed during 2014 (figure 5).
The regression analysis for the AGB and WSC estimation using CA ap as the predictor variable was highly significant with R 2 values of 0.82 for Q. castanea trees and 0.91 for grouped trees (figure 6). These results indicate that it is possible to use CA ap as a proxy variable for the estimation of the total AGB and WSC (table  3). Similar to the results obtained for the relationship between CA ap and AB, biomass estimations for individual Q. castanea trees showed adjusted R 2 values lower than the grouped trees.

Discussion
Although the idea of estimating biomass from the crown area measured in aerial photographs is not novel (Ilvessalo 1950), this study explored the generation of a cost-effective framework for estimating biomass by using an intensively managed oak forest as a case study. No previous approaches of this kind have been applied before to these ecosystems and such information is urgently needed due to its implications for forest management and woodfuel security. The only accurate equations available to estimate the total AGB and WSC in our study area were developed by Aguilar et al (2012). However, they used DBH as the predictive variable. In this study, we present an approach to generate sound biomass estimates using the tree crown area measured in high-resolution aerial photographs, instead of field-measured DBH, as the predictor variable. Although tree height is considered an important variable in forest inventories (Korpela and Anttila, 2004) as well as in AGB estimation as clearly stated by Chave et al (2014), when aerial photographs at detailed scales (1:10 000−1:15 000) are digitized with very high-resolution scanners (∼20 m) yielding a submeter nominal pixel resolution and processed with photogrammetrical techniques, mean tree height can be estimated. Nevertheless, with this approach, in forests stands, height is normally underestimated (Naesset 2002, Korpela andAnttila 2004) with errors ranging from 2.47−5.42 m (Naesset 2002, Korpela andAnttila 2004). Correction of these errors implies field measurement of tree height at sampling points to better predict mean tree height at the stand scale using a regression between aerial image-derived mean stand height and field-measured mean stand height (Naesset 2002). With this approach, the predicted mean heights overestimated true heights with less precision in very young forest stands, and higher precision in young and mature stands (Naesset 2002). Due to the associated uncertainty when tree height is to be estimated using photogrammetrical methods coupled with small-frame digital images, we decided to focus on the CA ap versus BA relation and employ the AGB and WSC models that did not need tree height as an input. Biomass estimates at a regional scale are useful for generating information across large areas (Jenkins et al 2003, Baccini et al 2004, Ralevic et al 2010, Castillo-Santiago et al 2013. However, at a local scale, field measurements still constitute the most accurate method to ensure that the allometric equations reflect the actual growing stock volume (Nogueira et al 2008, Návar 2009, Picard et al 2012. For instance, trees from intensively managed (i.e. coppiced) oak forests exhibit different architectural and growth patterns to those of primary oak forests, mainly due to the reduced growing space available to coppice shoots, which develop typical stool structure (Ducrey and Turrel 1992, Bellingham and Sparrow 2000, Logli and Joffre 2001, Espelta et al 2003, Konstantinidis et al 2005, Sands and Abrams 2009. In contrast, mature open-growth trees present very large crowns in proportion to their DBH and species (Johnson 2000). As a result, the use of standard allometric equations developed for primary oak forests, which exhibit different tree proportionality ratios, can alter biomass parameters leading to gross errors in biomass estimation (Picard et al 2012).
The equations developed in this contribution are relatively simple. However, there are several aspects that need to be addressed before applying them in other regions. First, the equations were developed for coppice oak trees with DBH ≤ 30 cm, so attention should be paid if they are used to estimate the biomass of non-resprouting oak trees with larger DBH values. A second aspect is related to image acquisition date and time. In this regard, varying phenological and lighting conditions can affect the spectral responses of the Quercus species and the discernibility of trees, leading to biased crown area estimates Leckie 2006, Barbier et al 2010). For Quercus tree species, particularly Q. castanea, the spectral response of leaves varies seasonally, facilitating the identification of oak trees in aerial photographs acquired before spring (from mid-February until mid-March). Furthermore, it is important to consider that parameters of individual trees measured in aerial photographs are indirect, and the allometric equations are inherently uncertain and susceptible to systematic errors (Korpela 2004). Thus, very high-estimation accuracy at individual tree level cannot be expected (Gougeon and Leckie 2006, Vastaranta et al 2011, Kaartinen et al 2012. As illustrated in this study, the allometric models developed to test the relationships between AC ap versus AGB and AC ap versus WSC showed in all cases a higher R 2 for grouped trees. Nonetheless, the analysis of this work at the individual tree and stand level (grouped trees) provides useful tools for obtaining accurate measurements of the crown area and precise biomass estimations. This makes the proposed framework less site-specific, allowing its replicability in other regions, subject to similar management practices.
Further research must evaluate the applicability of the methodology presented in this paper and the possibility to process historical aerial photographs in order to estimate the AGB and WSC at different dates. This would make it possible to analyze the AGB and WSC changes in time and evaluate the AGB productivity rates at different locations.
Finally, one aspect that should be considered to improve this framework applicability and make it fully operational is related to the degree of treecrown identification and delineation automation. Multi-scalar automated or semi-automated algorithms for identifying individual trees within forest stands still need to be systematically tested and validated. Although very promising software and approaches to perform automated tree-crown identification and delineation exist, e.g. eCognition by Trimble Inc. (www.ecognition.com/) or SPRING by Brazil's National Institute for Space Research (www.springgis.org/languages/english/index.html), the research on this subject highlights the need for a sound baseline for calibrating and validating automatic tree-crown identification and delineation (Gao et al 2006, Gao andMas 2008). In this respect, the information generated by visual interpretation (i.e. manual tree-crown identification and delineation) is crucial for calibration and validation purposes, which provide feedback to further improve the results of the automated techniques.

Conclusions
The aim of this work was to develop a framework for the estimation of the AGB and WSC at the individual tree and stand level using information from high-resolution aerial photographs coupled with field measurements and species-specific allometric equations previously developed for the study site. We used as a case study an intensively managed oak forest used to produce charcoal. We identified strong and highly significant relationships between the photo-measured crown area and BA, AGB and WSC. The advantages of this framework are: (1) the reduction of field measurements and the associated costs, (2) it can be easily replicated, since it requires a minimum input of parameters (i.e. crown area measurements), (3) it is flexible, since it works at different levels of analysis (i.e. individual trees and grouped trees), (4) it allows us to predict values for the BA, total AGB and WSC in terms of kilograms of dry matter and volume and, (5) it allows us to combine different sources of information to support forest management programs focusing on sustainable energy development. Future research is needed on costbenefit comparisons between aerial photographs and high-resolution satellite imagery for estimating WSC in woody landscapes. This information could eventually be compiled into a good guidance practice of 'best' methods according to a variety of environmental and socioeconomic conditions.