熊貓分佈密度製圖（ArcPy實現）

阿新 • • 發佈：2021-10-10

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

一、背景

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

二、目的

通過練習，熟悉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 fields.name.lower() in ["area", "power"]:
        arcpy.DeleteField_management(Xmpoint, fields.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 "執行完畢~~~"

注意看註釋提示！