The main purpose of this study was modeling the suitability of recreational services in order to identify the high-priority areas for the development of recreational activities. To identify the spatial patterns of recreational services, the landscape attractiveness index, visit rate estimation, and access distribution were combined using a set of sub-indices. Moreover, the weight of each index was calculated using CRITIC weighting method and the importance of each index in the formation of recreational services spatial suitability was determined. The results showed that the presentation of highly recreational services is mainly done in the areas with forest ecosystems, protected areas, and aquatic ecosystems such as rivers and wetlands that have diverse landscapes. On the contrary, urban and human-built areas have the least potential for the provision of recreational services. Getis-Ord Gi* statistic was used to model the recreational suitability, identify recreational suitability clusters, and zoning. With its very low suitability, Zone 1 covers 19.9 percent of the area under study. This zone is mainly filled with agricultural lands (54.7%), grasslands (36.1%), and urban lands (6.8%). In contrast, with its very high recreational suitability, Zone 5 covers 16.5 percent of the area under study. The main uses of the lands in this zone are forest lands (55.5%) and aquatic layers including rivers and wetlands (17.2%).  Zone 3, which has an average suitability rate, covers more than half of the area under study. The approach presented in this study helps identify the development priorities in the light of the recreational potentials of the ecosystems existing in an area. Thus, it can be used to improve the spatial planning of recreation and tourism in order to preserve ecosystem services and sustainable use of them.