Strategic expansion of existing forest monitoring plots: a case study using a stratified GIS-based modelling approach Thai Son Le1,* and Justin Morgenroth2 1Vietnam National University of Forestry, Hanoi, Vietnam 2New Zealand School of Forestry, University of Canterbury, Christchurch 8140, New Zealand *Corresponding author: Thaisonfuv@gmail.com (Received for publication 10 February 2019; accepted in revised form 14 October 2020) Abstract Background: Understanding the relationship between sites and the plant species they support is essential for effective vegetation management. Site-species matching requires knowledge of the growth response of a given species to the full range of environmental conditions in potential planting sites. This can be achieved by repeatedly measuring species growth at a comprehensive network of sample plots that cover a range of environmental conditions, including topography, climate, and soil factors. The New Zealand Dryland Forests Initiative has established permanent sample plots (PSPs) of a plantation species, Eucalyptus bosistoana F.Muell., across New Zealand. However, these PSPs do not cover the entire range of environmental conditions available for the species and hence there is a need to expand the network of sites. The aim of this study was to determine optimal locations for new PSPs to provide more unique information to support site-species matching studies for Eucalyptus bosistoana in New Zealand. Methods: A geographic information system (GIS) and stratified random sampling method were used to generate a model to identify optimal locations for E. bosistoana PSP establishment. The variables used in this study included topography, climate, and soil data. Redundancy between the initial set of potential explanatory variables was reduced by a multi-collinearity analysis. The potential habitat for the species was restricted to land with environmental conditions that could support E. bosistoana. All environmental variables were stratified and an initial priority index for each stratum in each variable was calculated. Then a weighted-overlay analysis was conducted to create the final priority index, which was mapped to identify high-priority areas for targeted PSP expansion. Results: The existing PSP network for E. bosistoana generally covers the environmental conditions in low-elevation New Zealand dry lands, which are located alongside the east coast of the South Island, and the southern part of the North Island. The model identified high priority areas for PSP expansion, including several large regions in the North Island, especially in Rangitikei and Taupo Districts. Conclusions: The model successfully allowed identification of areas for a strategic expansion of permanent sample plots for E. bosistoana. Newly identified areas expand upon the topographic, climatic, and soil conditions represented by the existing PSP network. The new area for PSP expansion has potential to provide valuable information for further site-species matching studies. The methodology in this paper has potential to be used for other plot networks of a different species, or even natural forests. New Zealand Journal of Forestry Science Le and Morgonroth New Zealand Journal of Forestry Science (2020) 50:9 https://doi.org/10.33494/nzjfs502020x41x E-ISSN: 1179-5395 published on-line: 05/11/2020 © The Author(s). 2020 Open Access This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. Research Article Open Access Keywords: Eucalyptus bosistoana; GIS; habitat modelling; NZDFI; permanent sample plot network; site-species matching; stratified random sampling. http://creativecommons.org/licenses/by/4.0/), Le and Morgonroth New Zealand Journal of Forestry Science (2020) 50:9 Page 2 Introduction The use of forest monitoring plots is important for forest resource protection and management. For plantation forest development, permanent sample plots (PSP) provide information on how trees grow, when and which silvicultural practices are required, and how effective they are if applied. This information is crucial for the decision making process of any plantation management plan and can be used to support decisions about whether more plantings of a given species should be established (Millen et al. 2016). Understanding the relationship between a species and its site requirements is essential to the success of the plantation. The New Zealand Dryland Forest Initiative (NZDFI) has used permanent sample plots to monitor the survival and growth of durable eucalypts not previously planted in large numbers in New Zealand. One particular species being trialed by the NZDFI is Eucalyptus bosistoana F.Muell., which has its natural habitat in the south- eastern coastal areas of Australia (Boland et al. 2006). This species has good potential as a plantation species in New Zealand due to its abilities to produce highly durable timber, to coppice vigorously after fire and harvesting, and to provide nectar/pollen for native biodiversity (Apiolaza et al. 2011; Millen et al. 2016; Nicholas & Millen 2012). Despite many positive characteristics, there are some concerns about its susceptibility to pests, namely defoliation from the Eucalyptus variegated beetle (Paropsisterna variicollis (Chapuis)) (Lin et al. 2017). With its hardness, wood density, straightness, and durability, the timber of this species has been used for numerous purposes such as farming fences, building materials, boat masts and railway ties. Of particular relevance is that this species has high drought tolerance that satisfies the prerequisite for plantations in New Zealand drylands (Apiolaza et al. 2011; Millen et al. 2016; Nicholas & Millen 2012). To determine whether this species can be commercially successful in New Zealand, the NZDFI aims to ensure its trials are well distributed throughout New Zealand’s drylands such that they can be used for growth and yield modelling. The growth and yield models will help determine the sites with the greatest growth rates and survival. As of 2015 the NZDFI had established and monitored 84 permanent sample plots (PSPs) of E. bosistoana in 30 sites across the North and South Islands of New Zealand (NZDFI 2015). However, there is a desire to expand the PSP network to span a more complete range of environmental conditions suitable for eucalypt plantations in New Zealand. The extended PSP network will provide a greater understanding of how E. bosistoana grows in a variety of environmental situations, especially in marginal conditions (NZDFI 2015, 2017). To fill the gaps in the existing PSP network, new PSPs should be established at sites with environmental conditions that are not well represented in the current network. Habitat modelling is often used to examine the set of conditions within the habitat of a particular species or a group of species. In this approach, a habitat is considered as a group of separate factors (i.e. terrain, climate and soil). This approach was introduced by James and Shugart (1970) for bird habitat modelling. In the early days of habitat modelling, the approach was often used to model the habitat or distribution of birds (e.g. Johnston and Temple (1986) or mammals (e.g. Pereira (1989); Pausas et al. (1995). These habitat models were typically created by collecting site description data in terms of site properties (i.e. model variables) to determine their statistical relationships with the distribution of a species. Such statistical assessments usually require large amounts of empirical data that are not easily collected. As computing technology became more widely available, habitat modelling took advantage of these processing advances. A computer-based model, CLIMEX, was developed and used to describe the climatic favourability of a given location for a particular animal species (Sutherst & Maywald 1985). CLIMEX modelling requires climatic data and population data for the target species (Sutherst & Maywald 1985; Sutherst et al. 2007). CLIMEX and its applications have enabled ecologists to utilise computing power for habitat modelling (Pattison & Mack 2008; Shabani et al. 2012; Taylor et al. 2012). Recently, CLIMEX has been applied to examining the potential distribution of invasive plants and insects (Aljaryian et al. 2016; Hill et al. 2016; Xuezhen et al. 2018). In New Zealand, a form of habitat modelling has been widely used to assess the potential productivity of plantation forest species across a range of environmental conditions. Two common indices used to predict the productivity of Pinus radiata include the 300 Index (defined as the stem volume mean annual increment at age 30 years (Kimberley et al. 2005)) and Site Index (defined as the mean top height at age 20 years (Goulding 2005)). Watt et al. (2009) modelled the spatial distribution of Cupressus lusitanica productivity. These indices reflect the potential growth of the target species in their habitat based on surrounding environmental conditions. Environmental site descriptions are costly and time-consuming to undertake, especially in the case of large habitat areas. However, with the development of computer science and geographic information systems (GIS), GIS-based habitat models have become more widespread relative to traditional statistical habitat models because of their applicability on large spatial extents (Basir 2014; Reisinger & Kennedy 1990; Store & Kangas 2001). In GIS-based habitat modelling, each environmental variable is represented by a GIS layer (i.e. either vector or raster). As a result, an environmental description of every spatial location in the habitat can be obtained (Dettmers & Bart 1999; Shaw & Atkinson 1988; Wadge et al. 1993). GIS-based habitat modelling has good potential to meet the NZDFI’s objective of strategically expanding their PSP network for E. bosistoana. A challenge involved with this objective is how to allocate new sample plots optimally across a large study area to cover the wide range of environmental variability. Simple random sampling tends to collect samples throughout the range of values in the population if those values appear with similar frequencies (Green 1979; Royall 1970). However, this method often involves a high risk of bias when applied on naturally distributed population, such as environmental objects on a large geographic surface (Cawsey et al. 2002; Danz et al. 2003; Kohl et al. 2006). Systematic sampling may solve the problem of the simple random sampling method but it cannot guarantee that samples from all important value gradients are taken into account unless the sampling density (i.e. the number of samples over a unit of area) is sufficiently high (J. E. Austin et al. 2001; Scott 1998). Due to these limitations, neither random, nor systematic sampling are appropriate for expanding the E. bosistoana PSP network. Alternatively, to ensure the samples are distributed representatively across the entire range of values, the stratified random sampling method has been widely used, especially for low-density samples over an large areas (e.g. national survey) or an isolated natural area that restricts the ability to collect samples (Esfahani & Dougherty 2014; Tomppo et al. 2014; Yves & Ecker 2014). Indeed, if the environmental conditions within a specific area are stratified into several separate groups, taking samples from these groups will provide more representative information about the population (M. P. Austin & Heyligers 1991). Stratified random sampling has been increasingly adopted in ecological and forestry studies (Danz et al. 2003; Knollova et al. 2005; Wallenius et al. 2011; Yves & Ecker 2014). This study takes advantage of the merits of GIS-based habitat modelling and stratified random sampling to find priority locations for a strategic expansion of the existing permanent sample plot network for E. bosistoana in New Zealand. The resulting priority map highlights the environmental gaps in the existing plot network that need to be filled by the strategic expansion. Although this study is specific to E. bosistoana, the methodology has general applicability as it involves a GIS-based model of habitat that is easy to modify and adjust for other species. Methods Study area The study area includes the North and South Islands of New Zealand, which cover a total of approximately 268,000 km2 (Fig. 1). New Zealand has a wide range of topographic and climatic conditions, ranging from sea level to 3737 m above sea level (Barringer et al. 2002). The annual average temperature ranges from -2.55 to 16.79 °C and annual precipitation ranges from 392 to 6807 mm year-1 (Fick & Hijmans 2017). Tree species and Permanent Sample Plots “Coast grey box” or “Gippsland grey box” (Eucalyptus bosistoana) is a species naturally found on the southeast coast of Australia (Boland et al. 2006). It can reach a height of approximately 40 to 60 m at maturity. The habitat of this species is coastal mixed forests located in areas below 500 m of elevation within a range of latitudes from 33 to 37.5°S (Apiolaza et al. 2011; Boland et al. 2006). The climatic range is warm humid to cool, with a mean maximum monthly temperature of 24-29 °C in the hot season and a mean minimum temperature of 1-6 °C in the winter months. It can withstand 5-40 frost days per year. The mean average precipitation in its natural range is 700-1200 mm a year (Boland et al. 2006), though it is currently planted on sites in NZ with less than 700 mm of annual rainfall (NZDFI 2019). E. bosistoana has been planted at 30 sites across New Zealand to support the NZDFI research programme. Within these sites, 84 PSPs have been established, and together, these PSPs are considered as a PSP network. The PSPs are further divided into 1095 sub-plots (Fig. 1). Description of data used in the study The environmental factors available for use in this study include climate, soil, and topography (Table 1). Climate data are from WorldClim and include average temperatures and precipitation based on the period from 1970-2000 (Fick & Hijmans 2017). Soil data are from the New Zealand Land Resource Inventory (Newsome et al. 2008), while the digital elevation model (DEM) was sourced from Landcare Research who interpolated the surface from Land Information New Zealand’s 1:50,000 topographic spot heights and contours. Le and Morgonroth New Zealand Journal of Forestry Science (2020) 50:9 Page 3 FIGURE 1: The study area (source: Statistics NZ 2016) and existing permanent sample plots for E. bosistoana. General modelling approach The method applied in this study consisted of five steps, each of which is detailed below: (1) identification of variables for GIS-based modelling, (2) building the dataset, (3) variable restriction, (4) variable stratification, and (5) cartographic modelling. Identification of variables for GIS-modelling According to studies by Apiolaza et al. (2011) and Prober et al. (2016) and based on the available sources of data, 17 variables were selected as potentially influencing the growth and distribution of E. bosistoana (Table 2). The multi-collinearity analysis, which involved an assessment of variation inflation factors (VIF), was undertaken to test the correlations between the variables to minimise information redundancy. VIF values near 1 indicate that the variables were independent, while VIFs exceeding 10 were indicative of multi-collinearity requiring correction (García et al. 2014; Kutner et al. 2003). The next step was to determine whether a weighting should be applied to the remaining variables to highlight certain variables as having greater impact upon the growth and development of E. bosistoana. Information from previous studies (Apiolaza et al. 2011; Boland et al. 2006; NZDFI 2015) and expert advice (EG Mason pers. comm) contributed to the decision-making process. Four Le and Morgonroth New Zealand Journal of Forestry Science (2020) 50:9 Page 4 variables were deemed as the most influential factors on the target species, including annual average temperature, precipitation, soil pH and elevation. These variables were given a weight coefficient of 1.5 whereas all other variables were assigned a weight coefficient of 1. Building the dataset In this study, data were derived from existing published data or through interpolation of field data in cases of large areas (Table 2). All acquired data were processed into raster layers (one layer corresponding to a single model variable). All the environmental variables (Table 2) were processed to raster-format layers with cell size of 25 × 25 m in the same projected coordinate system (i.e. New Zealand Transverse Mercator) before being clipped to the study area boundary. The data processing stage was conducted using ArcGIS version 10.4 (ESRI, Redlands, CA, USA) and SAGA-GIS 4.2 (Conrad et al. 2015). Variable restriction Potential trial planting sites in New Zealand should have similar environmental conditions to the native habitat of E. bosistoana in Australia, such that the species has a reasonable chance of surviving and providing growth and yield. Variable value restriction was used to restrict Category Raw data Unit Source Data type Resolution Climate Monthly mean temperature °C WorldClim Version2 Raster 1 × 1 km Monthly minimum temperature °C Monthly maximum temperature °C Monthly precipitation mm Soil Potential rooting depth m New Zealand Land Resource Inventory (NZLRI) Vector N/A Soil pH N/A Soil salinity % Soil temperature regime Classified (*) Profile available water in soil mm Topography Digital Elevation Model (DEM) m Landcare Research Raster 25 × 25 m Boundaries Territorial Authority (2016 Generalized Version) N/A Statistics New Zealand Vector N/A NZ Area Units (2015 Yearly Pattern) N/A Statistics New Zealand Vector N/A Land Cover Database version 4.1 N/A Landcare Research Vector N/A TABLE 1: GIS data used for habitat modelling including their sources, type of data and spatial resolution (*) Soil temperature regime classes include: T - thermic; WM - warm mesic; MM - mild mesic; CM - cool mesic; DM - cold mesic and C – cryic (Webb & Wilson 1995). the potential planting sites in New Zealand to those having similar environmental conditions to those in the natura range of E. bosistoana. There were two reasons for restricting the site availability. First, some areas had current land uses that were not consistent with plantation forestry (e.g. settlement area, horticulture, and indigenous forests). These areas were identified using the New Zealand Land Cover Database (LCDB v4.1) and subsequently excluded from the analysis. Second, areas were also excluded if they had environmental conditions that deviated considerably from the native habitat of the target species (e.g. permanent ice or extremely low precipitation that the species could not survive). Restricting the study area based on environmental conditions considered the environmental conditions appearing in the natural habitat of E. bosistonana in south-eastern Australia, those of existing trial plantations throughout NZ, values found in previous studies (Boland et al. 2006; Grieve et al. 1999; Webb & Wilson 1995) and, where necessary, expert knowledge (EG Mason, pers. comm.). Ranges for each of the 17 environmental variables were inferred from one or more of these four sources; the specifics are detailed below. We obtained data, including 1596 recorded points of E. bosistoana occurrence from its natural range in south- eastern Australia (ALA 2017). The occurrence locations were overlaid with environmental attribute data from the Atlas of Living Australia (ALA), WorldClim, and CSIRO Ecosystem Sciences (ALA 2017) to identify the ranges of environmental conditions (e.g. temperature, precipitation) that E. bosistoana tolerates in its native habitat. In addition to the natural range data, the species was successfully established in 30 sites throughout New Zealand. These trial plots provided additional information about the range of conditions in which the species could survive. Le and Morgonroth New Zealand Journal of Forestry Science (2020) 50:9 Page 5 Variable category Variable ID Variable Unit Source 1. Climate V1.1 Annual average temperature °C WorldClim Version2 (the determination of cold, hot, dry, and wet seasons based on the publication by Leathwick et al. (1998)) V1.2 Average monthly minimum temperature of the cold season °C V1.3 Average monthly maximum temperature of the hot season °C V1.4 Annual precipitation mm/year V1.5 Monthly precipitation of the dry season mm 2. Soil V2.1 Potential rooting depth m New Zealand Land Resource Inventory (NZLRI) V2.2 Soil pH None V2.3 Soil salinity % V2.4 Soil temperature regime Classified V2.5 Profile available water in soil mm 3. Topography V3.1 Elevation m Landcare Research V3.2 Slope degree Derived from DEM (Burrough & McDonell 1998) V3.3 Aspect degree Derived from DEM (Burrough & McDonell 1998) V3.4 Curvature None Derived from DEM (Zeverbergen & Thorne 1987) V3.5 Terrain Ruggedness Index - TRI None Derived from DEM (Riley et al. 1999) V3.6 Topographic Wetness Index - TWI None Derived from DEM (Beven & Kirkby 1979) V3.7 Wind Exposition Index None Derived from DEM (Gerlitz et al. 2015) TABLE 2: The environmental variables of the habitat model Finally, a number of previous studies on eucalypts provided guidance to help determine the thresholds for environmental conditions that the species can survive (Boland et al. 1984; Grieve et al. 1999). In brief, the three sources of data together contributed to determining the recorded value range of each variable describing the habitat of E. bosistoana. To avoid being overly restrictive (i.e. causing loss of potential habitat area) and to create chances for the species to adapt to conditions in new areas that are marginally outside the current habitat conditions, the recorded value range was generally buffered by positive and negative 10% to create the “restriction value range” for each variable. Some exceptions to buffering the recorded value range existed. The study did not apply an upper limit for the restriction value range of the variable Average monthly min temperature of the cold season as higher values of this variable would be better for the development of the target species (Millen et al. 2016). In practice, this restriction value range also contained values in Australia’s conditions that did not exist within the study area (i.e. New Zealand). An interference of the restriction value range in association with the range of values available in the study area was used to produce the “value range of interest” for each variable. Once the value range of interest was defined for each variable layer, all pixels with values outside the range of interest, or those with unsuitable land covers, were replaced by no-data pixels, such that they would be excluded from subsequent analyses. Variable stratification and standardisation The environmental raster layers used in this study contained continuous values. Even though the full value range for any given environmental variable was restricted by the value ranges of interest, no sample plot system could effectively provide sufficient samples to cover every value in those ranges. The stratified random sampling approach allowed the study to conduct an analysis on a relatively small number of groups of values instead of a wide range containing continuous values. Following the stratified random sampling approach, the value range of interest for each variable was stratified into non-overlapping compartments, also called strata. A stratum should include values reflecting similar conditions (e.g. cold weather or uneven surface) or a particular level of characteristics of the conditions at the corresponding locations (e.g. low slope or high soil pH). Defining strata requires an understanding of the full range of values present for each variable, and then specifying breakpoints within the data; the breakpoints act simultaneously as an upper limit for one stratum and a lower limit for the next stratum. Breakpoints for some of the variables used in this study were based on existing published values (Webb and Wilson 1995; Newsome et al. 2008), while others were based on an assessment of the frequency distribution of a given variable. The stratification of soil variables, slope and aspect were carried out mainly based on descriptions of New Zealand landscapes by Webb and Wilson (1995) and Newsome et al. (2008). For example, the soil pH value range, in general including values from 0 to 14, was divided into six strata: very low (<4.9), low (4.9–5.4), moderately low (5.5–5.7), near neutral (5.8– 6.4), moderately high (6.5–7.5) and high (>7.5). For variables where published strata were not available, an alternative approach was required. Climatic variables and topographic variables (other than slope and aspect), were stratified following the Jenks Natural Breaks algorithm (Jenks & Caspall 1971). This stratification uses the frequency distribution of data for a particular variable and identifies natural groupings inherent in the data. Breakpoints are identified to achieve groups with similar numbers of observations and to maximise the differences between strata. In this way, the variables were divided into strata that have relatively big differences in the data values (De Smith et al. 2015). After stratification, values of variables from existing plot locations were extracted and then assigned into strata to count the frequency of occurrence within each stratum (i.e. the number of times that values from existing plots appear in each stratum). The study used this frequency as an indicator for priority. The priority indicator increases with decreasing frequency in a stratum. The lower the frequency, the greater the need for new plots in that stratum. The next step was to calculate a normalised priority indicator for each stratum. This ensured that all priorities were normalised on a scale of 0–100, using Equation (1): pj = 100 × (1 – fj fmax -1) Equation (1) where pj and fj are the normalised priority and the frequency of the stratum j respectively, fmax is the highest frequency in the variable. For example, assuming a variable was divided into three strata with the corresponding frequencies: f1 = 200, f2 = 160 and f3 = 20. Thus, f1 = 200 is the maximum frequency and fmax = f1 = 200. Then, the normalised priority for the strata was calculated following Equation (1): p1 = 0, p2 = 20 and p3 = 90 respectively. This result was interpreted such that stratum 3, with normalised priority value of 90, had the highest priority. The example demonstrated that this normalisation enables a quantitative evaluation of priority in which the higher normalised priority values meant the higher priority to establish forest plots on locations with attribute values within the stratum. After the calculation of normalised priorities for all strata in all variables, the map layers were reclassified such that every pixel received the normalised priority value of the stratum which contained the pixel value. In other words, each map layer of a variable was converted to a priority layer that highlighted high-priority pixels on the map in respect to the particular environmental variable. Cartographic modelling The last step was to produce the map of priority index (i.e. final priority point) for each pixel within the feasible area. This index was calculated on each pixel as the weighted sum of normalised priority values Le and Morgonroth New Zealand Journal of Forestry Science (2020) 50:9 Page 6 from all variables by the priority function as shown in Equation (2): Equation (2) where P was the priority index, n was the number of environmental variables, pi was the normalised priority value of the variable i, and wi was the weight coefficient of the variable i. In particular, the priority index was calculated for each pixel by a weighted sum of, for example, 17 normalised priority values corresponding to the 17 selected variables (i.e. n = 17). The weighted sum was calculated by multiplying each normalised priority layer by its weight coefficient and then summing the weighted normalised priority values from all variable grids using an algebraic overlay analysis. In this stage, the coefficient of 1.5 was used for the four variables (i.e. annual average temperature, precipitation, soil pH and elevation) and a coefficient of 1 was used for all other variables. With each normalised priority value ranging between 0–100, the priority index in this example could have values between 0–1,900. Le and Morgonroth New Zealand Journal of Forestry Science (2020) 50:9 Page 7 Results Model structure After all the variable layers were created, we estimated the dependence between model variables in each category to detect any multi-collinearity. VIF values of all variables to identify multi-collinearity are shown in Table 3. The VIFs of variables 3.2 (Slope) and 3.5 (Terrain Ruggedness Index – TRI) considerably exceeded 10.0 (i.e. approximately 31.2 and 32.2, respectively). This was interpreted as there being a high correlation between these variables, which could lead to information redundancy and overlapping in the model. To solve this multi-collinearity problem, TRI was removed from further consideration because Slope is easier to measure/calculate and is more common in forestry research as compared to TRI. Estimation of the priority index Variable restriction Value ranges of interest of the variables are presented in Table 4. Most of the variables, except topographic aspect and soil salinity, had ranges of interest narrower than TABLE 3: Multi-collinearity analysis of each variable category Variable ID Detection- Tolerance (%) Variance Inflation Factor (VIF) Multi-collinearity Present? (A) Climatic variables Annual average temperature 82.66 1.21 FALSE Average monthly min temperature of the cold season 76.12 1.31 FALSE Average monthly max temperature of the hot season 72.91 1.37 FALSE Annual precipitation 35.50 2.82 FALSE Monthly precipitation of the dry season 41.56 2.40 FALSE (B) Soil variables Potential rooting depth 89.90 1.11 FALSE Soil pH 88.98 1.12 FALSE Soil salinity 96.11 1.04 FALSE Soil temperature regime 89.59 1.12 FALSE Profile available water in soil 97.64 1.02 FALSE (C) Topographic variables Elevation 29.29 3.41 FALSE Slope 3.21 31.18 TRUE Aspect 54.80 1.82 FALSE Curvature 66.67 1.50 FALSE Terrain Ruggedness Index - TRI 3.20 31.24 TRUE Topographic Wetness Index - TWI 85.19 1.17 FALSE Wind Exposition Index 42.43 2.36 FALSE Le and Morgonroth New Zealand Journal of Forestry Science (2020) 50:9 Page 8 Va ri ab le R an ge fr om n at iv e ha bi ta t R an ge fr om ex is ti n g P SP s Fr om re fe re n ce s R ec or de d va lu e ra n ge R es tr ic ti on va lu e ra n ge St ud y ar ea r an ge R an ge o f in te re st A nn ua l a ve ra ge te m pe ra tu re 9. 00 – 1 8. 00 10 .1 0 – 14 .1 0 9. 00 – 1 8. 00 8. 10 – 1 8. 90 -2 .5 5 – 16 .7 9 8. 10 – 1 6. 79 Av er ag e m on th ly m in te m pe ra tu re of th e co ld s ea so n -2 .0 – 8 .0 0. 5 – 5. 7 N o up pe r lim it -2 .0 – 8 .0 ≥ -3 .0 -1 1. 93 – 1 2. 2 -3 .0 – 1 2. 2 Av er ag e m on th ly m ax te m pe ra tu re of th e ho t s ea so n 22 .2 0 – 29 .8 0 18 .9 0 – 23 .3 0 18 .9 0 – 29 .8 0 18 .1 4 – 30 .5 6 5. 97 – 2 4. 87 18 .1 4 – 24 .8 7 A nn ua l p re ci pi ta ti on 63 6 – 19 16 69 3 – 18 77 63 6 – 19 16 50 8 – 20 44 39 2 – 68 07 50 8 – 20 44 M on th ly p re ci pi ta ti on o f t he d ry se as on 39 .0 – 1 00 .0 42 .0 – 1 35 .0 39 .0 – 1 35 .0 32 .9 – 1 41 .1 28 .0 – 6 39 .0 32 .9 – 1 41 .1 Po te nt ia l r oo ti ng d ep th 0. 40 – 1 .4 0 0. 29 – 1 .3 5 0. 29 – 1 .4 0 0. 19 – 1 .5 0 – 1. 35 0. 19 – 1 .3 5 So il pH 4. 0 – 5. 0 5. 2 – 8. 0 4. 0 – 8. 0 3. 9 – 8. 1 4. 7 – 8. 0 4. 7 – 8. 0 So il sa lin it y N /A = 0. 02 ≤ 2. 90 * 0. 02 – 2 .9 0 0. 02 – 2 .9 0 0. 02 – 0 .8 5 0. 02 – 0 .8 5 So il te m pe ra tu re r eg im e N /A CM – W M ≥ CM * * CM – W M CM – T C – T CM – T Pr of ile a va ila bl e w at er in s oi l 32 .0 – 1 99 .0 35 .0 – 3 00 .0 32 .0 – 3 00 .0 15 .3 – 3 16 .7 0 – 30 0. 0 15 .3 – 3 00 .0 El ev at io n 0 – 11 96 .0 15 .0 – 6 40 .0 0 – 11 96 .0 -1 19 .6 – 1 31 5. 6 -6 0. 0 – 37 37 .0 -6 0 – 13 15 .6 Sl op e 0 – 25 .5 7 0 – 33 .6 2 0 – 33 .6 2 -2 .5 6 – 36 .1 7 0 – 87 .3 2 0 – 36 .1 7 A sp ec t -1 – 3 59 .9 6 -1 – 3 51 .5 3 -1 – 3 59 .9 6 -3 7. 10 – 3 60 .0 0 -1 – 3 59 .9 4 -1 – 3 59 .9 4 Cu rv at ur e -7 .5 4 – 36 .8 2 -3 .8 4 – 7. 52 -7 .5 4 – 36 .8 2 -1 1. 98 – 4 1. 26 -3 36 .0 0 – 13 50 .4 0 -1 1. 98 – 4 1. 26 To po gr ap hi c W et ne ss In de x - T W I 1. 00 2 – 15 .8 05 7. 30 4 – 16 .0 12 1. 00 2 – 16 .0 12 -0 .4 78 – 1 7. 49 2 3. 37 7 – 26 .4 21 3. 37 7 – 17 .4 92 W in d Ex po si ti on In de x 0. 72 2 – 1. 31 3 0. 81 5 – 1. 26 5 0. 72 2 – 1. 31 3 0. 66 3 – 1. 37 2 0. 75 0 – 1. 35 3 0. 75 0 – 1. 35 3 TA B LE 4 : R es ul ts o f v ar ia bl e re st ri ct io n *( Gr ie ve e t a l. 19 99 ; W eb b & W ils on 1 99 5) . ** (B ol an d et a l. 19 84 ; N ZD FI 2 01 5) . the corresponding value ranges for the study area (i.e. all the values appeared in the study area). These conditions together determined available areas for the expansion of the existing PSPs (Fig. 2). In the variable value restriction process, areas to be excluded from the variables were either areas with values out of the value range of interest or areas under undesirable land covers (e.g. urban area, indigenous for- est and mangrove area). Through this process, approxi- mately 71.25% of the study area (i.e. 189,055.7 km2), primarily in the South Island, was excluded, including high mountainous areas, indigenous forests, and a large natural reserve area in the North Island (Fig. 2). Variable stratification and standardisation The results of variable stratification are in Table 5. From the counted frequency of each stratum, the normalised priority values of strata in each variable were calculated as in Table 6. Overlay analysis The final result was a map representing priority index values over the study area (Fig. 3). The maximum and minimum of the priority index were recorded at 1,506.5 and 20.5 respectively, with high values indicating areas with under-represented environmental conditions amongst the existing PSP network. The locations of the 1095 existing PSPs were also added to the map. Obviously, these locations distributed among low-priority areas as their environmental characteristics were well covered. In general, high index areas were mainly in high elevation areas where the conditions in terms of the three types of variables were quite different than those in existing PSPs. The highest priority areas for E. bosistoana were in Rangitikei District (i.e. Moawhango) and Taupo District (i.e. Broadlands, Rangitaiki, and Tongariro). The other high priority zones included Northland, Auckland and Gisborne regions, and southeast-facing hillsides of the mountain chains in central South Island. Discussion The main objective of the model in this paper was to detect locations for new PSPs, which would enhance the effectiveness of the current PSP network for Eucalyptus bosistoana in New Zealand. The modelling results successfully allowed us to identify areas for a strategic forest plot expansion in New Zealand with the greatest potential to provide valuable information for site-species interaction. The capabilities of GIS were used in this study in association with the stratified random sampling method to create a flexible habitat model that was used to identify Le and Morgonroth New Zealand Journal of Forestry Science (2020) 50:9 Page 9 FIGURE 2: Map of the study area BEFORE (A) and AFTER (B) variable restriction (value pixel: feasible for plot expansion; non-value pixel: excluded by variable restriction). The map inset in Figure 2A provides an example of the large number of sub-plots within any given point on the map. Le and Morgonroth New Zealand Journal of Forestry Science (2020) 50:9 Page 10 Va ri ab le R an ge o f i n te re st M in . va lu e B re ak po in t ( up pe r li m it o f t he s tr at um ) R ef er en ce s fo r st ra ti fi ca ti on S1 S2 S3 S4 S5 S6 S7 S8 S9 A nn ua l a ve ra ge te m pe ra tu re 8. 1 – 16 .7 9 8. 1 8. 93 9. 68 10 .4 6 11 .2 4 12 .0 0 12 .7 7 13 .6 1 14 .5 6 16 .7 9 N at ur al B re ak s Av er ag e m on th ly m in te m pe ra tu re of th e co ld s ea so n -3 – 1 2. 2 -3 -1 .6 7 -0 .4 0 0. 80 2. 03 3. 23 4. 47 5. 83 7. 57 12 .2 0 Av er ag e m on th ly m ax te m pe ra tu re of th e ho t s ea so n 18 .1 4 – 24 .8 7 18 .1 4 18 .8 0 19 .4 7 20 .1 0 20 .7 3 21 .3 3 21 .9 0 22 .5 0 23 .1 3 24 .8 7 A nn ua l p re ci pi ta ti on 50 8 – 20 44 50 8 73 0 90 1 10 67 12 25 13 77 15 22 16 78 18 50 20 44 M on th ly p re ci pi ta ti on o f t he d ry se as on 32 .9 – 1 41 .1 32 .9 51 .0 61 .0 71 .0 82 .0 93 .0 10 4. 0 11 5. 0 12 7. 0 14 1. 1 Po te nt ia l r oo ti ng d ep th 0. 19 – 1 .3 5 0. 19 0. 25 0. 45 0. 6 0. 9 1. 2 1. 35 (N ew so m e et a l. 20 08 ; W eb b & W ils on 1 99 5) So il pH 4. 7 – 8 4. 7 4. 9 5. 5 5. 8 6. 5 7. 6 8. 0 So il sa lin it y 0. 02 – 0 .8 5 0. 02 0. 05 0. 15 0. 3 0. 7 0. 85 So il te m pe ra tu re r eg im e CM – T CM CM M M W M T Pr of ile a va ila bl e w at er in s oi l 15 .3 – 3 00 15 .3 30 .0 0 60 90 15 0 25 0 30 0 El ev at io n -6 0 – 13 15 .6 -6 0 11 0 23 5 36 5 50 0 64 1 79 1 95 1 11 25 13 15 N at ur al B re ak s Sl op e 0 – 36 .1 74 27 0 4. 00 8. 00 12 .0 0 16 .0 0 21 .0 0 26 .0 0 36 .1 8 (W eb b & W ils on 19 95 ) A sp ec t -1 – 3 59 .9 4 -1 Fl at N N E E SE S SW W N W (* ) Cu rv at ur e -1 1. 97 6 – 41 .2 56 -1 1. 97 6 -3 .5 2 -1 .1 2 0. 32 1. 76 3. 36 5. 12 7. 04 10 .0 8 41 .2 6 N at ur al B re ak s To po gr ap hi c W et ne ss In de x - T W I 3. 37 7 – 17 .4 92 3. 37 7 7. 25 8. 06 8. 87 9. 77 10 .7 6 12 .0 2 13 .4 6 15 .1 7 17 .5 0 W in d Ex po si ti on In de x 0. 75 1 – 1. 35 0 0. 75 1 0. 82 0. 88 0. 95 1. 02 1. 09 1. 15 1. 22 1. 29 1. 36 TA B LE 5 : R es ul ts o f v ar ia bl e st ra ti fic at io n (* ) St ra ti fie d as fl at s ur fa ce a nd 8 c ar di na l-i nt er ca rd in al d ir ec ti on s: F la t ( -1 ); N ( 0 – 22 .5 ; 3 37 .5 – 3 59 .9 9) ; N E (2 2. 5 – 67 .5 ); E ( 67 .5 – 1 12 .5 ); S E (1 12 .5 – 1 57 .5 ); S ( 15 7. 5 – 20 2. 5) ; S W ( 20 2. 5 – 24 7. 5) ; W ( 24 7. 5 – 29 2. 5) ; N W (2 92 .5 – 3 37 .5 ). Le and Morgonroth New Zealand Journal of Forestry Science (2020) 50:9 Page 11 Va ri ab le St ra tu m 1 St ra tu m 2 St ra tu m 3 St ra tu m 4 St ra tu m 5 St ra tu m 6 St ra tu m 7 St ra tu m 8 St ra tu m 9 f p f p f p f p f p f p f p f p f p A nn ua l a ve ra ge te m pe ra tu re 0 10 0 0 10 0 2 10 0 0 10 0 38 4 24 50 7 0 19 9 61 3 99 0 10 0 Av er ag e m on th ly m in te m pe ra tu re o f t he c ol d se as on 0 10 0 0 10 0 2 99 37 0 2 37 8 0 20 6 46 13 9 63 0 10 0 0 10 0 Av er ag e m on th ly m ax te m pe ra tu re o f t he h ot s ea so n 0 10 0 4 99 79 78 85 76 35 4 0 20 2 43 35 4 0 11 97 6 98 A nn ua l p re ci pi ta ti on 27 3 24 35 9 0 27 6 23 15 6 57 16 96 2 99 8 98 3 99 2 99 M on th ly p re ci pi ta ti on o f t he dr y se as on 25 0 44 36 7 17 44 4 0 2 10 0 12 97 9 98 6 99 0 10 0 5 99 Po te nt ia l r oo ti ng d ep th 0 10 0 51 8 0 15 0 71 37 9 27 9 98 39 92 - - - So il pH 0 10 0 44 8 9 23 95 49 5 0 12 7 74 2 10 0 - - - So il sa lin it y 10 95 0 0 10 0 0 10 0 0 10 0 0 10 0 - - - - So il te m pe ra tu re r eg im e 13 0 84 80 9 0 15 6 81 0 10 0 - - - - - Pr of ile a va ila bl e w at er in s oi l 0 10 0 35 2 39 12 3 79 57 5 0 43 93 2 10 0 - - - El ev at io n 48 7 0 21 7 55 37 6 23 13 97 2 10 0 0 10 0 0 10 0 0 10 0 0 10 0 Sl op e 46 7 0 15 0 68 18 2 61 14 6 69 11 7 75 31 93 2 10 0 - - A sp ec t 10 7 57 52 79 24 6 0 12 6 49 23 91 16 3 34 96 61 22 7 8 55 78 Cu rv at ur e 1 10 0 16 98 87 2 0 18 7 79 13 99 5 99 0 10 0 1 10 0 0 10 0 To po gr ap hi c W et ne ss In de x - TW I 0 10 0 10 7 53 22 7 0 17 5 23 19 0 16 12 9 43 21 8 4 48 79 1 10 0 W in d Ex po si ti on In de x 1 10 0 14 97 29 7 26 40 2 0 21 4 47 51 87 11 2 72 4 99 0 10 0 TA B LE 6 : F re qu en cy ( f) a nd n or m al is ed p ri or it y va lu es ( p) fo r ea ch c om bi na ti on o f v ar ia bl e an d st ra tu m . priority locations for PSP expansion. GIS was used here as a platform for contributing to the decision making process via variable data management, variable layer production, normalised priority calculation by means of spatial analysis, a weighted-overlaying combination of variables, and finally, mapping by means of cartographic modelling. The data processing and output of the spatial results prove the potential of GIS-based modelling in generating information for large-scale studies. The most obvious advantage of the method used in this study is the possibility to generate priority indices for forest expansion in such a huge area as a whole country (i.e. New Zealand nationwide) consuming reasonable time and labor resources. This time and money-saving method enabled the study to adopt available environmental data in combination with expert knowledge to build a habitat model for a species to be strategically planted. Performing analyses in the habitat model allowed the study to spatially utilise normalised priority functions when dealing with different species and different existing forest plot systems in the same study area even if the area was huge with diverse environmental conditions. In the context of a new habitat for the species, the study set the habitat restriction based on actual occurrences of the species in its natural habitat and in the established sample plots within the study area. However, every species has its own tolerance ability to adapt to a wide range of environmental conditions. In this case, there were insufficient studies on the plasticity of the species, especially in respect to each environmental factor. We extended the range of each variable by 10% to reflect the potential for the species to adapt to a new habitat. However, the example of radiata pine (i.e. Pinus radiata) raised the issue of how much adaptable E. bosistoana could be. Indeed, radiata pine was introduced to New Zealand in the twentieth century and it has developed well in a much broader range of conditions than those present in its natural habitat (Mead 2013; MPI 2016; Weston 1957). Using an inappropriate value for the assumption of species’ plasticity may reduce the probability of identifying correct ecological thresholds. As a consequence, it is possible that variable buffering could result in either omitting significant suitable areas of potential new habitat (i.e. too small a buffer), or in adding unsuitable areas (i.e. too large a buffer) to the model. In this case study, the approach involved an assumption that the model input data was error-free. However, the evaluation processes used types of information that are often uncertain and imprecise. These problems may arise from errors from data measurement and processing (Fischer & Wang 2011; Karger et al. 2017; Pielou 1984). Most of the information to create variable layers was collected from various providers who have different ways of managing data, and the rest resulted from spatial analysis. These sources of information had the most significant influence on the model quality because their errors were almost systematic and propagated all over the study area (M. P. Austin & Heyligers 1991; Pielou 1984; Sellars & Jolls 2007). For example, the study by Pearse et al. (2015) highlighted inaccuracies in the soil data used in the present study. However, the model in this paper had to use that soil data because it is the best dataset in existence for the study area. For GIS-based modelling, it is possible to analyze the sensitivity of the model’s results to the uncertainties and uncertainty propagation of data. The solution should consider the analytical error propagation method suggested by Store and Kangas (2001) and the Monte Carlo simulation described by Burrough and McDonell (1998). The information from existing plots may also be affected by GPS inaccuracy, but this type of errors were less than 10 m (NZDFI 2015) that were not crucial at the large scale of the study area. In general, the case study used a new flexible technique that can be applied to a wide range of contexts. The habitat model was built as a description of environmental conditions across New Zealand. Similarly, with the development of global environmental data by remote sensing, such habitat models can be easily built for other parts of the world. Obviously, these models can be widely used for the expansion of any plot network of any species to be introduced to a new habitat. Although Le and Morgonroth New Zealand Journal of Forestry Science (2020) 50:9 Page 12 FIGURE 3: The priority index map following the same procedure, each study subject requires a specific set of variables and its own criteria for variable restriction as well as particular methods for variable stratification. Another potential use of this technique is to identify locations for the expansion of forest inventory plots, especially in natural forests. Indeed, environmental conditions among natural forests were unevenly distributed which leads to heterogeneous forest compositions and structures. Following the habitat modelling approach, an inventory plot system after expansion can better cover the whole range of environmental habitats of the forest statuses. This is crucial to achieve more precise and more accurate inventory result over the region of interest. Conclusions The results indicate that it was inappropriate to plant E. bosistoana in some parts of the South Island, such as the rainforest areas in Westland, areas too high and too cold deep in the south of the island. The existing PSP network for the species is over-represented with environmental conditions present in low-elevation New Zealand dry lands, which are located alongside the east coast of the South Island, and the southern part of the North Island. Moreover, high priority areas for further trial plots of this species included several large regions in the North Island, such as Northland, Auckland and Gisborne regions and especially some smaller parts in the center of the North island in Taupo and Rangitikei Districts. Plantations of this species should also be tested at higher elevation (e.g. in mountainous areas of the Canterbury, Marlborough and Tasman regions). We have developed a new approach involving a priority index to determine strategic expansion of forest monitoring plots that uses GIS-based models in association with stratified random sampling. The technique used in this study was illustrated by a case study that successfully identified optimal sites for new establishments of forest plots of Eucalyptus bosistoana in New Zealand. New plots in these sites, if established, will provide crucial information for the site-species matching programme of the NZDFI. We hope that, by describing the methodology in this study, a broader base will help forest resource professionals and researchers able to build GIS-based habitat models and apply these models in creating more adequate and efficient plot network designs to monitor and assess forests and the surrounding conditions as well as relationships between them. For further studies, we recommend using more variables to better describe the environment surrounding PSPs, if possible, especially in terms of soil conditions that critically influence a plant’s growth. However, it is important to find references determining the restriction ranges regarding the new variables to be added. Moreover, trial plantings of the target species should be established in areas with marginal environmental conditions to identify the species’ plasticity. This could provide more significant information to enhance the performance of the model. The model in this paper can also be supported by statistical methods to determine appropriate weight coefficients for the model variables in consideration. It is recommended to use the analytical hierarchy process method, which is one of the most popular methods to obtain variables’ weights in GIS- based modelling (Carver 1991; Chen et al. 2010; Marinoni et al. 2009). Competing interests The authors declare that they have no competing interests. Authors’ contributions TSL analysed the data and wrote the manuscript under the supervision of JM. Both authors read and approved the final version of the manuscript. Acknowledgements The authors are thankful to the NZDFI for the permission to use the data of PSPs. They also appreciate the useful discussions and knowledge provided by Professor Euan Mason. Finally, the authors are grateful to the reviewers and editor whose comments improved the manuscript. References Aljaryian, R., Kumar, L., & Taylor, S. (2016). Modelling the current and potential future distributions of the sunn pest Eurygaster integriceps (Hemiptera: Scutelleridae) using CLIMEX. Pest Management Science, 72(10), 1989-2000. https://doi. org/10.1002/ps.4247 Apiolaza, L., McConnochie, R., Millen, P., Van Ballekom, S., & Walker, J.C.F. (2011). Introducing durable species to New Zealand drylands: Genetics of early adaption of Eucalyptus bosistoana. Paper presented at the Developing a Eucalypt Resource-Learning from Australia and elsewhere, Blenheim Workshop proceedings, Wood Technology Research Centre. Austin, J.E., Buhl, T.K., Guntenspergen, G.R., Norling, W., & Sklebar, H.T. (2001). Duck populations as indicators of landscape condition in the prairie pothole region. Environmental Monitoring and Assessment, 69(1), 29-47. https://doi. org/10.1023/A:1010748527667 Austin, M.P., & Heyligers, P.C. (1991). Vegetation survey design, a new approach: Gradsect sampling. In: C. R. Margules & M. P. Austin (Eds.), Nature Conservation: Cost Effective Biological Surveys and Data Analysis (pp. 31-36). Melbourne, Australia: CSIRO. Barringer, J.R.F., Pairman, D., & McNeill, S.J. (2002). Development of a high-resolution Digital Elevation Model for New Zealand. Retrieved from https:// lris.scinfo.org.nz/document/9213/download/ Le and Morgonroth New Zealand Journal of Forestry Science (2020) 50:9 Page 13 https://doi.org/10.1002/ps.4247 https://doi.org/10.1002/ps.4247 https://doi.org/10.1023/A:1010748527667 https://doi.org/10.1023/A:1010748527667 https://lris.scinfo.org.nz/document/9213/download/ https://lris.scinfo.org.nz/document/9213/download/ Basir. (2014). GIS-based approach to participatory land suitability analysis for tree plantations. (Doctor of Philosophy in Natural Resources and Environmental Sciences), University of Illinois, Urbana-Champaign. Beven, K.J., & Kirkby, M.J. (1979). A physically based, variable contributing area model of basin hydrology. Hydrological Sciences Bulletin, 24(1), 43-69. https://doi.org/10.1080/02626667909491834 Boland, D.J., Brooker, M.I.H., Chippendale, G.M., Hall, N., Hyland, B.P.M., Johnston, R.D., Kleinig, D.A., McDonald, M.W., & Turner, J.D. (Eds.) (2006). Forest Trees of Australia (5th ed). Melbourne, Australia: CSIRO. https://doi.org/10.1071/9780643069701 Boland, D.J., Brooker, M.I.H., Chippendale, G.M., Hall, N., Hyland, B.P.M., Johnston, R.D., Kleinig, D.A., & Turner, J.D. (Eds.) (1984). Forest Trees of Australia. Melbourne, Australia: CSIRO. Burrough, P.A., & McDonell, R.A. (1998). Principles of Geographical Information Systems. New York: Oxford University Press. Carver, S.J. (1991). Integrating multi-criteria evaluation with geographical information systems. International Journal of Geographical Information System, 5(3), 321-339. https://doi. org/10.1080/02693799108927858 Cawsey, E.M., Austin, M.P., & Baker, B.L. (2002). Regional vegetation mapping in Australia: a case study in the practical use of statistical modelling. Biodiversity and Conservation, 11(12), 2239-2274. https://doi. org/10.1023/A:1021350813586 Chen, Y., Yu, J., & Khan, S. (2010). Spatial sensitivity analysis of multi-criteria weights in GIS-based land suitability evaluation. Environmental Modelling & Software, 25(12), 1582-1591. https://doi. org/10.1016/j.envsoft.2010.06.001 Conrad, O., Bechtel, B., Bock, M., Dietrich, H., Fischer, E., Gerlitz, L., Wehberg, J., Wichmann, V., & Böhner, J. (2015). System for Automated Geoscientific Analyses (SAGA) v. 2.1.4. Geoscientific Model Development, 8(7), 1991-2007. https://doi. org/10.5194/gmd-8-1991-2015 Danz, N.P., Regal, R.R., Niemi, G.J., Brady, V.J., Hollenhorst, T., Johnson, L.B., Host, G.E., Hanowski, J.M., Johnston, C.A., Brown, T., Kingston, J., & Kelly, J.R. (2003). Environmentally stratified sampling design for the development of Great Lakes environmental Indicators. Environmental Monitoring and Assessment, 102(1), 41-65. https:// doi.org/10.1007/s10661-005-1594-8 De Smith, M.J., Goodchild, M.F., & Longley, P.A. (2015). Geospatial Analysis. A Comprehensive Guide to Principles, Techniques and Software Tools (5th ed). Kibworth Harcourt, UK: Troubador Publishing Ltd. Dettmers, R., & Bart, J. (1999). A GIS modeling method applied to predicting forest songbird habitat. Ecological Applications, 9(1), 152-163. https://doi. org/10.1890/1051-0761(1999)009[0152:AGMM AT]2.0.CO;2 Esfahani, M.S., & Dougherty, E.R. (2014). Effect of separate sampling on classification accuracy. Bioinformatics, 30(2), 242-250. https://doi. org/10.1093/bioinformatics/btt662 Fick, S.E., & Hijmans, R.J. (2017). WorldClim 2: new 1-km spatial resolution climate surfaces for global land areas. International Journal of Climatology, 37(12), 4302-4315. https://doi.org/10.1002/joc.5086 Fischer, M.M., & Wang, J. (2011). Spatial Data Analysis: Models, Methods and Techniques. Berlin/ Heidelberg, Germany: Springer Science & Business Media. García, C., García, J., López Martín, M., & Salmerón, R. (2014). Collinearity: revisiting the variance inflation factor in ridge regression. Journal of Applied Statistics, 42(3), 648-661. https://doi.org/ 10.1080/02664763.2014.980789 Gerlitz, L., Conrad, O., & Bohner, J. (2015). Large-scale atmospheric forcing and topographic modification of precipitation rates over High Asia - a neural- network-based approach. Earth System Dynamics, 6(1), 61-81. https://doi.org/10.5194/esd-6-61- 2015 Goulding, C.J. (2005). Measurement of trees. In M. Colley (Ed.), The NZIF Forestry Handbook (4th ed.). Christchurch, New Zealand: New Zealand Institute of Forestry. Green, R.H. (1979). Sampling Design and Statistical Methods for Environmental Biologists. Hoboken, New Jersey, USA: John Wiley & Son. Grieve, C.M., Guzy, M.R., Poss, J.A., & Shannon, M.C. (1999). Screening Eucalyptus clones for salt tolerance. Horticulture Science, 34(5), 867-870. https://doi. org/10.21273/HORTSCI.34.5.867 Hill, M.P., Bertelsmeier, C., Clusella-Trullas, S., Garnas, J., Robertson, M.P., & Terblanche, J.S. (2016). Predicted decrease in global climate suitability masks regional complexity of invasive fruit fly species response to climate change. Biological Invasions, 18(1), 1105-1119. https://doi.org/10.1007/ s10530-016-1078-5 James, F.C., & Shugart, H.H. (1970). A quantitative method of habitat description. Audubon Field Notes by the National Audubon Society, 24(6), 727-736. Jenks, G.F., & Caspall, F.C. (1971). Error on Choroplethic Maps: Definition, Measurement, Reduction. Annals of the Association of American Geographers, 61(2), 217-244. https://doi. org/10.1111/j.1467-8306.1971.tb00779.x Johnston, R., & Temple, S. (1986). Assessing habitat quality for birds nesting in fragmented tallgrass prairies. In: Wildlife 2000: Modelling Habitat Relationships of Terrestrial Vertebrates, J. Verner, Le and Morgonroth New Zealand Journal of Forestry Science (2020) 50:9 Page 14 https://doi.org/10.1080/02626667909491834 https://doi.org/10.1071/9780643069701 https://doi.org/10.1080/02693799108927858 https://doi.org/10.1080/02693799108927858 https://doi.org/10.1023/A:1021350813586 https://doi.org/10.1023/A:1021350813586 https://doi.org/10.1016/j.envsoft.2010.06.001 https://doi.org/10.1016/j.envsoft.2010.06.001 https://doi.org/10.5194/gmd-8-1991-2015 https://doi.org/10.5194/gmd-8-1991-2015 https://doi.org/10.1007/s10661-005-1594-8 https://doi.org/10.1007/s10661-005-1594-8 https://doi.org/10.1890/1051-0761(1999)009%5b0152:AGMMAT%5d2.0.CO;2 https://doi.org/10.1890/1051-0761(1999)009%5b0152:AGMMAT%5d2.0.CO;2 https://doi.org/10.1890/1051-0761(1999)009%5b0152:AGMMAT%5d2.0.CO;2 https://doi.org/10.1093/bioinformatics/btt662 https://doi.org/10.1093/bioinformatics/btt662 https://doi.org/10.1002/joc.5086 https://doi.org/10.1080/02664763.2014.980789 https://doi.org/10.1080/02664763.2014.980789 https://doi.org/10.5194/esd-6-61-2015 https://doi.org/10.5194/esd-6-61-2015 https://doi.org/10.21273/HORTSCI.34.5.867 https://doi.org/10.21273/HORTSCI.34.5.867 https://doi.org/10.1007/s10530-016-1078-5 https://doi.org/10.1007/s10530-016-1078-5 https://doi.org/10.1111/j.1467-8306.1971.tb00779.x https://doi.org/10.1111/j.1467-8306.1971.tb00779.x M.L. Morrison, and C.J. Ralph, Eds. (pp. 245-249). Madison, Wisconsin, USA: University of Wisconsin Press. Karger, D.N., Conrad, O., Bohner, J., Kawohl, T., Kreft, H., Soria-Auza, R.B., Zimmermann, N.E., Linder, H.P., & Kessler, M. (2017). Climatologies at high resolution for the earth’s land surface areas. Scientific Data, 4(170122). https://doi.org/10.1038/ sdata.2017.122 Kimberley, M., West, G., Dean, M., & Knowles, L. (2005). The 300 index - A volume productivity index for radiata pine. New Zealand Journal of Forestry, 50(2), 13-18. Knollova, I., Chytry, M., Tichy, L., & Hajek, O. (2005). Stratified resampling of phytosociological databases: some strategies for obtaining more representative data sets for classification studies. Journal of Vegetation Science, 16(4), 363-372. h t t p s : / / d o i . o rg / 1 0 . 1 1 1 1 / j . 1 6 5 4 - 1 1 0 3 . 2 0 0 5 . tb02388.x Kohl, M., Magnussen, S., & Marchetti, M. (2006). Sampling methods, remote sensing and GIS multiresource forest inventory. Heidelberg, Germany: Springer- Verlag. https://doi.org/10.1007/978-3-540- 32572-7 Kutner, H.M., Nachtsheim, C.J., & Neter, J. (2003). Applied Linear Regression Models. (4th ed.) New York, USA: McGraw-Hill Companies. Leathwick, J.R., Wilson, G., & Stephens, R.T.T. (1998). Climate Surfaces for New Zealand. Landcare Research Contract Report: LC9798/126. 26 p. Retrieved from https://www.landcareresearch.co.nz/uploads/ public/Tools-And-Resources/Maps/LENZ/Climate_ Surfaces_for_New_Zealand.pdf Lin, H., Murray, T., & Mason, E. (2017). Incidence of and defoliation by a newly introduced pest, Paropsisterna variicollis (Coleoptera: Chrysomelidae), on eleven durable Eucalyptus species in Hawke’s Bay, New Zealand. New Zealand Plant Protection, 70, 45-51. https://doi.org/10.30843/nzpp.2017.70.26 Marinoni, O., Higgins, A., Hajkowics, S., & Collins, K. (2009). The multiple criteria analysis tool (MCAT): a new software tool to support environmental investment decision making. Environmental Modelling & Software, 24(2), 549-562. https://doi. org/10.1016/j.envsoft.2008.06.015 Mead, D. (2013). Sustainable management of Pinus radiata. [FAO Forestry Paper No. 170]. Rome, Italy: Food and Agriculture Organization of the United Nations. Millen, P., Burgess, J., McConnochie, R., May, R., & Buck, K. (2016). Early height growth of durable eucalypt species. The measurement and analysis of NZDFI demonstration trials planted 2010-2014. Retrieved from http://nzdfi.org.nz/wp-content/ uploads/2014/09/NZDFI-ReportonAssessmentof- 2 0 1 0 - 1 4 - D u r a b l e E u c a l y p t D e m o T r i a l s - 30June2016.pdf MPI. (2016). National Exotic Forest Description as at 1 April 2016. Retrieved from http://www. n o r t h l a n d w o o d c o u n c i l . c o . n z / d o w n l o a d s / national-exotic-forest-description-april-2016.pdf Newsome, P.F.J., Wilde, R.H., & Willoughby, E.J. (Eds.). (2008). Land Resource Information System of Spatial Data Layers. Palmerston North, New Zealand: Landcare Research New Zealand Ltd. Nicholas, I., & Millen, P. (2012). Durable Eucalypt Leaflet Series. [NZDFI project publication]. Retrieved from http://nzdfi.org.nz/wp-content/ uploads/2015/01/E-bosistoana-information- leaflet.pdf NZDFI. (2015). NZDFI Site-Species Matching Research Plan. Retrieved from http://nzdfi.org.nz/wp- content/uploads/2014/08/Section-3-NZDFI-Site_ Species-Matching-Research-Plan-.pdf NZDFI. (2017). Breeding tomorrow’s trees today. Retrieved from http://nzdfi.org.nz NZDFI. (2019). E. bosistoana: Information for Growers. Retrieved from https://nzdfi.org.nz/grower- information/growing-ground-durable-eucalypts/ growing-regimes/choosing-the-right-species/e- bosistoana-grower-information/ Pattison, R.R., & Mack, R.N. (2008). Potential distribution of the invasive tree Triadica sebifera (Euphorbiaceae) in the United States: evaluating Climex predictions with field trials. Global Change Biology, 14(4), 813-826. https://doi.org/10.1111/ j.1365-2486.2007.01528.x Pausas, J., Braithwaite, M., & Austin, M. (1995). Modelling habitat quality for arboreal marsupials in the South Coastal forests of New South Wales, Australia. Forest Ecology Management, 78(1-3), 39-49. https://doi.org/10.1016/0378-1127(95)03598-5 Pearse, G., Moltchanova, E., & Bloomberg, M. (2015). Assessment of the accuracy of profile available water and potential rooting depth estimates held within New Zealand’s fundamental soil layers geo- database. Soil Research, 53(7), 737-744. https:// doi.org/10.1071/SR14012 Pereira, J.M.O.D. (1989). A Spatial Approach to Statistical Habitat Suitability Modeling: the Mt. Graham Red Squirrel Case Study. USA: University of Arizona. Pielou, E.C. (1984). The Interpretation of Ecological Data: a Primer of Classification and Ordination. John Wiley & Sons. Prober, S.M., Potts, B.M., Bailey, T., Byrne, M., Dillon, S., Harrison, P.A., Hoffmann, A.A., Jordan, R., Mclean E.H., Steane, D.A., Stock, W.D., & Vaillancourt, R.E. (2016). Climate adaptation and ecological restoration in Eucalypts. Proceedings of the Royal Society of Victoria, 128(1), 40-53. https://doi.org/10.1071/RS16004 Le and Morgonroth New Zealand Journal of Forestry Science (2020) 50:9 Page 15 https://doi.org/10.1038/sdata.2017.122 https://doi.org/10.1038/sdata.2017.122 https://doi.org/10.1111/j.1654-1103.2005.tb02388.x https://doi.org/10.1111/j.1654-1103.2005.tb02388.x https://doi.org/10.1007/978-3-540-32572-7 https://doi.org/10.1007/978-3-540-32572-7 https://www.landcareresearch.co.nz/uploads/public/Tools-And-Resources/Maps/LENZ/Climate_Surfaces_for_New_Zealand.pdf https://www.landcareresearch.co.nz/uploads/public/Tools-And-Resources/Maps/LENZ/Climate_Surfaces_for_New_Zealand.pdf https://www.landcareresearch.co.nz/uploads/public/Tools-And-Resources/Maps/LENZ/Climate_Surfaces_for_New_Zealand.pdf https://doi.org/10.30843/nzpp.2017.70.26 https://doi.org/10.1016/j.envsoft.2008.06.015 https://doi.org/10.1016/j.envsoft.2008.06.015 http://nzdfi.org.nz/wp-content/uploads/2014/09/NZDFI-ReportonAssessmentof-2010-14-DurableEucalyptDemoTrials-30June2016.pdf http://nzdfi.org.nz/wp-content/uploads/2014/09/NZDFI-ReportonAssessmentof-2010-14-DurableEucalyptDemoTrials-30June2016.pdf http://nzdfi.org.nz/wp-content/uploads/2014/09/NZDFI-ReportonAssessmentof-2010-14-DurableEucalyptDemoTrials-30June2016.pdf http://nzdfi.org.nz/wp-content/uploads/2014/09/NZDFI-ReportonAssessmentof-2010-14-DurableEucalyptDemoTrials-30June2016.pdf http://www.northlandwoodcouncil.co.nz/downloads/national-exotic-forest-description-april-2016.pdf http://www.northlandwoodcouncil.co.nz/downloads/national-exotic-forest-description-april-2016.pdf http://www.northlandwoodcouncil.co.nz/downloads/national-exotic-forest-description-april-2016.pdf http://nzdfi.org.nz/wp-content/uploads/2015/01/E-bosistoana-information-leaflet.pdf http://nzdfi.org.nz/wp-content/uploads/2015/01/E-bosistoana-information-leaflet.pdf http://nzdfi.org.nz/wp-content/uploads/2015/01/E-bosistoana-information-leaflet.pdf http://nzdfi.org.nz/wp-content/uploads/2014/08/Section-3-NZDFI-Site_Species-Matching-Research-Plan-.pdf http://nzdfi.org.nz/wp-content/uploads/2014/08/Section-3-NZDFI-Site_Species-Matching-Research-Plan-.pdf http://nzdfi.org.nz/wp-content/uploads/2014/08/Section-3-NZDFI-Site_Species-Matching-Research-Plan-.pdf http://nzdfi.org.nz https://nzdfi.org.nz/grower-information/growing-ground-durable-eucalypts/growing-regimes/choosing-the-right-species/e-bosistoana-grower-information/ https://nzdfi.org.nz/grower-information/growing-ground-durable-eucalypts/growing-regimes/choosing-the-right-species/e-bosistoana-grower-information/ https://nzdfi.org.nz/grower-information/growing-ground-durable-eucalypts/growing-regimes/choosing-the-right-species/e-bosistoana-grower-information/ https://nzdfi.org.nz/grower-information/growing-ground-durable-eucalypts/growing-regimes/choosing-the-right-species/e-bosistoana-grower-information/ https://doi.org/10.1111/j.1365-2486.2007.01528.x https://doi.org/10.1111/j.1365-2486.2007.01528.x https://doi.org/10.1016/0378-1127(95)03598-5 https://doi.org/10.1071/SR14012 https://doi.org/10.1071/SR14012 https://doi.org/10.1071/RS16004 Reisinger, T., & Kennedy, D. (1990). A spatial decision support system for opportunity area analysis on the Jefferson national forest. Paper presented at the Annual Conference and Exposition GIS/LIS ‘90, Anaheim, Canada 5-10 November. Riley, S.J., Degloria, S.D., & Elliot, R. (1999). A Terrain Ruggedness Index that quantifies topographic heterogeneity. Intermountain Journal of Sciences, 5(1-4), 23-27. Royall, R.M. (1970). On finite population sampling theory under certain linear regression models. Biometrika, 57(2), 377-387. https://doi.org/10.1093/ biomet/57.2.377 Scott, C.T. (1998). Sampling methods for estimating change in forest resources. Ecological Applications, 8(2), 228-233. https://doi.org/10.1890/1051- 0761(1998)008[0228:SMFECI]2.0.CO;2 Sellars, J.D., & Jolls, C.L. (2007). Habitat modeling for Amaranthus pumilus: An application of Light Detection and Ranging (LIDAR) Data. Journal of Coastal Research, 23(5), 1193-1202. https://doi. org/10.2112/04-0334.1 Shabani, F., Kumar, L., & Taylor, S. (2012). Climate change impacts on the future distribution of date palms: a modeling exercise using CLIMEX. PLoS ONE, 7(10). https://doi.org/10.1371/journal.pone.0048021 Shaw, D., & Atkinson, S. (1988). GIS applications for golden-cheeked warbler habitat description. Paper presented at the Annual Conference and Exposition GIS/LIS ‘88, Falls Church, Virginia, USA. Store, R., & Kangas, J. (2001). Integrating spatial multi- criteria evaluation and expert knowledge for GIS- based habitat suitability modeling. Landscape and Urban Planning, 55(2), 79-93. https://doi. org/10.1016/S0169-2046(01)00120-7 Sutherst, R.W., & Maywald, G.F. (1985). A computerized system for matching climates in ecology. Agriculture, Ecosystems & Environment, 13(3- 4), 281-299. https://doi.org/10.1016/0167- 8809(85)90016-7 Sutherst, R.W., Maywald, G.F., & Kristicos, D.J. (2007). CLIMEX Version 3: User’s Guide. South Yarra, Australia: Hearne Scientific Software. Taylor, S., Kumar, L., Reid, N., & Kriticos, D.J. (2012). Climate change and the potential distribution of an invasive shrub, Lantana camara L. PLoS ONE, 7(10). https://doi.org/10.1371/journal.pone.0035565 Tomppo, E., Malimbwi, R., Katila, M., Mäkisara, K., Henttonen, H.M., Chamuya, N., Zahabu, E., & Otieno, J. (2014). A sampling design for a large area forest inventory: case Tanzania. Canadian Journal of Forest Research, 44(8), 931-948. https://doi. org/10.1139/cjfr-2013-0490 Wadge, G., Wislocki, A., & Pearson, E. (1993). Spatial Analysis in GIS for natural hazard assessment. In M. Goodchild, B. Parks & L. Steyaert (Eds.), Environmental Modeling with GIS (pp. 332-338). New York: Oxford University Press. Wallenius, K., Niemi, R.M., & Rita, H. (2011). Using stratified sampling based on pre-characterisation of samples in soil microbiological studies. Applied Soil Ecology, 51(1), 111-113. https://doi. org/10.1016/j.apsoil.2011.09.006 Webb, T.H., & Wilson, A.D. (1995). A Manual of Land Characteristics for Evaluation of Rural Land. [Landcare Research Science Series 10]. Lincoln, New Zealand: Landcare Research. Weston, G. (1957). Exotic Forest Trees in New Zealand. Statement prepared for the 7th British Commonwealth Forestry Conference. [Australia & New Zealand NZ Forest Service Bulletin]. Xuezhen, G., Shanyong, H., Chenyi, Z., Tao, W., Zhichun, X., & Shixiang, Z. (2018). Projecting the current and future potential global distribution of Hyphantria cunea (Lepidoptera: Arctiidae) using CLIMEX. Pest Management Science, 75(1), 160-169. https://doi. org/10.1002/ps.5083 Yves, T., & Ecker, K. (2014). Complex national sampling design for long-term monitoring of protected dry grasslands in Switzerland. Environmental and Ecological Statistics, 21(3), 453-476. https://doi. org/10.1007/s10651-013-0263-2 Zeverbergen, L.W., & Thorne, C.R. (1987). Quantitative analysis of land surface topography. Earth Surface Processes and Landforms, 12(1), 47-56. https://doi. org/10.1002/esp.3290120107 Le and Morgonroth New Zealand Journal of Forestry Science (2020) 50:9 Page 16 https://doi.org/10.1093/biomet/57.2.377 https://doi.org/10.1093/biomet/57.2.377 https://doi.org/10.1890/1051-0761(1998)008%5b0228:SMFECI%5d2.0.CO;2 https://doi.org/10.1890/1051-0761(1998)008%5b0228:SMFECI%5d2.0.CO;2 https://doi.org/10.2112/04-0334.1 https://doi.org/10.2112/04-0334.1 https://doi.org/10.1371/journal.pone.0048021 https://doi.org/10.1016/S0169-2046(01)00120-7 https://doi.org/10.1016/S0169-2046(01)00120-7 https://doi.org/10.1016/0167-8809(85)90016-7 https://doi.org/10.1016/0167-8809(85)90016-7 https://doi.org/10.1371/journal.pone.0035565 https://doi.org/10.1139/cjfr-2013-0490 https://doi.org/10.1139/cjfr-2013-0490 https://doi.org/10.1016/j.apsoil.2011.09.006 https://doi.org/10.1016/j.apsoil.2011.09.006 https://doi.org/10.1002/ps.5083 https://doi.org/10.1002/ps.5083 https://doi.org/10.1007/s10651-013-0263-2 https://doi.org/10.1007/s10651-013-0263-2 https://doi.org/10.1002/esp.3290120107 https://doi.org/10.1002/esp.3290120107