Multi Criteria Evaluation for Tourism

Study area

Jahan-Nama Protected area was taken as an investigation site. It is situated in the Northern piece of Iran and falls within Golestan Province with a territory of 31747 ha of which more than 15000 ha is occupied by rangelands and mixed forest-rangelands and 11000 ha by dense and semi-dense forest (Fig 1). Gravely soils are predominant in the area. Average yearly precipitation in the western and northern pieces on the investigation site is around 640 mm and 350 mm for eastern part. Elevation ranges between 590 m and 3033 m. Average slope of the area is about 21.5 percent. It enjoys a mid-climate during spring and summer. Natural interesting landscapes such as springs, rough mountains, forest routs; permanent and seasonal waterways and waterfalls make it as one of the main recreational destinations in the Golestan Province. Numerous valuable species of mammal and birds use this region as a habitat. Large number of travelers visits the area every year. In light of most-recent information recorded by Statistical Center of Iran (SCI) the number of inhabitants at the study site is 2819.

Figure 1: Map of the study area and observed points

Waterways, track roads in forest, waterfalls, springs, green landscapes, caves, histo-cultural heritages and other natural attractions are considered as the leisure destinations by travelers. Two destinations in Jahan-Nama protected area namely Ziarat and Jahan-Nama villages are the main tourist residential places.

MCE process

Multi criteria evaluation is a technique wherein multiple discrete and continues criterion could be combined to achieve a single composite basis for decision making (Eastman, 2012). This technique consists of series of consecutive stages included setting the goals, criteria identification, database establishment, calculation of criteria weights, standardization of maps and finally weighted aggregation of the layers. Figure 2 depicts succession stages of multi criteria evaluation.

Criteria identification proceeds according to past works within the field of tourism and additionally through expert knowledge. Some biophysical attributes determine the suitability of landscape for tourism, including, among others, climatic condition, topography of landscapes, flora and fauna, attractions and so on. Criteria which considered during present work were proposed during previous studies and are included: Visibility (Bunruamkaew and Murayam (2011);(Suryabhagavan et al., 2015)) Vegetation (Bali et al., 2015; Chhetri and Arrowsmith, 2008; Nino et al., 2017; Samanta and Baitalik, 2015), Biodiversity (Bunruamkaew and Murayam, 2011; Dhami et al., 2014), Elevation (Aklıbaşında and Bulut, 2014; Bunruamkaew and Murayam, 2011; Chhetri and Arrowsmith, 2008; Nahuelhual et al., 2013; Nino et al., 2017; Pareta, 2013; Samanta and Baitalik, 2015), Slope (Aklıbaşında and Bulut, 2014; Bunruamkaew and Murayam, 2011; Chhetri and Arrowsmith, 2008; Suryabhagavan et al., 2015), Distance to cultural heritages (Bunruamkaew and Murayam, 2011; Pareta, 2013; Suryabhagavan et al., 2015), Distance to roads ((Nahuelhual et al., 2013);(Samanta and Baitalik, 2015); (Aklıbaşında and Bulut, 2014); Distance to rural area ((Nino et al., 2017), Mahiny et al., 2009), Soil erosion (Bali et al., 2015), Temperature (Suryabhagavan et al., 2015; (Aklıbaşında and Bulut, 2014); Faults (Shojaee et al., 2012); Distance to water resources ((Chhetri and Arrowsmith, 2008); (Bunruamkaew and Murayam, 2011); Aspect (Mahiny et al., 2009; (Bali et al., 2015)), Distance to rivers (Mahiny et al., 2009), Precipitation (Aklıbaşında and Bulut, 2014); (Suryabhagavan et al., 2015), Scenic beauty ((Nahuelhual et al., 2013); Soil texture (Bali et al., 2015) and Number of sunny days ((Makhdoom, 2013). Moreover, Geology and distance to attractive points were also considered based on expert knowledge. All considered criteria and sub-criteria are shown in Table 1.

Visibility was acquired from the viewshed map which created from attractive points and digital elevation model. All cells visible to defined points over the surface were determined. Possible attractive points and tracks (for example, springs, rivers, and histo-cultural points) were imported as the source layer in viewshed. Visible cells were considered suitable and nonvisible ranked unsuitable in the final layer.

Greenery is one of the main tourist attractions all around the world. NDVI is one of the key indicators that implies on the density of vegetation. This index was calculated from Landsat 8 images band 4 (red) and 5 (near-infrared) as Eq. 1 (Carlson and Ripley, 1997).

NDVI= (NIR-Red)/ (NIR+Red)         (Eq.1)

Unlike green spaces, bare lands are not so much appealing for travellers; consequently, vegetation density lower than 0.05 ranked as unsuitable. Although vegetation can play a decisive role in attracting tourists, they often avoid areas with highly dense cover (Pérez-Maqueo et al., 2017); thereby density higher than 0.8 was also considered as unsuitable.

Presence of wildlife resources may increase the suitability of area for tourism activities. Wild animals are so interesting not only for visitors seeking natural habitats but also for those scientists coming to study some particular species (Dye and Shaw, 2007). Proximity to wildlife was recommended by experts to take into account fauna bio diversity in this study. In this way, distance operator was applied on wildlife hotspots map. A distance of 250 m was considered as buffer for natural habitats to prevent human impacts on the valuable species.

Topographies of landscapes along with flora and fauna are of the main resources in attracting tourists (Bahaire and Elliott-White, 1999). Slope, aspect and elevation are three primary components of topography, which consider in this study. These components were derived from Digital Elevation Model (DEM). Area with slight slopes is safer for tourists and the impact of tourist trampling on soil erosion may be diminished there (Yang et al., 2014). Slight slopes are also more feasible for all travellers from different ages. Steep slopes could increase the risk of soil erosion; hence slopes greater than 50 percent were classified as constraint. Aspect map was classified based upon expert knowledge thereby north – and west-faced slopes ranked as the most suitable followed by south- and east-faced slopes. Elevation was also classifies based upon expert knowledge and previous literatures. Due to possible physiological problems like respiratory disorders which occur in some visitors, altitude greater than 2500 m was considered as limited range.

Cultural heritages include known historical buildings, holy shrines and ancient hills were considered as the source point layer to calculate distance to cultural heritages. Maximum distance of 15 km regarded as the minimum suitability.

Accessibility refers to the ease with which a site or service may be reached or obtained (Nicholls, 2001). Accessibility of landscapes to visitors can affect the suitability of these areas for tourism activities so that the location which provides easy access to attractions may have the greater potential than sites with difficult access (Chhetri and Arrowsmith, 2008). The plans of tourists are also limited by accessibility of destinations (Ólafsdóttir and Runnström, 2009). To calculate accessibility of landscape, all roads (both asphalts and dirt roads) were extracted on Google’s earth and transferred to GIS to create the road layer. Distance operator was applied on road layer to estimate distance from the extracted roads. Distance greater than 5000 m from roads was considered as unsuitable.

Villages serve as service providers to visitors in many natural areas. Many of the histo-cultural destinations are located near the villages. On the one hand, the life-style of the indigenous people inhabited in rural areas is appealing for many of the travellers. On the other hand, tourism activities close to the villages might support local economy. Therefore, proximity to villages ought to incorporate in land allocation for tourism. Distance map was created around main villages located within the study area. Distance greater than 5000 m from villages thought-about as unsuitable based on expert knowledge and literature.

Where the soil erosion potential (erodibility) is in critical condition, tourism activities could result in enhanced soil loss. So taking into consideration soil erodibility in land evaluation for tourism is of a great importance. In the present work, soil erosion potential at the study site was predicted by RUSLE  following the Eq. 2 (Renard et al., 1997).

A=R.K.L.S.C.P               (Eq.2)

Where A is mean annual soil loss (ton/hectare/year),

R is rainfall and runoff erosivity factor (MJ/ha/mm/yr.),

K is soil erodibility factor (ton/MJ/mm),

L is slope length,

S is slope steepness,

C is Cropping-management factor,

P is erosion control factor

Soil erosion greater than 100 (ton/ha/year) was considered as unsuitable.

Climate plays important role in tourism (Day et al., 2013). Tourism demand is considerably influenced by the climatical factors, so that the selection of destinations strongly governed by weather conditions (Liu, 2016; Ridderstaat et al., 2014). Accordingly, temperature, rainfall and sunshine duration are three climatic parameters that incorporate during this study. Temperature is one of the explanatory factors in tourism, in order that the number of tourists in the future can be predicted by temperature (Liu, 2016). To create temperature map, monthly average was calculated over the tourism period (March to October). Temperature 37 C and 4 C was deemed as the thresholds for suitability. Sunshine duration, also, conjointly plays a decisive role in selecting destinations. Location within which the number of sunny days was lower than 7 days thought of unsuitable.

Tourism significantly depends on water resources, both in providing basic human needs and as an asset for a wide range of tourist activities (Gössling et al., 2012). Accessibility to high quality drinking and potable water is vital for tourism development (Koç et al., 2017).  Future Viability and sustainability of tourist destinations depends upon adequate water supply of sufficient quality and quantity (Essex et al., 2004). Therefore, considering the water resources in land allocation for tourism has a crucial importance. All water resources including springs, rivers and wells were choose as source point. A distance of 2000 m defined as the final suitability based upon the expert team. To protect the water resources in term of quality and surrounding ecosystem, a distance of 30 m from the source point was considered as the buffer.

Rivers are used for touristic purposes in many parts of the world (van Balen et al., 2014). Main perennial and ephemeral streams serve as scenic landscapes by a large number of visitors. Accordingly, proximity to rivers takes into account in land planning for tourism. A distance of 30 m was reserved as the buffer to ensure the protection of riparian ecosystems and river physical stability.

One of the most important attributes of the landscapes which determine the suitability of tourism is landscape scenery (Nahuelhual et al., 2013; Weyland and Laterra, 2014). Scenic beauty layer was  adopted from (Sakieh et al., 2017). This layer is a MCE output of biophysical attributes including spot with high topographic variability, vegetation density, vegetation ecotone, summit visibility, river visibility, vegetation type and diversity of vegetation.

To prevent landscape degradation in term of stability, soil texture and geology was classified primarily based on expert knowledge. Loamy texture has been considered as those with best suitability followed by loamy clay, clay loam, sandy clay loam, sandy and clay.

We also considered distance from main faults as a conservative sub criterion in this work. A minimum distance of 150 m from the main faults served as buffer.

Proximity to possible attractions including, among others, springs, rivers, forest roads and vegetation ecotons were also calculated using distance operator. As distance to attraction points increased, the suitability of landscape for tourism decreased.

Criteria weighting: Calculating the weight of each factor which affects the land suitability for tourism is a critical and essential step in land allocation. Different criteria have unequal degree of importance in land allocation. Analytical hierarchy process (AHP) as superior and effective method (Zhang et al., 2015) employed in MCE to weight multiple factors. AHP provides the weights of each factor through pairwise comparison of multiple factors based on matrix calculations. To implement AHP a two-layer structure was constructed in which in the first layer five criteria and the 20 sub-criteria in the second layer. AHP questionnaires fill out by 30 experts of tour operators, academic members and tourism students. Finally, questionnaires were imported into Expert Choice 11 to calculate the factors weight and consistency ratio. All criteria and their relative weights are shown in Table 1.

Standardization: As each factor has its own dimension with different value, using fuzzy sets the input layers are standardized in a 0 – 255 range. Variety of fuzzy membership functions employed including monotonically linear increasing, monotonically linear decreasing, linear symmetric and user defined. All constraints (the range in which there is no suitability for targeted land use) were prepared using Boolean logic; whereby all suitable values gives 1 and values exceeds the suitability thresholds gives 0. All fuzzy membership functions and thresholds are shown in Table 1.

Weighted Linear Combination (WLC): Weighted linear combination is an aggregation method whereby all layers should be standardized to a common numeric range (0-255). Then, fuzzified layer regarding to their relative weights, which are calculated by AHP, and constraints are integrated using Eq. 3.

SI=(∑i=1nWiXi)*Ci

(3)

Where:

SI is the suitability index,

W_i is the weight of factor i,

X_i is the fuzzified factor i,

C_i is the constraint i,

Table 1: Criteria, sub-criteria and their relative fuzzy membership functions, thresholds and weighting score

Logistic Regression (LR) process

Regresion analysis is an approach for estimating the relationships between a response variable and one or more independent variables. The Logistic regression is a non-linear statistical model for binary dependent variables (Shu et al., 2014) such as presence/absence point samples. It estimates the probability of ouccurance of the dependent variable as Eq. 4. This probability varies between 0 and 1.

P=11+e-y

(4)

Where P is the probability of the occurance of dependent variable, e is the basis of natural logarithm and y is the linearised regression function expressed as Eq. 5.

y=α+β1X1+β2X2+…+βnXn

(5)

Where α is the model intercept,

β₁, β_{2 …}β_n are the coefficeint of n independent variables X₁, X₂ … X_n.

LR employs Maximum Likelihood Estimation (MLE) to find the best fitting set of parameters given observation samples. LR employs two sampling approaches including stratified random sampling and systematic sampling. We adopt stratified random sampling as it provides reasonable balance between spatial dependency and population characteristics (Xie et al., 2005). Dependent variable should be in binary form (0, 1). In this study, point samples representing observed tourist destinations were used as dependent variable (Fig 1). Fig 2 depicts 10 biophysical factors which used as independent variables. The output of LR is a probability layer of dependent occurrence that indicates tourism potential, Pseudo R_square and Relative Operating Characteristic (ROC) (which measures the goodness of fit of logistic regression). R_Square and ROC 1 indicated the perfect fit while 0 denotes the lack of relationships between dependent and exploratory parameters. Pseudo R_square expressed as follow:

Pseudo RSquare=1- (logLikelihoodlogL0)

Where Log_likelihood is the value of the likelihood function for the full model as fitted and log_L0 is the value of likelihood function if all coefficients except the intercept are 0.

Fig 2: Independent layers used in logistic regression and multi-layer perceptron

Multi-Layer Perceptron (MLP) process

An ANN consists of several nodes and layers called input, hidden and output. One of the main advantages of artificial neural network approach is that it is independent of statistical distribution of data and there is no need for specific statistical variables (Pradhan and Lee, 2010) . MLP, a kind of feedforward networks, is one of the most frequently used artificial neural networks. MLP architecture consists of three layers including input, hidden and output layer and thus can recognize relationships that are non-linear in nature (Pijanowski et al., 2002). This type of network trained using Back Propagation (BP) learning algorithm. In BP two main steps namely forward values and backward propagation are applied to carry out its modification of the neural state. In the forward step, an input pattern presents to the input layers and its effects propagates through the network till an output produced. Produced output then compared with expected one and an error signal computed for each output nodes. In the back propagation, the error signal backwards to the all contributed nodes of output nodes. This stage continues till all nodes in the network receive an error signal. Then each node updates the value for each connection weights. In the back propagation the weights between input and hidden layer, and between hidden and output layer are calculated through modifying the number of hidden nodes and learning rate (Pradhan and Lee, 2010). Back propagation algorithm tries for a minimum value of global error through delta rule. Global error calculates as the average squared error between model output and targeted output.

One of the main steps in implementing MLP is to set the training and testing sample size. Training sub set of two small size may led to an ANN with limited or bad generalization (Basheer and Hajmeer, 2000). Too few samples may not represent the population for each category, while too many samples may cause samples to overlap, leading to a possible over training of the network (Eastman, 2012). There is no clear mathematical approach for determination of minimum size of training and testing sub set. Some studies proposed a ratio between training and testing sub sets. Nelson and Illingworth (1990), for example, propose 20-30% of database to be used for testing. In this study, 70% of data was considered for training and the rest for testing.

Network topology explains the number of input and output layer nodes and also the number of hidden layers. A three-layer neural back propagation network was used in this study (Fig. 3) in which input layer nodes = 10 (equals to input parameters), output layer nodes = 1, hidden layer = 1 and the hidden layers = 20.

Dynamic learning rate was used as training mode, during which start learning rate and end learning rate was set to 0.01 and 0.001, respectively; Momentum factor and sigmoid constant A was also set to 0.5 and 1, respectively. RMS = 0.01 and iterations = 10000 was set as stopping criteria. Finally, output function was set on sigmoidal mode. The number of input layer in MLP was same as LR inputs.

Fig 3: MLP Network used for tourism suitability assessment

Preparation of training and validation subsets for MLP and LR

Ground-based samples

Field observations were carried out from early spring to late summer 2016 during which 265 points were recorded as tourist destinations. Recorded points were transformed into QGIS 2.16 (Fig 1). 25% of samples were randomly selected and exported in a separate layer as a validation set. Validation samples were removed from the observation and the remaining 75% served as training set to develop the models.

MCE-based training samples

On the MCE output, pixels with the suitability values of 200 to 255 were recast as a new Boolean layer. Random sampling was implemented on reclassified layer to extract training samples.

Comparing accuraccy and landscape metrics of ground-based and MCE-based results

Accuraccy analysis

To perform the accuracy analysis of MCE-based and ground-based results, Receiving Operator Characteristic (ROC) method was employed. ROC analysis is useful specially, for cases in which the researcher wants to see how well the suitability map portrays the location of a particular category but does not have an approximation of the quantity of the category (Eastman, 2012). To perform ROC analysis, true positive were ploted against false positive of each thresholds. Thresholds refer to the precentage of pixel in model result for comparision with the refrence image (validation set). Area Under Curve (AUC) explain the ROC statistic. AUC calculates as follow:

AUC= ∑i=1nxi+1-xiyi+yi+1-yi2

Where:

AUC is area under curve,

X_i is the rate of false positive for threshold i,

Y_i is the true positive for threshold i,

n is the number of thresholds

AUC of 1 denotes the perfect performance of the model and AUC of 0.5 shows the random performance.

Calculating landscape metrics

Number of landscape metrics which explain heteroginity of landscapes in term of compactness and connectivity were calculated for patches with the high suitability for tourism development to compare MCE-based and ground-based results. This indices measures the basic features of landscape structure (area, size, clumpiness and the connectivity) (Sakieh et al., 2017). Hence, Zonal Land Suitability (ZLS) (Eastman, 2012) was performed for each output (MCE-based and ground-based outputs).

ZLS expressed as follow:

ZLSz= ∑Liznz

Where:

ZLS_z is the zonal land suitability of zone z,

(L_i)z is the local suitability of cells i belonging to the zone z,

n_z is the number of cells in zone z

Table 2 depicts selected landscape metrics, their description and their rages.

	CA	NP	MESH	SPLIT	ENN-MN
logistic	2472.57	6	45.1227	1345.987	2611.749
logistic_mce_based	2823.84	5	58.8828	1031.447	2406.743
mlp	2909.16	7	60.2256	1008.451	2979.351
mlp_mce_based	2902.05	5	65.5059	927.1613	3290.93

Aklıbaşında M, Bulut Y. Analysis of terrains suitable for tourism and recreation by using geographic information system (GIS). Environmental Monitoring and Assessment 2014; 186: 5711-5719.

Bahaire T, Elliott-White M. The Application of Geographical Information Systems (GIS) in Sustainable Tourism Planning: A Review. Journal of Sustainable Tourism 1999; 7: 159-174.

Bali ALI, Monavari SM, Riazi B, Khorasani N, Zarkesh MMK. A SPATIAL DECISION SUPPORT SYSTEM FOR ECOTOURISM DEVELOPMENT IN CASPIAN HYRCANIAN MIXED FORESTS ECOREGION. Boletim de Ciências Geodésicas 2015; 21: 340-353.

Basheer IA, Hajmeer M. Artificial neural networks: fundamentals, computing, design, and application. Journal of Microbiological Methods 2000; 43: 3-31.

Bunruamkaew K, Murayam Y. Site Suitability Evaluation for Ecotourism Using GIS & AHP: A Case Study of Surat Thani Province, Thailand. Procedia – Social and Behavioral Sciences 2011; 21: 269-278.

Carlson TN, Ripley DA. On the relation between NDVI, fractional vegetation cover, and leaf area index. Remote Sensing of Environment 1997; 62: 241-252.

Chhetri P, Arrowsmith C. GIS-based Modelling of Recreational Potential of Nature-Based Tourist Destinations. Tourism Geographies 2008; 10: 233-257.

Day J, Chin N, Sydnor S, Cherkauer K. Weather, climate, and tourism performance: A quantitative analysis. Tourism Management Perspectives 2013; 5: 51-56.

Dhami I, Deng J, Burns RC, Pierskalla C. Identifying and mapping forest-based ecotourism areas in West Virginia – Incorporating visitors’ preferences. Tourism Management 2014; 42: 165-176.

Dye AS, Shaw S-L. A GIS-based spatial decision support system for tourists of Great Smoky Mountains National Park. Journal of Retailing and Consumer Services 2007; 14: 269-278.

Eastman RJ. IDRISI Selva. Clark Labs, Clark University, Worcester MA, 2012.

Essex S, Kent M, Newnham R. Tourism Development in Mallorca: Is Water Supply a Constraint? Journal of Sustainable Tourism 2004; 12: 4-28.

Gössling S, Peeters P, Hall CM, Ceron J-P, Dubois G, Lehmann LV, et al. Tourism and water use: Supply, demand, and security. An international review. Tourism Management 2012; 33: 1-15.

Koç C, Bakış R, Bayazıt Y. A study on assessing the domestic water resources, demands and its quality in holiday region of Bodrum Peninsula, Turkey. Tourism Management 2017; 62: 10-19.

Liu T-M. The influence of climate change on tourism demand in Taiwan national parks. Tourism Management Perspectives 2016; 20: 269-275.

Makhdoom M. Foundation of Land Use Planning: University of Tehran 2013.

Nahuelhual L, Carmona A, Lozada P, Jaramillo A, Aguayo M. Mapping recreation and ecotourism as a cultural ecosystem service: An application at the local level in Southern Chile. Applied Geography 2013; 40: 71-82.

Nicholls S. Measuring the accessibility and equity of public parks: a case study using GIS. Managing Leisure 2001; 6: 201-219.

Nino K, Mamo Y, Mengesha G, Kibret KS. GIS based ecotourism potential assessment in Munessa Shashemene Concession Forest and its surrounding area, Ethiopia. Applied Geography 2017; 82: 48-58.

Ólafsdóttir R, Runnström MC. A GIS Approach to Evaluating Ecological Sensitivity for Tourism Development in Fragile Environments. A Case Study from SE Iceland. Scandinavian Journal of Hospitality and Tourism 2009; 9: 22-38.

Pareta K. REMOTE SENSING AND GIS BASED SITE SUITABILITY ANALYSIS FOR TOURISM DEVELOPMENT. nternational Journal of Advanced Research in Engineering and Applied Sciences 2013; 2: 43-58.

Pérez-Maqueo O, Martínez ML, Cóscatl Nahuacatl R. Is the protection of beach and dune vegetation compatible with tourism? Tourism Management 2017; 58: 175-183.

Pijanowski BC, Brown DG, Shellito BA, Manik GA. Using neural networks and GIS to forecast land use changes: a Land Transformation Model. Computers, Environment and Urban Systems 2002; 26: 553-575.

Pradhan B, Lee S. Landslide susceptibility assessment and factor effect analysis: backpropagation artificial neural networks and their comparison with frequency ratio and bivariate logistic regression modelling. Environmental Modelling & Software 2010; 25: 747-759.

Renard K, Foster G, Weesies G, McCool D, Yoder D. Predicting soil erosion by water: a guide to conservation planning with the Revised Universal Soil Loss Equation (RUSLE). 703. U.S. Departement of Agriculture,, 1997, pp. 404.

Ridderstaat J, Oduber M, Croes R, Nijkamp P, Martens P. Impacts of seasonal patterns of climate on recurrent fluctuations in tourism demand: Evidence from Aruba. Tourism Management 2014; 41: 245-256.

Sakieh Y, Salmanmahiny A, Mirkarimi SH, Saeidi S. Measuring the relationships between landscape aesthetics suitability and spatial patterns of urbanized lands: an informed modelling framework for developing urban growth scenarios. Geocarto International 2017; 32: 853-873.

Samanta S, Baitalik A. Potential Site Selection for Eco-Tourism : A Case Study of Four Blocks in Bankura District Using Remote Sensing and GIS Technology, West Bengal. International Journal of Advanced Research 2015; 3: 978-989.

Shojaee H, Kowalczyk A, Einafshar A. A tourism demand based method of geosites assessment on geotourism prioritization modeling:  The case of Razavi Khorasan Province Journal of Hospitality Management and Tourism 2012; 3: 82-94.

Shu B, Zhang H, Li Y, Qu Y, Chen L. Spatiotemporal variation analysis of driving forces of urban land spatial expansion using logistic regression: A case study of port towns in Taicang City, China. Habitat International 2014; 43: 181-190.

Suryabhagavan KV, Tamirat H, Balakrishinan M. Multi criteria evaluationn in identification of potential ecotourism sites in Hawassa town and its surroundings, Ethiopia. Journal of Geomatics 2015; 1: 86-92.

van Balen M, Dooms M, Haezendonck E. River tourism development: The case of the port of Brussels. Research in Transportation Business & Management 2014; 13: 71-79.

Weyland F, Laterra P. Recreation potential assessment at large spatial scales: A method based in the ecosystem services approach and landscape metrics. Ecological Indicators 2014; 39: 34-43.

Xie C, Huang B, Claramunt C, Chandramouli M. Spatial Logistic Regression and GIS to Model Rural-Urban Land Conversion Second International Colloquium on the Behavioral Foundations of Integrated Land-use and Transportation Models: Frameworks, Models and Applications. University of Toronto, Canada, 2005.

Yang M-y, Van Coillie F, Hens L, De Wulf R, Ou X-k, Zhang Z-m. Nature conservation versus scenic quality: A GIS approach towards optimized tourist tracks in a protected area of Northwest Yunnan, China. Journal of Mountain Science 2014; 11: 142-155.

Zhang J, Su Y, Wu J, Liang H. GIS based land suitability assessment for tobacco production using AHP and fuzzy set in Shandong province of China. Computers and Electronics in Agriculture 2015; 114: 202-211.