學校選址（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)

八、結果

實驗結束 byebye~