1 Introduction

In many developing countries, beekeeping activities have an importance to rural economic development through the derived products such as honey, pollen, beeswax, royal jelly, bee venom, and propolis which are very important for human health and crop pollination (Estoque and Murayama 2010, 2011; Ceylan 2004; Damián 2016). Because honeybees are the key pollinator of 33% of crop species, there is a high amount of indirect economic income that involved in agricultural activities (Maris et al. 2008; Oldroyd and Nanork 2009). Turkey has considerable potential in beekeeping with her rich flora, proper ecological conditions and existence of colony. According to the 2015 beekeeping statistics, Turkey has a rapidly increasing honey production with 107.665 tons and 7.709.636 beehives. However, Turkish beekeeping sector has not utilized the rich natural resources sufficiently. Thus, management and monitoring beekeeping activities are being more important to provide efficient and sustainable productivity. Furthermore, determining suitable locations for beekeeping should be evaluated in the field of land use planning considering economical, ecological, environmental and social aspects within spatiotemporal perspective (Awad et al. 2019).

Land suitability analysis (LSA) can be assessed on the basis of physical environmental, social, and economic data (FAO 1976; Jafari and Zaredar 2010; Zhang et al. 2015). Land use should be planned to meet human needs and ensure the sustainability of ecosystems (Amiri and Shariff 2012) and optimum use of the resources for sustainable land management by identifying the most appropriate future land planning according to the requirements and preferences (Ahamed et al. 2000; Collins et al. 2001; Malczewski 2004; Zolekar and Bhagat 2015).

Multicriteria decision analysis (MCDA) techniques are widely used for LSA. MCDA of land suitability involves multiple criteria like elevation, slope, atmospheric conditions and land use, etc. as well as environmental and socio-economic approaches to find best solutions within multiple alternatives (Wang et al. 1990; Joerin et al. 2001; Yu et al. 2011; Zolekar and Bhagat 2015). One of the most applied MCDA approaches is the analytical hierarchy process (AHP) which calculates the weights of criteria among the factors that affect the total suitability (Saaty 1977, 1980, 1994, 2001; Saaty and Vargas 1991). AHP refers to the applications which are used to determine the most suitable solutions to the real problems by providing a selection of different data clusters (Arentze and Timmermans 2000) and calculates the weights associated with criteria via pairwise preference matrix where all criteria are compared against each other (Chen et al. 2010). The calculated weights represent the importance of criteria relatively which will contribute to the generation of suitability map. The AHP is widely used MCDA method compared to others such as the technique for order of preference by similarity to ideal solution (TOPSIS) and vise kriterijumska optimizacija I kompromisno resenje (VIKOR) due to the high applicability rate of AHP to a wide range of disciplines for spatial and non-spatial data. TOPSIS is based on determining the best alternative which has the shortest distance to positive ideal solution and longest distance from negative ideal solution (Hwang and Yoon 1981). The positive ideal solution represents the maximized benefit criteria and minimized cost criteria. In other words, the negative ideal solution represents the maximized cost criteria and minimized benefit criteria (Sakthivel et al. 2015; Wang and Elhag 2006; Ho et al. 2010). Finally, the VIKOR method is a MCDA method to determine the compromise ranking and compromise solution via given criteria weights. The VIKOR focuses on ranking and selecting from a set of alternatives by multi-criteria ranking index based on a measuring the distances to the ideal solution. The compromise ranking list can be determined by calculating the closeness of alternatives to the ideal solution (Opricovic 1998).

In literature, there are quite a few studies focused on beekeeping suitability via MCDA. (Abou-Shaara et al. 2013), used AHP to determine beekeeping suitability considering maximum temperature, relative humidity, summer crop area, water resources, and land cover criteria (Maris et al. 2008), determined suitability via AHP including rainfall, topographic, hydrology, road network, nectar and pollen classes criteria. (Amiri and Shariff 2012), used geographical information system (GIS) to determine suitability based on road and water availability, temperature, and precipitation rating, land use, and vegetation. Similarly, Estoque and Murayama (2010) considered distance to water and roads, elevation, nectar, and pollen class criteria within AHP. Camargo et al. (2014) detected the suitability by calculating 3 km buffer for each located beekeeper and evaluated land use, flora, and honey productivity. Land use, flora, solar radiation, distance to water resources, distance to electromagnetic radiation, climatic conditions, distance to urban areas and road networks criteria are considered within AHP by (Fernandez et al. 2016). Nor (2007), used AHP method to generate beekeeping suitability with weighted overlay function for Malaysia. All the studies that focused on locating and evaluating suitable areas involve similar criteria and characteristic features of beekeeping requirements. Although considered criteria are similar to recent studies, in this study, TOPSIS and VIKOR methods are applied in addition to AHP to make a comparison of methods for beekeeping suitability at Konya in Turkey. The main aim of this study is not only improve beekeeping in Konya province but also generating a conceptual model for beekeeping suitability assessment which can be applied to any region in the world. The diverse topographic, climatic, and environmental characteristics of Konya can simulate all regions in Turkey, and because of this, the study area and application will establish a conceptual model for beekeeping suitability assessment for all around Turkey and other countries in the world by proving the reliability and applicability of MCDA models. Slope, elevation, aspect, flora, precipitation and distance to water resources, roads and settlements criteria are considered. The suitability maps for beekeeping were generated with AHP, TOPSIS, and VIKOR methods and GIS integration.

2 Material and method

2.1 Study area

The study was done at Konya and it is the largest city (40,814 km2) in Turkey with its 31 districts (Fig. 1). Konya has the largest agricultural lands in Turkey with 19,239,667 decare according to the 2015 statistics (URL 1). Due to the large area (bigger than some countries in the world), climatic conditions, flora and topographic features are being different and could be suitable or not to beekeeping.

Figure 1.
figure 1

Study area Konya map and boundaries.

Full size image

Konya is one of the main centers of grain farming with its plain and large agricultural lands. According to the Turkish Statistical Institute, Konya has 18,618,142 decare for grain farming, 207,665 decare for vegetable gardens, and 412,918 for fruits, beverage, and spices within total 19,239,667 decare arable lands. The forests are found mostly in the mountainous parts of the province and consist of black pine, oak, red pine, juniper, cedar, and fir, respectively. The topography of Konya is mostly plain in the middle of the city and has high mountains (Toros Mountains) in the south of the city that compose the boundary with the Mediterranean region and its climate. The city has average 1020 height above sea level and has water resources with 2127 m2. The city has cold semi-arid climates which typically found in continental interiors and have distance from large bodies of water. Due to the different characteristic features of Konya topography, it has rich flora that suitable for beekeeping.

2.2 Methodology

The AHP, TOPSIS, and VIKOR methods were applied to determine the beekeeping suitability for study area. The methods and calculations are explained step by step and given in Fig. 2.

Figure 2.
figure 2

Implementation model of beekeeping suitability analysis.

Full size image

2.2.1 Criteria selection

The criteria selection reflects the requirements, expectations, and restrictions of beekeeping activities when locating beehives in the field of topographic, environmental, meteorological, and economical perspective. Advanced beekeeping activities require being in ideal intervals for each criterion which are defined by bee experts. The selected criteria, scale/resolutions, and data sources are given in Table I.

Table I Spatial data resolutions, scales and related institutes
Full size table

Aspect

Aspect criterion is included to be able to determine the direction effect (Fig. 3a). Considering beehives locations and directions, beekeepers prefer south, south-east, and south-west directions when locating beehives to benefit from the daylight. These directions are also important to protect them from north winds. The aspect map is derived from ASTER GDEM data downloaded from web site at 30 × 30 m resolution (URL 2).

Figure 3.
figure 3

a Aspect map for directions. b Elevation map. c Flora map for land use. d Map for distances to roads. e Map for distances to water resources. f Map for distances to settlements. g Slope map to define topography. h Annual precipitation map of study area.

Full size image

Elevation

Elevation criterion related to flora and defines seasonal start of the beekeeping activities. For study area, honey production yield and efficiency is decreasing above 2000 m due to the meteorological conditions and winds. The elevation is varying from 591 to 3419 in Konya city boundaries (Fig. 3b). Elevation maps are downloaded from ASTER GDEM web site (URL 2) at 30 × 30 m resolution.

Flora

Flora of the study area defines the honey production quality and quantity addition to honey type (Amiri and Shariff 2012; Abou-Shaara et al. 2013). Thus, the most important criterion should be flora and weighted higher values than others (Fig. 3c). Forests and natural plant areas are preferred to benefit from natural plant diversity. Urban settlements and industrial areas are not included to avoid disadvantages and effects of urbanization on honey production. Although agricultural lands have an important role on honey production, pesticide using is one of the main risks for bees. Thus, agricultural lands are weighted as non-important.

Distance to roads and settlements

Beekeepers prefer to locate beehives outside of urban places and roads to decrease greenhouse gases, air and noise pollution, exhaust emissions, urban and industrial contaminants, and human-related factors (Maris et al. 2008). Thus, distance from settlements and distance from highways criteria are included in suitability analysis (Fig. 3d, f).

Settlement and road data are digitized from Google Earth 2016 imageries. The distances are calculated by using buffer analysis in ArcGIS software and converted to a raster image.

Distance to waters

Water resources are important for bees to provide enough water that will be used for cooling the colonies and honey production (Amiri et al. 2011). The city has average 1020 elevation above sea level and has water resources with 2127 m2 (Fig. 3e).

Slope

Similar to elevation, slope criterion has a close relationship with flora due to rapidly changing topography, meteorological conditions, and directions.

Slope map is generated by using ArcGIS software with using ASTER GDEM data which available at (URL 2). The ASTER GDEM data have elevation value for each pixel at 30 × 30 m resolution. Slope map is derived from elevation data by calculating the elevation difference between each pixel and vary from 0 to 71.2% (Fig. 3g).

Precipitation

Precipitation has a close relationship with flora and defines the characteristic features of study area. Precipitation expected to be between 1275 and 1800 mm annual rainfall (Maris et al. 2008) and related with elevation, flora, and its flowering season (Fig. 3h).

Each criterion is mapped and then reclassified with the ArcGIS software according to the defined classes. In each figure, the suitability value is illustrated from highly suitable (green) to none suitable (red) relatively.

2.2.2 AHP

The AHP, proposed by Saaty (1980), is a complex decision-making tool that considers the priorities of each criterion. For this purpose, AHP establishes an importance scale from 1 to 9 (1 = equal, 3 = moderately, 5 = strongly, 7 = very, 9 = extremely). Moreover, the AHP provides a consistency rate concept to be able to calculate the consistency of overall weights and priorities. The AHP method provides a consistency ratio which should be less than 0.1 to prove that the weights and priorities are consistent.

2.2.3 TOPSIS

The TOPSIS technique refers to a distance calculation which assigns shortest distance from positive ideal solution and longest distance from negative ideal solution as best alternative. These distances are included in a similarity index concept which will be ranked to find best solutions. The ranking values are used to calculate R and V matrices that represent normalized decision matrices via the W criterion weights that are calculated with the AHP (Hwang and Yoon 1981). The number of positive ideal solutions (Di*) and negative ideal solutions (Di−) will be equal to the number of alternatives (Peters and Zelewski 2007; Triantaphyllou 2000). Finally, best solutions are determined by relative closeness to the ideal solution Ci* (1 > Ci* > 0).

2.2.4 VIKOR

Within the VIKOR method L1,j = Sj and L∞j = Rj are used to model a measure for ranking. The minimum Sj refers to the maximum group utility and the minimum Rj refers to the individual regret of the opponent (Yu 1973; Zeleny 1982; Opricovic 1998). The VIKOR method calculation requires to determine the best \( {f}_i^{\ast } \) and \( {f}_i^{-} \) values of all criteria functions. These values are then used to calculate the S and R values to be able to calculate the Q values. The minimum Q value is assigned the highest suitable and the maximum Q value is none suitable.

3 Beekeeping suitability application

Generating suitability maps require calculating the weights of each criterion to determine the importance of criteria to each other. AHP pairwise matrix is used to calculate the weights of criteria by using ranking values from 1 to 9.

In the first stage, criteria weights are calculated with a pairwise matrix via AHP by specifying the importance of each criterion to another. A 0.081 consistency ratio value means the weights are consistent. The calculated weights and pairwise comparison matrix is given in Table II.

Table II AHP beekeeping suitability pairwise matrix
Full size table

Total AHP suitability is calculated by multiplying the criteria maps with the following formula:

$$ \boldsymbol{TS}={\sum}_{i=1}^n{w}_i.{r}_i={W}_{AS}. AS+{W}_{EL}. EL+{W}_{SL}. SL+{W}_{FL}. FL+{W}_{DtR}. DtR+{W}_{PR}. PR+{W}_{DtW}. DtW+{W}_{DtS}. DtS $$

The values of the calculation will be between 1 and 9. According to the FAO (1976), land suitability classes are divided into 5 classes as highly suitable (S1), moderately suitable (S2), marginally suitable (S3), currently not suitable (N1), and permanently not suitable (N2) to be able to classify land suitability from none suitable to highly suitable. Then, suitability values are classified considering AHP suitability rates as follows:

  • 7 < TS < 9: highly suitable (S1)

  • 6 < TS < 7: moderately suitable (S2)

  • 5 < TS < 6: slightly suitable (S3)

  • 4 < TS < 5: moderately not suitable (N2)

  • < TS < 4: none suitable (N1)

In the second stage, the calculated weights will be used in TOPSIS and VIKOR to calculate suitability with an evaluation matrix. Once the evaluation matrix is prepared, it will be used by TOPSIS and VIKOR methods to calculate the suitability. Despite the methods use the same evaluation matrix, the suitability calculation methods which differ from each other.

For the purpose of constituting evaluating matrix, there are 600 test points specified in study area. The test points are specified by selecting the suitable areas for beekeeping considering water resources, agricultural sites, urban areas, industrial sites, and flora. Thus, the evaluation matrix includes 600 test points and related ranking values for 8 criteria. Each value for test points is assigned with extracting raster values to points. Test points distribution is given in Fig. 4.

Figure 4.
figure 4

Six hundred test points distribution for VIKOR and TOPSIS.

Full size image

The ranking values for each criterion are assigned between 1 and 9 considering beekeeping requirements (Table III).

Table III Evaluation matrix for TOPSIS and VIKOR
Full size table

After the evaluation matrix is prepared, TOPSIS technique is used to determine the most suitable test points. The TOPSIS technique is based on a model that the selected alternative should have the shortest distance from the positive ideal solution and the longest from the negative ideal solution. The respective distances to positive and negative ideal solutions are defined as a similarity index.

The ranking values are used to calculate R and V matrices via W criterion weights that calculated with AHP. The positive ideal solution A+ and the negative ideal solution A−, which are the maximum and minimum values of the V matrix, are calculated. Based on the A+ and A− values, distance to positive ideal solutions Di* and distance to negative ideal solution Di− values are calculated for each test point. Finally, relative closeness to ideal solution Ci*values are calculated (Table IV) to determine the beekeeping suitability ranking definition. The Ci* values are classified as follows;

  • 0.80 < Ci* < 1.00: highly suitable (S1)

  • 0.65 < Ci* < 0.80: moderately suitable (S2)

  • 0.50 < Ci* < 0.65: slightly suitable (S3)

  • 0.40 < Ci* < 0.50: moderately not suitable (N1)

  • 0.00 < Ci* < 0.40: none suitable (N2)

Table IV TOPSIS and VIKOR solution and rankings
Full size table

In addition to TOPSIS, VIKOR method calculation requires determining the best \( {f}_i^{\ast } \)and \( {f}_i^{-} \) values of all criteria functions. These values then used to determine the maximum group utility S and minimum individual regret to the opponent R values to be able to calculate the Q values. Determining the suitability of each test point depends to ranking alternatives sorting by the values S, R, and Q from the minimum value that represent the suitability order. The minimum Q value is assigned the highest suitable and the maximum Q value is none suitable. However, stability of the decision should be evaluated. Proposing Q value as a compromise solution, conditions 1 and 2 are evaluated. For acceptable advantage (C1), Q(a′′) − Q(a′) ≥ 1/(1 − 600) equation is not satisfied for this calculation. Thus, the condition 2 is satisfied for v = 0.5 value by concensus. According to the evaluation, the suitability is determined by ranking Q values and classified as follows;

  • 0.20 > Q > 0.00: highly suitable (S1)

  • 0.35 > Q > 0.20: moderately suitable (S2)

  • 0.50 > Q > 0.35: slightly suitable (S3)

  • 0.60 > Q > 0.50: moderately not suitable (N1)

  • 1.00 > Q > 0.60: none suitable (N2)

Suitability calculations for TOPSIS and VIKOR according to the Ci* and Q values are given in Table IV.

4 Results and discussion

The results indicate that 48% of the study area is assigned as suitable and that 52% of the study area is not suitable according to the AHP calculation. As can be seen in Table II, flora criterion have 44%, distance to waters 14.60%, and aspect have 10% weights in total weight ranking. It is possible to say that approximately 70% of suitability is defined by these classes. Because distance from settlements and distance from roads criteria does not have an effect on beekeeping suitability directly, these classes have 3% and 4% weights in total weight ranking. The suitability index maps are produced for each method respectively. The suitability index maps are generated by using TS, Ci* and Q values which are calculated with AHP, TOPSIS, and VIKOR. The suitability index maps are given in Fig. 5.

Figure 5.
figure 5

a AHP suitability. b TOPSIS suitability. c VIKOR suitability.

Full size image

The land suitability classes (in m2) are compared according to the 600 test points of the study area. The S1, S2, and S3 are assigned as suitable areas and N1, N2 as unsuitable. Considering the total suitability index values, the highest suitability rates are calculated by AHP and TOPSIS, respectively. The most suitable ranking (S1) is mostly assigned by AHP method and none suitable ranking (N2) by VIKOR. The comparisons are given in Fig. 6.

Figure 6.
figure 6

S1, S2, S3, N1, and N2 suitability (% of 600 test point) comparisons of AHP, TOPSIS, and VIKOR.

Full size image

The suitability ranking rates are compared according to the test points. Considering the S1, S2, and S3 are suitable rankings, 235 test points are assigned suitable by TOPSIS and 303 test points by VIKOR. For N1, N2 rankings, 365 test points by TOPSIS and 297 test points by VIKOR are assigned unsuitable (Table V).

Table V Beekeeping statistics for correlation analysis
Full size table

The effectiveness and reliability of the determined suitability can be verified in several ways such as considering existing beekeeping locations, evaluating with experts and testing suitable locations over the next year. The most reliable and rapid results can be obtained through a correlation analysis of the existing beekeeping statistics and determined suitability values. Thus, Turkish Statistical Institute 2015 beekeeping statistics are used to calculate correlation. According to the statistics, total honey production, total beehives and total beekeeper counts are available at district level and Bozkır, Hadim, Seydişehir, Beyşehir and Meram districts have the highest honey production rate which are also overlapped with the suitability maps. The statistic thematic maps are given in Fig. 7.

Figure 7.
figure 7

Thematic maps of Konya bee statistics.

Full size image

For the purpose of evaluating reliability and making a comparison, three different correlation analyses were generated. First correlation analysis was generated between beekeeping statistics (total honey production, total beekeepers and total hive) and determined suitability values for all districts. Suitability values of each district were calculated by zonal statistics of ArcGIS software. Only S1 and S2 suitability classes were included and % rate of the suitability were calculated. For instance, 60% of the Derbent district area (359 km2) were calculated as highly suitable for beekeeping.

According to the r values of correlation analysis, there is a good correlation with 0.70 r value between the total beekeeper count and VIKOR suitability values (AHP = 0.66 and TOPSIS = 0.67). The reason that total beekeeper count correlations are higher than others, the beekeepers must be registered to Republic of Turkey Ministry of Food, Agriculture and Livestock provincial directorates to be able to locate beehives. Thus, beekeeper counts represent the most real values. Total beehive correlations were calculated less than others, because beehive counts are not reflecting the real active beehive counts in production. It is impossible to know that how many beehives are being used actively in production by a beekeeper. The correlation graphics are given in Fig. 8.

Figure 8.
figure 8

Correlation graphics between MCDA methods and beekeeping statistics of Konya.

Full size image

Another correlation analysis was determined between MCDA methods by using 600 test point values. For this purpose, 600 test point AHP, TOPSIS, and VIKOR (TS, Ci, and Q) suitability values were assigned to each test point by extract values to points process and compared to each other. According to the correlation analysis of each method, there is a strong negative correlation between TOPSIS and VIKOR method. It is possible to say that all the calculations with AHP, TOPSIS, and VIKOR are consistent, and one of these methods can be used in beekeeping suitability. The correlation graphics are given in Fig. 9.

Figure 9.
figure 9

Correlation graphics of MCDA methods.

Full size image

As a different validation, for the purpose of determining the accuracy and reliability rate of these methods, existing beekeeper locations are retrieved from the Konya - Seydişehir, Beyşehir, Çumra, Hadim, and Taşkent Directorate of Provincial Food Agriculture and Livestock. Existing beekeeper locations are recorded between May and September 2016 with their attribute data such as beehive count, beekeeper name, address and honey type. One hundred seventeen existing beekeeper location coordinates are integrated to suitability maps to visualize the intersections in total. Distribution of the locations and intersection with the AHP suitability map are given in Fig. 10.

Figure 10.
figure 10

Existing beekeeper locations and AHP intersection.

Full size image

The values were retrieved from suitability map by using the intersection of locations and pixel values via ArcGIS software. Then, suitability values were converted to ranking values to determine the existing locations rankings for AHP, TOPSIS, and VIKOR suitability maps. Considering the S1, S2, and S3 suitability index values, AHP has 82%, VIKOR 88%, and TOPSIS 91% overlapping rates with the suitability maps. The rates are given in Fig. 11.

Figure 11.
figure 11

Comparison of the rankings with overlapping values.

Full size image

The intersection of existing beekeepers and ranking values indicated 2 and 11 beekeeper locations were intersected with non-suitable areas with VIKOR and TOPSIS method respectively. Best intersection for S1 class was determined by AHP method with 21 existing beekeeper locations.

5 Conclusions

Determining bee requirements and setting the optimum intervals to make decisions from complex alternatives are very difficult and inevitable process. This is the main reason of some limitations that encountered in this study. Increasing accuracy of this project can be succeed by involving additional criteria such as meteorological conditions, wind directions, flowering, foraging area, electromagnetic fields, and pesticide usage in agricultural lands. Due to the unavailability of temporal flora map and plant density information which are important for bees, there is a particular limitation when making decision accurately. Despite of some studies on biodiversity mapping projects started in Turkey, foraging and flowering maps are still unavailable because of unstable structure and difficulties to update and monitor annually.

For the purpose of eliminating limitations in the field of flowering season and foraging area, remote sensing methods could be used to detect the plant diversity and landscape. However, very short flowering season and small size of natural plants required high spatial resolution (even 5 cm) to be able to detect the flowers and plants. Moreover, remote sensing data must be retrieved temporarily within flowering season to detect all plants in the region. This can be possible in a small region due to the high cost of remote sensing data collection. Addition to this, meteorological and climatic conditions are varying from year to year and retrieved flowering season and plant maps could be changed next year. Thus, a rough classification of land use was used in this study.

Nevertheless, the results and validation of the suitability are quite satisfactory considering the 82% intersection rate of existing locations with suitability maps and correlation analysis with beekeeping statistics. The results also indicated that the weight calculation, interval settings of each criterion and ranking each interval according to the bee requirements are quite successful considering the intersection of existing beekeeper locations.

It is possible to say that this study can provide valuable and significant experience for beekeeping suitability projects when designing not only provincial but also national projects. Because Konya province contains all the characteristic features of Turkey in the field of topographic, climatic, meteorological and environmental features, this study can be accepted as a conceptual model of beekeeping suitability in Turkey. For the purpose of benefiting from this study as much as possible, sharing, visualizing and querying of suitability maps should be provided with interactive maps via Web Based Geographical Information Systems or GeoPortal systems. All the criteria, spatial data and suitability maps should be accessible to experts, beekeepers, researches, institutes and organizations to increase the interoperability and provide a platform for multi-disciplinary projects.