熊貓分布密度制圖（ArcPy實現）

一、背景

大熊貓是我國國家級珍惜保護動物，熊貓的生存必須滿足一定槽域（獨占的獵食與活動范圍）條件。因此，科學准確的分析熊貓的分布情況，對合理制定保護措施和評價保護成效具有重要意義。

二、目的

通過練習，熟悉ArcGIS密度制圖函數的原理及差異性，掌握如何根據事跡采樣數據特點，結合ArcGIS提供的密度制圖功能和其它空間分析，制作符合要求的密度圖。

三、要求

1）熊貓活動具有一定的槽域范圍，一個槽域范圍只有一個或一對熊貓，假設熊貓槽域半徑為5000m，理論上最大槽域面積為3.1450005000

2）雖然一個采樣點代表一個熊貓，但由於熊貓的生存具有確定槽域特征，不同的采樣點具有不同的空間控制面積。假定熊貓活動范圍分布滿足以采樣點為中心的泰森多邊形，如何將這一信息加入密度分布圖是本練習的重點。

3）在野外實采的熊貓活動足跡數據的基礎上，以每個熊貓槽域范圍為權重，運用ArcGIS 中的區域分配功能和密度制圖功能制作該地區熊貓分布密度圖。

四、數據

野外實采的熊貓活動足跡數據，一個足跡代表一個熊貓曾在此處活動過，相同足跡只記載一次。（\Chp8\Ex3）

五、計算原理

首先利用柵格數據空間分析模塊提供的區域分配功能提取熊貓的槽域范圍，然后用理論最大域面積（假定是半徑為5km的圓，面積為3.141592755 km²）除以所提取的熊貓實際槽域面積，作為采樣點的加權值（記為Power字段），生成熊貓分布密度圖。

解題思路：
1）歐式分配，求出每個熊貓的實際控制范圍

2）添加字段area，計算每個熊貓的實際控制面積：柵格數柵格邊長柵格邊長

3）添加power字段，計算權重：理論上最大槽域面積/柵格數柵格邊長柵格邊長

4）利用power字段進行核密度分析

操作流程圖：

六、模型構建器

七、ArcPy實現

# -*- coding: utf-8 -*-
# ---------------------------------------------------------------------------
# 8-3熊貓密度圖.py
# Created on: 2021-10-10 15:53:34.00000
#   (generated by ArcGIS/ModelBuilder)
# Description: 
# ---------------------------------------------------------------------------

# Import arcpy module
import arcpy
import os


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

# Local variables:
Xmpoint = path + "\\Xmpoint.shp"  # 數據
Output_distance_grid = "Distance"  # 輸出距離柵格,格網基本名稱的長度不能超過了 13
Output_direction_grid = "Direction"  # 輸出方向柵格,格網基本名稱的長度不能超過了 13
Output_inverse_grid = "Inverse"  # 輸出反向柵格,格網基本名稱的長度不能超過了 13
EucAllo_shp = "EucAllo_shp"  # 歐式分配輸出要素
KernelD_shp = "KernelD_shp"  # 核密度分析輸出要素
constant = "10000000"
Distribution_map = "密度分布圖"  # 格網基本名稱的長度不能超過了 13

# Set Geoprocessing environments
print "Set Geoprocessing environments"
arcpy.env.workspace = paths  # 工作空間
arcpy.env.scratchWorkspace = paths  # 臨時工作空間
arcpy.env.extent = "15153765.390904 -140950.922884 15227181.528600 -92036.748014"  # 處理范圍,左下右上

# Process: 檢查屬性表字段是否含有"area", "power"，有則刪除
fields = arcpy.ListFields(Xmpoint)
for field in fields:
    if field.name.lower() in ["area", "power"]:
        arcpy.DeleteField_management(Xmpoint, field.name)

# Process: 歐氏分配
print "Process: 歐氏分配"
arcpy.gp.EucAllocation_sa(Xmpoint, EucAllo_shp, "5000", "", "500", "ID", Output_distance_grid, Output_direction_grid, "PLANAR", "", Output_inverse_grid)

# Process: 添加字段
print "Process: 添加字段"
arcpy.AddField_management(EucAllo_shp, "Area", "LONG", "", "", "", "", "NULLABLE", "NON_REQUIRED", "")

# Process: 計算字段
print "Process: 計算字段"
arcpy.CalculateField_management(EucAllo_shp, "Area", "[COUNT] * 500 *500", "VB", "")

# Process: 連接字段
print "Process: 連接字段"
arcpy.JoinField_management(Xmpoint, "ID", EucAllo_shp, "Value", "Area")

# Process: 添加字段 (2)
print "Process: 添加字段 (2)"
arcpy.AddField_management(Xmpoint, "power", "LONG", "", "", "", "", "NULLABLE", "NON_REQUIRED", "")

# Process: 計算字段 (2)
print "Process: 計算字段 (2)"
arcpy.CalculateField_management(Xmpoint, "power", "3.1415927*5000*5000 / [Area]", "VB", "")

# Process: 核密度分析
print "Process: 核密度分析"
arcpy.gp.KernelDensity_sa(Xmpoint, "power", KernelD_shp, "500", "", "SQUARE_MAP_UNITS", "DENSITIES", "PLANAR")

# Process: 乘
print "Process: 乘"
arcpy.gp.Times_sa(KernelD_shp, constant, Distribution_map)

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

注意看注釋提示！

八、結果

實驗結束 byebye~