在城市研究与空间规划中，POI（Point of Interest）数据是常见的基础资料。除了常规的核密度分析，我们往往还需要进一步量化不同区域的功能“丰富度”，以便更直观地比较哪些片区生活设施齐全，哪些区域功能单一。

本文介绍一种在 ArcGIS Pro Notebook 中一键运行的方法。通过一段脚本，可以自动完成以下步骤：

1. 构建分析单元

   • 支持生成规则网格（如 500m × 500m）

   • 或直接使用已有的多边形（如行政区、热点区）

2. 统计POI分类

   • 根据属性字段（如“中类”）统计每个分析单元内的设施分布

   • 可选过滤掉无意义类别（如“其他”）

3. 计算多样性指标

   • 种类数（Richness）：区域内设施类型数量

   • Shannon 指数：综合考虑类别数和均匀性

   • Simpson 指数：衡量类别占比是否过于集中

   • 均匀度（Evenness）：设施分布的均衡程度

4. 结果写回

   • 自动将结果字段写入分析单元要素类

   • 输出独立结果表，便于进一步统计或制图

运行完成后，用户即可在地图中直接以分级设色方式展示不同区域的设施丰富度，实现从数据到可视化的快速转换。

这种方法无需手动逐类统计，参数可灵活调整（网格大小、分析区域、过滤类别等），适合批量化、自动化的城市POI分析工作流。代码如下：

# -*- coding: utf-8 -*-
import arcpy, os, math

# ======= 路径（按需修改） =======
input_points = r"xx"   # 原始 POI 点（你的是 WGS84/4326）
grid_fc      = r"xxx"                         # 网格（现在是 Unknown）
out_points   = r"xxxx"          # 输出的新点

grid_id_field  = "grid_id"
metrics_fields = ["POI_TOT","RICHNESS","SHANNON","SIMPSON","EVENNESS"]

arcpy.env.overwriteOutput = True

def is_unknown(sr):
    try:
        name = (sr.name or "").strip().lower()
        wkid = sr.factoryCode
        return (not sr) or (name in ("", "unknown", "未知")) or (wkid in (None, 0, -1))
    except:
        return True

def guess_and_define_projection(fc):
    """当坐标系 Unknown 时，根据范围粗略推断 3857/4326 并 Define Projection"""
    sr = arcpy.Describe(fc).spatialReference
    if not is_unknown(sr):
        print(f"✔ 网格已有坐标系：{sr.name} (WKID={sr.factoryCode})，无需Define。")
        return sr

    ext = arcpy.Describe(fc).extent
    xs = [ext.XMin, ext.XMax]
    ys = [ext.YMin, ext.YMax]
    width  = abs(ext.XMax - ext.XMin)
    height = abs(ext.YMax - ext.YMin)

    # 简单推断：经纬度范围通常 |X|<=180, |Y|<=90；Web Mercator 在中国 |X|大约 1e6~2e7
    # 另外，鱼网是按“米”生成的：width/height 通常在几千~几万以上
    looks_like_projected = (max(map(abs, xs)) > 1000 or max(map(abs, ys)) > 1000) and (max(width, height) > 1000)

    if looks_like_projected:
        # 绝大多数情况下，我们用 3857（你的网格是由脚本用“米”建的）
        target = arcpy.SpatialReference(3857)  # WGS 1984 Web Mercator Auxiliary Sphere
        arcpy.management.DefineProjection(fc, target)
        print(f"⚙ 已为网格 Define Projection 为：{target.name} (WKID={target.factoryCode})")
        return target
    else:
        target = arcpy.SpatialReference(4326)  # WGS 1984 Geographic
        arcpy.management.DefineProjection(fc, target)
        print(f"⚙ 已为网格 Define Projection 为：{target.name} (WKID={target.factoryCode})")
        return target

def get_count(ds):
    try:
        return int(arcpy.management.GetCount(ds)[0])
    except:
        return -1

# 1) 修复网格坐标系（Define Projection）
grid_sr = guess_and_define_projection(grid_fc)

# 2) 若点与网格坐标系不同，则把“点”投影到“网格”坐标系
pts_sr = arcpy.Describe(input_points).spatialReference
scratch = arcpy.env.scratchGDB or arcpy.env.workspace
if not scratch:
    scratch = arcpy.CreateFileGDB_management(os.path.dirname(grid_fc), "tmp_scratch.gdb").getOutput(0)

pts_proj = os.path.join(scratch, "input_pts_proj_to_gridSR")
if arcpy.Exists(pts_proj):
    arcpy.management.Delete(pts_proj)

if (pts_sr.factoryCode or -1) != (grid_sr.factoryCode or -2) and (not is_unknown(pts_sr)) and (not is_unknown(grid_sr)):
    arcpy.management.Project(input_points, pts_proj, grid_sr)
    pts_for_join = pts_proj
    print(f"✔ 已将点投影到：{grid_sr.name} (WKID={grid_sr.factoryCode})")
else:
    pts_for_join = input_points
    print("ℹ 点与网格坐标系一致，或有一侧未知（但我们刚刚已为网格定义），直接相交。")

# 3) 空间连接：仅把 grid_id 拼到点
pts_with_gid = os.path.join(scratch, "poi_with_gridid_tmp")
if arcpy.Exists(pts_with_gid):
    arcpy.management.Delete(pts_with_gid)

# 只从网格带 grid_id
fm = arcpy.FieldMappings()
fm.addTable(pts_for_join)

fms_grid = arcpy.FieldMappings()
fms_grid.addTable(grid_fc)
idx = fms_grid.findFieldMapIndex(grid_id_field)
if idx == -1:
    raise RuntimeError(f"在网格中未找到字段 {grid_id_field}，请确认网格带有唯一ID。")
fm.addFieldMap(fms_grid.getFieldMap(idx))

arcpy.analysis.SpatialJoin(
    target_features=pts_for_join,
    join_features=grid_fc,
    out_feature_class=pts_with_gid,
    join_operation="JOIN_ONE_TO_ONE",
    join_type="KEEP_ALL",
    field_mapping=fm,
    match_option="INTERSECT"
)

cnt_all = get_count(pts_with_gid)
cnt_gid = get_count(arcpy.management.MakeFeatureLayer(pts_with_gid, "lyr_gid", f"{grid_id_field} IS NOT NULL"))
print(f"步骤A：共 {cnt_all} 个点，其中 {cnt_gid} 个点获得了 {grid_id_field}")

# 4) 按 grid_id 把多样性指标 JoinField 到点
if arcpy.Exists(out_points):
    arcpy.management.Delete(out_points)
arcpy.management.CopyFeatures(pts_with_gid, out_points)

# 确认网格上确有指标字段
grid_fields = [f.name for f in arcpy.ListFields(grid_fc)]
missing = [f for f in metrics_fields if f not in grid_fields]
if missing:
    raise RuntimeError("网格缺少以下指标字段（请先运行多样性计算脚本写入）： " + ", ".join(missing))

arcpy.management.JoinField(
    in_data=out_points, in_field=grid_id_field,
    join_table=grid_fc, join_field=grid_id_field,
    fields=metrics_fields
)

cnt_nz = get_count(arcpy.management.MakeFeatureLayer(out_points, "lyr_nz", "POI_TOT IS NOT NULL"))
print(f"步骤B：输出 {out_points}，其中 POI_TOT 非空点数 = {cnt_nz}/{cnt_all}")
print("✅ 完成。打开属性表检查 grid_id 与 POI_TOT / RICHNESS / SHANNON / SIMPSON / EVENNESS。")

• 想用你已有的热点区/行政区作单元：把 use_existing_polygons=True，并把 existing_polygon_fc 指向你的多边形数据；脚本会自动投影到米制坐标后进行统计。

• 想改变网格尺度：改 cell_size_m（如 300、1000 等）。

• 想保留“其他”类别：把 exclude_categories = []。

• 生成结果后，在 ArcGIS Pro 里对 SHANNON 或 EVENNESS 做分级设色即可直观看到“设施丰富性”空间格局。