学校选址（ArcPy实现）

一、背景

合理的学校空间位置布局，有利于学生的上课与生活。学校的选址问题需要考虑地理位置、学生娱乐场所配套、与现有学校的距离间隔等因素，从总体上把握这些因素能够确定出适宜性比较好的学校选址区。

二、目的

通过练习，熟悉ArcGIS栅格数据距离制图、成本距离加权、数据重分类、多层面合并等空间分析功能﹔熟练掌握利用ArcGIS空间分析功能，分析和结果类似学校选址的实际应用问题。

三、数据

(1) landuse (土地利用图);
(2) dem(地面高程图);
(3) rec _ sites(娱乐场所分布图);
(4) school(现有学校分布图)。
（\Chp8 \ Ex1 \）

四、要求

(1）新学校选址需注意如下几点:①新学校应位于地势较平坦处;②新学校的建立应结合现有土地利用类型综合考虑，选择成本不高的区域;③新学校应该与现有娱乐设施相配套，学校距离这些设施越近越好﹔④新学校应避开现有学校，合理分布。
(2）各数据层权重比为:距离娱乐设施占0.5，距离学校占0.25，土地利用类型和地势位置因素各占0.125。
(3）实现过程运用ArcGIS的扩展模块( extension）中的空间分析（ spatial ana-lyst）部分功能，具体包括:坡度计算、直线距离制图功能、重分类及栅格计算器等功能完成。
(4)给出适合新建学校的适宜地区图，并作简要分析。

五、实现流程图

ArcGIS中实现学校选址分析，首先利用现有学校数据集、现有娱乐场所数据集和高程数据派生出坡度数据，以及到现有学校、娱乐场所距离数据集。然后重分类数据集到相同的等级范围，再按照上述数据集在学校选址中的影响率赋权重值，最后合并这些数据即可创建显示新学校适宜位置分布的地图。学校选址的逻辑过程主要包括四个部分(如下图所示，结果中深色部分为学校候选区)。
(1）数据准备，确定需要哪些数据作为输人，包括高程数据（dem)、土地利用数据(landuse)、现有学校数据（school)娱乐场所数据（rec_ sites) ;
(2）派生数据集，从现存数据派生出能提供学校选址的原始成本数据，包括坡度数据、到现有学校距离数据集和到娱乐场所数据集;
(3）重分类各种数据集，消除各成本数据集的量纲影响，使各成本数据具有大致相同的可比分类体系。各成本数据均按等间距分类原则分为1~10级，级数越高适宜性越好﹔
(4）给各数据集赋权重。必要的话在适宜性模型中影响比较大的数据集赋比较高权重，然后合并各数据集合以寻找适宜位置。

六、模型构建器

七、ArcPy实现

# -*- coding: utf-8 -*-
# ---------------------------------------------------------------------------
# 8-1 学校选址.py
# Created on: 2021-10-12 13:19:17.00000
#   (generated by ArcGIS/ModelBuilder)
# Description: 
# ---------------------------------------------------------------------------

# Import arcpy module
import arcpy
from arcpy import Raster
import os
import shutil
import time

print time.asctime()
path = raw_input("请输入数据所在文件夹的绝对路径(如果是数据库，则带上数据库名)：").decode("utf-8")
if "db" in path:
    paths = os.path.split(path)[0] + '\\result'
else:
    paths = path + "\\result"
# 开始计时
time_start = time.time()
if not os.path.exists(paths):
    os.mkdir(paths)
else:
    shutil.rmtree(paths)
    os.mkdir(paths)

# Local variables:
rec_sites = path + "\\rec_sites"
dem = path + "\\dem"
landuse = path + "\\landuse"
school = path + "\\school"
# 温馨建议：输出文件的名称不能超过13个字符
Direct_grid1 = "Direct_grid1"
Reverse_grid1 = "Reverse_grid1"

Direct_grid2 = "Direct_grid2"
Reverse_grid2 = "Reverse_grid2"
Dist_recsites = "Dist_recsites"
Re_recsite = "Re_recsite"
Dist_schools = "Dist_schools"
Re_schools = "Re_schools"

Slope_dem = "Slope_dem"
Re_slope = "Re_slope"

land_use = "land_use"
Re_landuse = "Re_landuse"
site_result = "学校选址结果"

# Set Geoprocessing environments
print "Set Geoprocessing environments"
arcpy.env.scratchWorkspace = paths
arcpy.env.workspace = paths
arcpy.env.extent = landuse
arcpy.env.cellSize = dem
arcpy.env.mask = landuse

# Process: 要素转栅格
print "Process: 要素转栅格"
arcpy.FeatureToRaster_conversion(landuse, "Land_use", land_use, "5")

# Process: 重分类
print "Process: 重分类"
arcpy.gp.Reclassify_sa(land_use, "Land_use", "'city center' 3;transnational 7;public 10;vegetable 4;agriculture 5;forest 2;'barren land' 1", Re_landuse, "DATA")

# Process: 坡度
print "Process: 坡度"
arcpy.gp.Slope_sa(dem, Slope_dem, "DEGREE", "1", "PLANAR", "METER")

# Process: 重分类 (2)
print "Process: 重分类 (2)"
arcpy.gp.Reclassify_sa(Slope_dem, "Value", "0.108518 6.291646 10;6.291646 10.751648 9;10.751648 15.131917 8;15.131917 19.819941 7;19.819941 24.931349 6;24.931349 30.529348 5;30.529348 36.494091 4;36.494091 42.837433 3;42.837433 50.637821 2;50.637821 71.992943 1", Re_slope, "DATA")

# Process: 欧氏距离
print "Process: 欧氏距离"
arcpy.gp.EucDistance_sa(rec_sites, Dist_recsites, "", dem, Direct_grid1, "PLANAR", "", Reverse_grid1)

# Process: 欧氏距离 (2)
print "Process: 欧氏距离 (2)"
arcpy.gp.EucDistance_sa(school, Dist_schools, "", dem, Direct_grid2, "PLANAR", "", Reverse_grid2)

# Process: 重分类 (3)
print "Process: 重分类 (3)"
arcpy.gp.Reclassify_sa(Dist_recsites, "Value", "0 61.846584 10;61.846584 111.018013 9;111.018013 158.113876 8;158.113876 206.155289 7;206.155289 256.320099 6;256.320099 309.273346 5;309.273346 366.230804 4;366.230804 427.463440 3;427.463440 492.544403 2;492.544403 614.837402 1", Re_recsite, "DATA")

# Process: 重分类 (4)
print "Process: 重分类 (4)"
arcpy.gp.Reclassify_sa(Dist_schools, "Value", "0 93.941475 10;93.941475 155.724121 9;155.724121 212.896683 8;212.896683 269.258240 7;269.258240 326.955658 6;326.955658 388.490662 5;388.490662 454.009918 4;454.009918 521.176575 3;521.176575 598.539917 2;598.539917 738.782776 1", Re_schools, "DATA")

# Process: 栅格计算器
# arcpy.gp.RasterCalculator_sa("\"%Re_recsite%\" * 0.5 + \"%Re_schools%\" * 0.25 + \"%Re_slope%\" * 0.125 + \"%Re_landuse%\" * 0.125", site_result)

# Process: 乘运算、加运算
print "Process: 乘运算、加运算"
(Raster('Re_recsite') * 0.5 + Raster('Re_schools') * 0.25 + Raster('Re_slope') * 0.125 + Raster('Re_landuse') * 0.125).save(site_result)

save = [u"学校选址结果"]
rasters = arcpy.ListRasters()
for raster in rasters:
    if raster.lower() not in save:
        print u"正在删除{}图层".format(raster)
        arcpy.Delete_management(raster)
# 结束计时
time_end = time.time()
# 计算所用时间
time_all = time_end - time_start
print time.asctime()
print "执行完毕！>>><<< 共耗时{:.0f}分{:.2f}秒".format(time_all // 60, time_all % 60)