Yambaru forests, which are located in the northern part of Okinawa Island, have a precious ecosystem with many endemics and rare species in the flora and fauna. Public expectations of proper management for the preservation of this precious habitat are rising. Therefore, in September 2016 Yambaru National Park was established as the 33rd national park in Japan to protect this precious biodiversity. However, this area is also the center of forestry on Okinawa, covering around 60% of the private forest area in Okinawa Prefecture. Therefore, it is essential that proper forest management practice is developed so that the balance between environmental conservation and local development such as agriculture and forestry can be maintained.
It is essential for sustainable forest management that the distribution of forest resources is elucidated. Therefore, it is important for resource management that basic information becomes available so the forest productivity can be predicted.
Tree height is commonly used for various tree species as an indicator of forest productivity, for example in North America (McNab, 1993; Chen et al., 2002) and Europe (Corona et al., 1998; Seynave et al., 2005). Forest productivity can be expressed in terms of tree height in the so-called site index; generally, the site index is given by the mean tree height of the dominant trees at a reference age (for example, 40 years old). A high site index value indicates a favorable location where better growth is obtained in a short time compared with a low value. The site index is strongly influenced by soil and water conditions. The site index is estimated directly by measuring the tree height or forest volume, by indirect estimation using environmental factors such as soil, terrain, and climatic condition, or using topographic factors derived from a digital elevation model (DEM). Recently, DEM has been widely for obtaining basic information for sustainable forest management. DEM has been used for not only site index estimation but also the ecological land classification for landscape level using vegetation-environmental relations (Kusbach et al., 2017), modeling of climate conditions in forest vegetation zones (Machar et al., 2017), and ecosite-based site productivity evaluation (Pokharel & Dech, 2011).
In previous studies, many environmental factors have been adapted to predict site index. For instance; elevation, temperature, and rainfall have been used as regional climate factors; slope, aspect, and slope position have been used as topographic factors; and collected soil physical and chemical properties have been used as soil-related factors (Corona et al., 1998; Curt et al., 2001; Chen et al., 2002; Seynave et al., 2005). However, expert knowledge is necessary to obtain information such as soil type and surface geology, and more labor is required for field surveys. Recently, because GIS has become widespread and geographical data can be collected easily without labor-consuming on-site surveys, it has become easy to acquire large amounts of environmental factors from DEM. Thus, many studies have predicted the site index using only topographic factors (Teraoka et al., 1991; Mitsuda et al., 2001; Mitsuda et al., 2007).
Moreover, several studies have reported site indices for prediction and mapping based on a limited number of factors that can be easily derived from remote sensing data, such as Landsat (Günlü et al., 2008; Huang et al., 2017) and Light Detection And Ranging (LiDAR) (Laamrani et al., 2014). Spatially continuous maps are needed for tactical planning and operations for forest management; therefore, high-resolution site index maps are very powerful tools for local forest management (Huang et al., 2017).
In Japan, it has become possible to estimate the site index using topographic factors derived from DEM, especially for conifer plantations of trees, such as the Japanese larch (Larix kaempferi) (Mitsuda et al., 2001), Japanese cedar (Cryptomeria japonica) (Zushi, 2006; Mitsuda et al., 2007). In addition, these studies reported that topographic factors could better explain the site index in a relatively narrow area that had similar climatic conditions, such as temperature and precipitation. These studies are good examples that the site index can be estimated without including direct soil information by assuming that soil conditions correlate with topographic factors.
While research on conifer plantations has proceeded in many regions, including Japan, hardly any studies have been made of the relationship between forest growth and topographic factors in broad-leaved forests. Only a few studies have reported the relationships between site index and environmental factors (For example, Asato, 1979 ). These studies suggested that topographic factors related to soil moisture conditions are also closely related with the site index in subtropical broad-leaved forests. However, there are no studies that comprehensively elucidate the relationship between the site index and topographic factors with fine-scale resolution, and these relationships remain unknown.
If it were possible to predict forest productivity for natural broad-leaved forests using topographical information that can easily be derived from DEM, this would contribute enormously to regional forest management. In other words, highly accurate resource management can be expected even in regions, where it is difficult to conduct a resource inventory due to the shortage of labor and time and due to the expenses caused by the complexity of the topography and concomitant access difficulties. Besides, if it were possible to conduct efficient resource estimation in natural broad-leaved forests, where the demand for environmental conservation and forestry tend to compete, this would be effective as a conflict management tool, and provide practical and useful information to the regional forest management.
In order to contribute to regional forest management, this study predicted the forest productivity at a fine-scale resolution on Okinawa Island using topographical factors that are easily obtained in subtropical broad-leaved forests. The purposes of this study are 1) to clarify whether the topographic factors calculated from DEM and used to estimate the site index for conifer plantations or broad-leaved forests in previous studies, affect the site index in our study area, and 2) to develop a site index model using multiple regression analysis.
2. Materials and Methods
2.1. Study Area
The study area was a forest area of approximately 3500 ha located in Kunigami Village, in the northern part of Okinawa Island (approximately from 26˚46'40'' to 26˚49'57''N latitude and 128˚14'31'' to 128˚19'42''E longitude) (Figure 1). Most forests in the study area are subtropical evergreen broad-leaved forest, primarily Castanopsis sieboldii mixed with various tree species such as Schima wallichii. Ryukyu pine (Pinus luchuensis), which was planted in the 1960s and 1970s, occupies part of the study area. According to the forest register, the age of the forest stands is maximally 82 years, and stands of more than 50 years old occupy roughly 60% of the forest total. The altitude ranges from 3 to 420 m above sea level. The soil types are mainly dry to moderately moist yellow soils. Other types include red soils, surface gleyed red and yellow soils, dark red soils, gley soils, and immature soils (Ohnuki, 2017). Meteorologically, the climate is maritime subtropical; it is mild and with heavy rainfall. The mean annual air temperature and mean precipitation (from 1981 to 2010) at Oku weather station (AMeDAS) in Kunigami village is 20.7˚C and 2502 mm, respectively (data from the Japan Meteorological Agency website). Precipitation is the highest in June, and air temperature is the highest in August. The typhoon attacks from summer to autumn and the seasonal wind in winter cause damage to the agriculture and forestry (Kunigami village, 2009). From World War II to the present, most forests were affected by human activities such as logging and forest development (Azuma et al., 1997; Saito, 2011).
Figure 1. (a) Location of the study area on Okinawa Island; (b) Light gray indicates the forested area. Dark gray indicates the 38 - 42-year-old forests.
2.2. Calculation of Site Index
In our study area, airborne LiDAR data for forest management were taken in March 2010 (point density 1/m2, vertical accuracy 25 cm); digital orthophotography with 50 cm resolution (January 2010) and 25 cm resolution (March of the same year) was also performed. Datasets for a digital terrain model (DTM), a digital surface model, and a digital canopy height model (DCHM) were generated from the airborne LiDAR data. They are offered by Okinawa Prefecture as raster data with 1-m resolution.
Site index was calculated using relationships between tree height and forest age, as given by Terazono et al. (2015) and Terazono (2017). First, Terazono et al. (2015) conducted a smoothing process with a Gaussian filter to the digital ortho-images and generated the tree crown area of each individual dominant tree (Castanopsis sieboldii) using the watershed method. Then, they calculated the xy coordinates of each tree’s position by particle analysis. The height of each individual dominant tree was calculated by overlaying the tree’s position information with the DCHM image. Next, they created 10-m vector grid data and overlaid this data on the tree position map, which has information regarding individual tree height. Terazono (2017) calculated the mean dominant tree height within each 10-m grid and converted to 10-m grid raster data. Forest ages were determined from a forest register record (vector format) published by Okinawa Prefecture. Forest register record consists of forest compartment maps, which consist of forest management units and their attribute data from forest inventory for creating a forest management plan. The forest inventory data provides tree species and age for each forest compartment. The forest register record is revised using aerial photographs, forest cutting, and afforestation records according to government surveys conducted every 5 years. Terazono (2017) converted the forest compartment map with forest age to a 10-m grid raster to adjust the units of the mean dominant tree height map. He selected 179 sample grids that did not contain large reserved trees or large gaps referenced by the afforestation records and digital ortho-images and obtained values for forest age and tree height in those sample grids. An S-curve estimation model [Equation (1)] was created to estimate tree height from forest age using SPSS Statistics ver. 20. Finally, site index was calculated at a reference age of 40 years using the site index Equation (2) that was developed by Nishizawa & Mashita (1966). As a result, although there were errors in some places where the data of the reforestation register were incorrect, Terazono (2017) concluded that site index performed well in the remaining area. Therefore, we considered the site index value has enough accuracy to be used as an actual site index value.
where Yt (m) is the mean dominant tree height at t forest age; SI is the site index (m); Hi (m) is the observed mean dominant tree height of each grid in the study area; and x is the forest age of each grid in the study area.
2.3. Topographic Factors
Takeshita et al. (1966) stated that forest productivity could be represented by the quantity of water supplied to the soil if the density of nutrients in the soil water solution is constant, because forest productivity is determined mainly by the quantity of water-soluble nutrients forest trees can utilize. He then estimated forest productivity by using topographic factors closely related to soil moisture. In addition, he asserted that evaporation conditions (drying conditions), which are related to wind velocity, have a negative influence on stand growth. In this study, based on these hypotheses, we selected eight topographic factors related to soil moisture conditions and topographic exposure factors based on their known significant impact on site index described in previous works (Takeshita et al., 1966; Teraoka et al., 1991; Mitsuda et al., 2001; Chen et al., 2002; Zushi, 2006) and those that might be important in a sub-tropical forest in an island region. We calculated these factors using DEM (10 m cell resolution) resampled from a 1-m resolution to adjust the units of the site index grid by ArcGIS 10.3.1 and SAGA GIS 4.0.1.
Soil moisture factors
These factors are related to the amount of nutrient water that could be absorbed by trees for growth.
• Effective relief―the difference in elevation (m) between the highest and the specific cells within a horizontal distance of 100 m from the specific cell (Takeshita et al., 1966).
• Topographic wetness index―This factor is calculated based on the specific catchment area (α) and slope gradient (β). It characterizes the amount of the water accumulation depending on the slope’s position (Beven & Kirkby, 1979).
• Topographic position index―The difference between the elevation of a specific cell and the average elevation calculated for all cells in a specific radius (Guisan et al., 1999). Positive values mean the cell is higher than its surroundings while negative values mean it is lower, and values near zero are either flat or areas of constant slope (Jenness, 2006).
where Z0 is the elevation of a specific cell, Zn is the elevation of a cell within a specific radius, and n is the total number of surrounding cells in a specific radius (Seif, 2014). The specific radius was set at 50 m to examine the small-neighborhood slope position.
Topographic exposure factors
These factors are related to the local exposure that regulates wind velocity and which in turn determines evaporation from trees.
• The degree of exposure―The angle with the horizontal, in a direction which is not enclosed by surrounding mountains, over the inclination to the objective point (Takeshita et al., 1966). It is well known that there is a negative correlation between the degree of exposure and forest growth (Takeshita et al., 1966; Fukushima et al., 1974). We used an algorithm to calculate the degree of exposure from DEM developed by Murakami et al. (2000), with the cell size set at 100 m, the elevation angle set to 3 degrees, and the search radius set to 1 km.
• Openness (Positive openness)―This can be defined as the degree of enclosure of a location in the landscape and is an angular measure of the relation between surface relief and horizontal distance (Yokoyama et al., 2002). Values express the openness above the surface and are high for convex forms. The search radius was set at 1 km.
• Solar radiation―This factor is related to soil-moisture shortage. Global radiation (direct solar radiation and diffuse solar radiation) (Wh/ha) is calculated for each position of an input cell (ESRI, 2018). If the amount of solar radiation is large, the site is considered to dry out easily. On the other hand, if moisture conditions are sufficient, photosynthesis will proceed and transpiration will become active, so this index will also be an indicator of transpiration volume. We calculated the cumulative value of solar radiation during July and August, which have long sunshine hours in study area.
• Distance from the sea―This factor was used as an environmental factor to estimate the site index and accounts for the strength of seasonal winds characteristic of Okinawa Main Island (Asato, 1979).
• Elevation―This factor was adopted as a surrogate of the amount of precipitation and temperature. Generally, annual precipitation increases and temperature decreases as the elevation increases.
2.4. Statistical Methods
Firstly, the Pearson’s correlation coefficient was calculated to examine the relationships between site index and topographic factors. Subsequently, multiple linear regression analysis was used to examine these relationships. To select independent variables, a backward stepwise procedure was used, with the statistical significance level of entries and removals set at p = 0.05. Multicollinearity was determined using variance inflation factors (VIFs), with values larger than 5 indicating the existence of a serious multicollinearity problem in the model. All statistical analyses were conducted using SAS Add-In 6.1 for Microsoft Office (SAS, Cary, NC, USA). The topographic wetness index was log transformed and used for regression analysis.
For modeling, we randomly selected 340 grid cells from 38 - 42-year-old forests (approximately 1% of grid cells of these forests) using AcGIS 10.3.1 (Table 1). The sample size was sufficient for multiple regression analysis (Free Statistics Calculators, 2018). We randomly selected 340 other cells for validation. In forest management, the site index is given as the mean height of the dominant tree species at a reference age of 40 years. Usually, age classes with a width of 5 years
Table 1. Summary of site index and topographic variables of the 340 sample grid cells and 1st to 99th percentile of terrain attributes.
aSD, standard deviation.
are used, and thus, 38 - 42-year-old trees are classified as the same age class (40 years of age). We thus also used 38 - 42-year-old trees (approximately 340 ha in the study area) in our analysis. In addition, we used field survey data (n = 8) for validation. We created eight plots of 20 × 20 m2 for field measurements in the study area. On each plot, tree species, tree height, and diameter at breast height for all the trees measuring ≥ 1.2 m in height were measured. The location of each plot was identified using a global positioning system and topographic map (1:25,000 scale).
3.1. Correlation Coefficient
Site index was positively correlated with effective relief (r = 0.749, p < 0.001) and topographic wetness index (r = 0.396, p < 0.001), and correlated negatively with openness (r = −0.729, p < 0.001), topographic position index (r = −0.689, p < 0.001), solar radiation (r = −0.473, p < 0.001), elevation (r = −0.307, p < 0.001) and degree of exposure (r = −0.289, p < 0.001) (Table 2). There was a correlation between most topographic factors, especially openness had a strong positive correlation with topographic position index, and a strong negative one with effective relief. In addition, elevation had a strong positive correlation with distance from the sea (r = 0.809, p < 0.001).
3.2. Multiple Regression Analysis Model
As the independent variables for estimating the site index, the effective relief, openness, and elevation were selected in the model [Equation (3)] (Table 3). The model was significant at the 0.1% level. The adjusted R2 was 0.677.
Table 2. Pearson’s correlation coefficient between site index and topographic factors (n = 340).
*p < 0.05; **p < 0.01; ***p < 0.001.
Table 3. Multiple regression model for predicting site index from topographic factors (n = 340).
a. Effective relief, b. Openness, c. Elevation.
Focusing on the standardized partial regression coefficient, all the variables adopted in the model equation were statistically significant (Table 4). Also, when the influence of each independent variable on the site index was evaluated based on the t-value of the partial regression coefficient, the effective relief (t = 10.670, p < 0.001) was shown to have a relatively large impact on the model. The remaining two variables have almost the same degree of influence. The VIF of all variables in the model ranged from 1.026 to 2.189, and did not exceed 5, indicating that there was no serious multicollinearity.
We estimated the site index using the validation data (Figure 2). The root mean square error (RMSE) was 2.3 m. The predictive spatial site index map is shown in Figure 3 using the multiple regression model. Effective relief, openness, and elevation raster (created for the entire study area) were used.
We created a multiple regression model to explain the site index using topographic factors (effective relief, openness, and elevation) for broad-leaved forests in northern Okinawa Main Island. As a result, we were able to explain about 68% of the variance, which compares favorably with other values in the literature for Norway spruce (Picea abies (L.) Karst.) (Seynave et al., 2005) and Douglas-Fir (Pseudotsuga menziesii) (Mirb) Franco) (Corona et al., 1998; Curt et al., 2001).
Figure 2. Comparison of actual and estimated site indices using validation data (n = 348). Actual site index indicates the value calculated from Equation (2). Estimated site index indicates the value predicted using multiple regression models [Equation (3)]. The dashed line represents perfect correlation. Triangle symbols indicate the estimated value using field measurement data (n = 8).
Figure 3. Map of the estimated site index of the entire study area based on prediction from the multiple regression model. The site index map is shown overlaid on the shading map to facilitate recognition of the terrain attributes.
Table 4. Standardized partial correlation coefficient variables.
a. Variance inflation factor.
Our results indicate that effective relief is the most important factor of the site index. The performance of effective relief has also been reported in previous studies. For example, the effective relief greatly affected the site index in Japanese cedar (Takeshita et al., 1966; Teraoka et al., 1991) and Japanese cypress (Chamaecyparis obtusa) (Teraoka et al., 1991).
Openness has also been reported to be an important factor by Zushi (2006). The degree of exposure, which is conceptually similar to openness, correlated negatively with site index in our study area. This agrees with other studies for Japanese cedar (Takeshita et al., 1966), Japanese cypress (Fukushima et al., 1974), and Japanese larch (Mitsuda et al., 2001).
Elevation has also been identified as an important factor for the prediction of the site index for trembling aspen (Populus tremuloides Michx.) (Chen et al., 1998; Chen et al., 2002) and Norway spruce (Seynave et al., 2005). We used this index as a surrogate for climatic variables, and it might be useful for the prediction of site indices in regions where it is difficult to obtain fine-scale climatic data, because it can easily be deduced from DEM.
Only few reports on forest productivity estimation for the subtropical broad-leaved forest have appeared until now. Therefore, this research is important for forest management in this region. Asato (1979) analyzed the relationship between the site index at 30 years and eleven environmental factors such as climate, topography, geology, and soil using quantification method I. This is a method for predicting the quantitative external criterion or criterion variable on the basis of the information concerning the qualitative attributes of each subject and for analyzing the influence of each attribute to the criterion variable (Tanaka, 1979). In his analysis, the multiple correlation coefficient was 0.808 (the coefficient of determination was 0.653), and he pointed out that this was sufficient to be of use for estimating the site index in that region. Compared with their results, since our model has a coefficient of determination of 0.680, it can be said that it has almost the same accuracy as the results of Asato (1979).
The model created by Asato (1979) was intended to be used for predicting average forest stand growth as basic information for regional forest planning. His study area covered the broader northern forest area, including our study area. Currently, estimation on a spatially fine-scale is required because large-scale forest logging is not conducted, and the unit of forest practice and the size of forest loggings are becoming smaller. Our model has a resolution surpassing the scale of forest stand units. Therefore it may contribute essential information for current forest management such as the selection of small-scale logging sites. Furthermore, by incorporating other factors such as soil condition, sufficient improvements can be expected in the estimation.
Approximately 30% of the variance in the site index remained unexplained. Other factors that have been reported to contribute to site index variance are climatic factors, former land use, genetics, and soil properties (Corona et al., 1998; Curt et al., 2001; Fontes et al., 2003; Laamrani et al., 2014). In particular, several studies have identified the importance of soil properties on forest productivity. For example, organic layer thickness was found to be the most important variable for the estimation of the black spruce (Picea mariana (Mill.)) site index (Laamrani et al., 2014). A soil property model including soil physical and chemical properties was found to be a better predictor for trembling aspen than climatic, topographic, and foliar nutrient models (Chen et al., 1998). Ohnuki (2017) investigated soil type and surface soil layer thickness at the study site. He reported that the soil type shifted from the ridge toward the valley and that the thickness of the surface soil layer changed depending on the microtopography. Thus, the topographic factors used in our analysis were considered insufficient to express the soil’s condition in island areas. Only soil maps with coarse resolution exist in this region now. However, if detailed soil maps are developed, the accuracy will be expected to be improved by incorporating such information into the analysis.
On Okinawa Island, it has traditionally been believed that winds strongly influence forest growth due to seasonal winds and strong typhoons. For example, gentle-to-moderate slopes (Purves et al., 2009; Zawawi et al., 2014), where the surrounding mountains serve to protect and enclose the area from wind exposure, exhibit high forest growth. Conversely, flat terrain, gentle slopes (Purves et al., 2009; Zawawi et al., 2014), and mountain ridges (Iguchi et al., 2008) are considered to be at the lowest grade for forest growth. In areas where wind exposure is high, there is no barrier for protection; therefore, the area has a high potential for damage from winds. Although some areas along the coastline had higher estimated values for site index due to the low value for elevation and high value for effective relief, overall traditional empirical knowledge of forest growth was well reflected in our study results (Figure 3). Site indices have high values in mountainous valleys and slope depressions, whereas values are low on mountain ridges or relatively flat areas. Thus, traditional knowledge of forest growth characteristics in Okinawa Island was well reflected by these study results.
Using LiDAR-derived DEM, we were able to estimate the site index in a wide area at a fine-scale. Even in places where acquisition of field measurement data is limited, the possibility of estimating the site index easily in a wide area using DEM was shown. We only confirmed relationships between LiDAR derived site index and filed measurement data at the limited number of field data. It may be possible to improve the estimation accuracy and be a substitute for the field survey in the future by further clarifying the relationships between the LiDAR-derived data and filed measurement.
The prediction of forest productivity is essential for sustainable forest management, particularly in areas where environmental conservation is also important. Our study predicts the site index at a fine-scale resolution using multiple regression analysis with easily calculated topographic factors for broad-leaved forests in northern Okinawa Main Island. Consequently, three topographic factors were selected as independent variables. Sixty-eight percent of the variance could be explained with the effective relief, openness, and elevation. Traditional empirical knowledge of forest growth related to strong winds was also well reflected in the estimated site index map. It was thus possible to predict forest productivity with resolution at a finer scale than ever before. Our model has a resolution surpassing the scale of forest stand units and may therefore contribute essential information for current forest management, such as the selection of small-scale logging sites. Furthermore, by incorporating other factors, such as soil condition, and using other modeling techniques, the accuracy of the model may be further improved. We expect that the results of our research can be used to improve future conservation and management plans.
A.M. and A.S. conceived and designed the research; A.M. and T.R. analyzed the data; S.M. helped on statistical analysis. A.M. wrote the manuscript with assistance from all co-authors. A.S. and A.M. performed the field measurement. Special thanks to Dr. Naoki Kabeya, Dr. Takanori Shimizu, and Takuya Arakaki for providing GPS data and field measurement data. This study was supported by the “Research Project with regard to Environmental Conservation and Resource Utilization of Forest Area in Nansei Islands (2012-2016),” which was funded by Okinawa Prefecture in Japan.
 Asato, I. (1979). Anettaiseitennennkouyoujurinbunn no segyokaizennnikannsurukenkyu (VII)-kankyoinshi to jyukoseityonituite (Study on Improving the Operation of Subtropical Broad-Leaved Forest (VII)—Relationships between Environmental Factors and Tree Height Growth). Transactions of the 32nd Annual Meeting of the Kyusyu Branch of the Japanese Forestry Society, 32, 39-40. (In Japanese)
 Azuma, S., Sasaki, T., & Ito, Y. (1997). Effects of Undergrowth Removal on the Species Diversity of Insects in Natural Forests of Okinawa Honto. Pacific Conservation Biology, 3, 156-160. https://doi.org/10.1071/PC970156
 Chen, H. Y. H., Klinla, K., & Kabzems, R. D. (1998). Site Index, Site Quality, and Foliar Nutrients of Trembling Aspen: Relationship and Predictions. Canadian Journal of Forest Research, 28, 1743-1755. https://doi.org/10.1139/x98-154
 Chen, H. Y. H., Krestov, P. V., & Klinka, K. (2002). Trembling Aspen Site Index in Relation to Environmental Measures of Site Quality at Two Spatial Scales. Canadian Journal of Forest Research, 32, 112-119. https://doi.org/10.1139/x01-179
 Corona, P., Scotti, R., & Tarchiani, N. (1998). Relationship between Environmental Factors and Site Index in Douglas-Fir Plantations in Central Italy. Forest Ecology and Management, 110, 195-207. https://doi.org/10.1016/S0378-1127(98)00281-3
 Curt, T., Bouchaud, M., & Agrech, G. (2001). Predicting Site Index of Douglas-Fir Plantations from Ecological Variables in the Massif Central Area of France. Forest Ecology and Management, 149, 61-74. https://doi.org/10.1016/S0378-1127(00)00545-4
 Fontes, L., Tomé, M., Thompson, F., Yeomans, A., Luis, J. S., & Savill, P. (2003). Modelling the Douglas-Fir (Pseudotsuga menziesii (Mirb.) Franco) Site Index from Site Factors in Portugal. Forestry, 76, 491-507.
 Fukushima, T., Takai, J., Takeshita, K., & Takata, M. (1974). Stand Analysis Concerning Height-Estimation of Forest-Land Productivity (for Chamaecypris obtuse Siebetzucc.). The Bulletin of Fukuoka-ken Forest Experiment Station, 23, 1-34.
 Günlü, A., Baskent, E. Z., Kadiogullari, A. I., & Ercanli, I. (2008). Classifying Oriental Beech (Fagus orientalis Lipsky.) Forest Sites Using Direct, Indirect and Remote Sensing Methods: A Case Study from Turkey. Sensors, 8, 2526-2540.
 Huang, S., Ramirez, C., Conway, S., Kennedy, K., Kohler, T., & Liuc, J. (2017). Mapping Site Index and Volume Increment from Forest Inventory, Landsat, and Ecological Variables in Tahoe National Forest, California, USA. Canadian Journal of Forest Research, 47, 113-124. https://doi.org/10.1139/cjfr-2016-0209
 Iguchi, T., Takashima, A., Yoshida, S., & Mizoue, N. (2008). Topographic Effects on the Forest Structure in Subtropical Secondary Forests on Northern Part of Okinawa Island. Kyusyu Journal of Forest Research, 62, 185-187.
 Kusbach, A., Friedl, M., Zouhar, V., Mikita, T., &Šebesta, J. (2017). Assessing Forest Classification in a Landscape-Level Framework: An Example from Central European Forests. Forests, 8, 461. https://doi.org/10.3390/f8120461
 Laamrani, A., Valeria, O., Bergeron, Y., Fenton, N., Cheng, L. Z., & Anyomi, K. (2014). Effects of Topography and Thickness of Organic Layer on Productivity of Black Spruce Boreal Forests of the Canadian Clay Belt Region. Forest Ecology and Management, 330, 144-157.
 Machar, I., Vlckova, V., Bucek, A., Vozenilek, V., Salek, L., & Jerabkova, L. (2017). Modelling of Climate Conditions in Forest Vegetation Zones as a Support Tool for Forest Management Strategy in European Beech Dominated Forests. Forests, 8, 82.
 Mitsuda, Y., Yoshida, S., & Imada, M. (2001). Use of GIS-Derived Environmental Factors in Predicting Site Indices in Japanese Larch Plantations in Hokkaido. Journal of Forest Research, 6, 87-93. https://doi.org/10.1007/BF02762493
 Mitsuda, Y., Ito, S., & Sakamoto, S. (2007). Predicting the Site Index of Sugi Plantation from GIS-Derived Environmental Factors in Miyazaki Prefecture. Journal of Forest Research, 12, 177-186. https://doi.org/10.1007/s10310-007-0004-1
 Murakami, T., Teraoka, Y., & Imada, Y. (2000). Calculation of Topographic Factors and Their Accuracy Evaluations with Digital Maps 50-M Grid Published by the Geographical Survey Institute. Japanese Journal of Forest Planning, 34, 12-26.
 Ohnuki, Y. (2017). Elucidation of the Influence of Forest Logging on Geographical Condition and Hydrological Climate Environment—Impact on Forest Environment. In Research Report on Environmental Conservation and Resource Utilization of Forest Area in Nansei Island (pp. 1-24). Okinawa: Okinawa Prefecture Forest Resources Research Centre.
 Purves, J. M., Chen, B., & Nakama, Y. (2009). The Secret of Forestry: An English Translation of the Sanrin Shinpi of Sai on. Scientific Bulletin of the Faculty of Agriculture, University of the Ryukyus, 56, 11-21.
 Seynave, I., Gégout, J. C., Hervé, J. C., Dhôte, J. F., Drapier, J., Bruno, é., & Dumé, G. (2005). Picea abies Site Index Prediction by Environmental Factors and Understorey Vegetation: A Two-Scale Approach Based on Survey Databases. Canadian Journal of Forest Research, 35, 1669-1678. https://doi.org/10.1139/x05-088
 Takeshita, K., Fukushima, T., Hagiwara, Y., & Saijoh, T. (1966). Stand Analysis Concerning Estimation of Forest-Land Productivity. The Bulletin of Fukuoka-ken Forest Experiment Station, 18, 41-76.
 Teraoka, Y., Masutani, T., & Imada, M. (1991). Estimating Site Index of Sugi and Hinoki from Topographical Factors on Maps for Forest Management. Science Bulletin of the Faculty of Agriculture, Kyushu University, 45, 125-133.
 Terazono, R., Arakaki, T., Washizaki, K., & Ikuzawa, H. (2015). Estimation of Forest Growing Stock and Mapping of Forest Logging Area in Resource Recycling Forests. In Research Report on Environmental Conservation and Resource Utilization of Forest Area in Nansei Island (pp. 109-119). Okinawa: Okinawa Prefecture.
 Terazono, R. (2017). Estimation of Forest Growing Stock with High-Resolution Data in Resource Recycling Forest. In Research Report on Environmental Conservation and Resource Utilization of Forest Area in Nansei Island (pp. 136-150). Okinawa: Okinawa Prefecture Forest Resources Research Centre.
 Yokoyama, R., Shirasawa, M., & Pike, R. J. (2002). Visualizing Topography by Openness: A New Application of Image Processing to Digital Elevation Models. Photogrammetric Engineering and Remote Sensing, 68, 257-265.
 Zawawi, A. A., Shiba, M., & Jemali, N. J. N. (2014). Landform Classification for Site Evaluation and Forest Planning: Integration between Scientific Approach and Traditional Concept. Sains Malaysiana, 43, 349-358.
 Zushi, K. (2006). Spatial Distribution of Soil Carbon and Nitrogen Storage and Forest Productivity in a Watershed Planted to Japanese Cedar (Cryptomeria japonica D. Don). Journal of Forest Research, 11, 351-358.