熊猫分布密度制图（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~