Strategic expansion of existing forest monitoring plots: a case study using a stratified GIS-based modelling approach

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


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 southeastern 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(NZDFI , 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.

Study area
The study area includes the North and South Islands of New Zealand, which cover a total of approximately 268,000 km 2 (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.

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) andProber 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 .   (Webb & Wilson 1995).

Variable restriction
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 southeastern 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.   (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): Equation (1) where p j and f j are the normalised priority and the frequency of the stratum j respectively, f max is the highest frequency in the variable. For example, assuming a variable was divided into three strata with the corresponding frequencies: f 1 = 200, f 2 = 160 and f 3 = 20. Thus, f 1 = 200 is the maximum frequency and f max = f 1 = 200. Then, the normalised priority for the strata was calculated following Equation (1): p 1 = 0, p 2 = 20 and p 3 = 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 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, p i was the normalised priority value of the variable i, and w i 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.

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.

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    (Grieve et al. 1999;Webb & Wilson 1995).
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 forest and mangrove area). Through this process, approximately 71.25% of the study area (i.e. 189,055.7 km 2 ), 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 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.

Range of interest
Min.     ; ).  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.

Variable
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 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 GISbased modelling (Carver 1991;Chen et al. 2010;Marinoni et al. 2009).