GDP區域分布圖的生成與對比（ArcPy實現）

一、背景

各地區經濟協調發展是保證國民經濟健康持續穩定增長的關鍵。GDP是反映各地區經濟發展狀況的重要指標。科學准確分析各地區GDP空間分布特征，對制定有效措施，指導經濟協調發展具有重要參考價值。

二、目的

ArcGis中提供了中提供了三種空間插值方法，每種插值方法在原理上和應用上都大不相同，在此通過具體實例練習如何利用內插方法和卸內插方法進行GDP空間分布特征的分析，以此來引導讀者對空間插值有一個更深刻的認識。

三、數據

某地區的統計GDP數據（GDP.shp）,數據范圍：4601萬~132630萬元。（\Chp8\Ex4\）

四、要求

*1) 經濟的發展具有一定的連帶效應和輻射作用。以該地區各區域年GDP數據為依據，采用IDW和Spline內插方法創建該地區GDP空間分異柵格圖。
*2）分析每種插值方法中主要參數的變化對內插結果的影響。
*3）分析兩種內插方法生成的GDP空間分布圖的差異性，簡單說明形成差異的主要原因。
*4）通過該練習，熟練掌握兩種插值方法的適用條件。

五、流程圖

• （1）反距離權重（IDW）內插方法流程：
  

• （2）樣條函數（Spline）內插方法流程：
  

• （3）兩種方法之間的對比，通過空間分析中【柵格計算器】命令來進行分析。選擇反距離權重(Power=2)，樣條函數（Regularized，Weight=0.01)。在【柵格計算器】中計算Abs (IDW-Spline)，比較二者產值變化空間分異特征。

六、模型構建器

七、ArcPy實現

注：柵格計算器工具專門用於 ArcGIS for Desktop 應用程序（僅作為 GP 工具對話框）或模型構建器。它不適用於腳本的編寫，而且也不能用於 ArcPy Spatial Analyst 模塊。所以給模型里的柵格計算器里的運算全部分成兩部，即相減並取絕對值。

# -*- coding: utf-8 -*-
# ---------------------------------------------------------------------------
# 8-4GDP區域分布圖的生成與對比.py
# Created on: 2021-10-10 20:21:39.00000
#   (generated by ArcGIS/ModelBuilder)
# Description: 
# ---------------------------------------------------------------------------

# Import arcpy module
import arcpy
from arcpy import Raster
from arcpy.sa import Abs
import os
import shutil

path = raw_input("請輸入數據所在的絕對路徑：").decode("utf-8")
paths = os.path.split(path)[0] + "\\result"
if not os.path.exists(paths):
    os.mkdir(paths)
else:
    shutil.rmtree(paths)
    os.mkdir(paths)

# Local variables:
GDP = path + "\\GDP.shp"
IDW2 = "IDW2"
IDW5 = "IDW5"
Spr0 = "Spr0"
Spr01 = "Spr01"
Spt0 = "Spt0"
Spt5 = "Spt5"


# Set Geoprocessing environments
print "Set Geoprocessing environments"
arcpy.env.workspace = paths
arcpy.env.scratchWorkspace = paths
arcpy.env.mask = path + "\\bound.shp"

# Process: 反距離權重法
print "Process: 反距離權重法"
arcpy.gp.Idw_sa(GDP, "GDP", IDW2, "500", "2", "VARIABLE 12", "")

# Process: 反距離權重法 (2)
print "Process: 反距離權重法 (2)"
arcpy.gp.Idw_sa(GDP, "GDP", IDW5, "500", "5", "VARIABLE 12", "")

# Process: Process: IDW2_IDW5相減並取絕對值
print "Process: IDW2_IDW5相減並取絕對值"
# arcpy.gp.RasterCalculator_sa("Abs(\"%IDW2%\" - \"%IDW5%\")", IDW2_IDW5)
Abs(Raster(IDW2) - Raster(IDW5)).save("IDW2_IDW5")

# Process: 樣條函數法
print "Process: 樣條函數法"
arcpy.gp.Spline_sa(GDP, "GDP", Spr0, "500", "REGULARIZED", "0", "12")

# Process: 樣條函數法 (2)
print "Process: 樣條函數法 (2)"
arcpy.gp.Spline_sa(GDP, "GDP", Spr01, "500", "REGULARIZED", "0.01", "12")

# Process: Process: Spr0_Spr01相減並取絕對值
print "Process: Spr0_Spr01相減並取絕對值"
# arcpy.gp.RasterCalculator_sa("Abs(\"%Spr0%\" - \"%Spr01%\")", Spr0_Spr01)
Abs(Raster(Spr0) - Raster(Spr01)).save("Spr0_Spr01")

# Process: 樣條函數法 (3)
print "Process: 樣條函數法 (3)"
arcpy.gp.Spline_sa(GDP, "GDP", Spt0, "500", "TENSION", "0", "12")

# Process: 樣條函數法 (4)
print "Process: 樣條函數法 (4)"
arcpy.gp.Spline_sa(GDP, "GDP", Spt5, "500", "TENSION", "5", "12")

# Process: Spt0_Spt5相減並取絕對值
print "Process: Spt0_Spt5相減並取絕對值"
# arcpy.gp.RasterCalculator_sa("Abs(\"%Spt0%\" - \"%Spt5%\")", Spt0_Spt5)
Abs(Raster(Spt0) - Raster(Spt5)).save("Spt0_Spt5")

# Process: Process: IDW2_Spr01相減並取絕對值
print "Process: IDW2_Spr01相減並取絕對值"
# arcpy.gp.RasterCalculator_sa("Abs(\"%IDW2%\" - \"%Spr01%\")", IDW2_Spr01)
Abs(Raster(IDW2) - Raster(Spr01)).save("IDW2_Spr01")

rasters = arcpy.ListRasters()
for raster in rasters:
    if u"_" not in raster:
        print u"正在刪除{}圖層".format(raster)
        arcpy.Delete_management(raster)
print "運行完畢~~~"

八、結果

實驗結束 byebye~