章节概述

本章围绕空间点数据，学习用 sf 整理位置与属性，并结合制图、聚合与空间统计识别分布模式。

曼哈顿 Airbnb：整理房源的每夜报价，学习价格分布、空间聚合，以及 Moran’s I 与 LISA 的全局和局部关联分析。
上海微博与公共交通：用核密度估计（KDE）提取高密度 AOI，比较 AOI 内外记录到最近交通设施的距离。

章节练习引入 DBSCAN，观察 eps 与 minPts 如何改变聚类；小组实践侧重参数敏感性与结果解释。

本章学习内容

学习导航

本章学习内容

欢迎开始 第 5 章 的学习。请点击 下方导航卡 进入相应小节：

阅读提示：通过返回主页或主页按钮，可回到本章学习导航。

1.0 空间点数据

1.0 空间点数据概论

基本概念

第四章介绍了 sf 矢量表：每行是一条要素记录，geometry 列保存几何，其余列保存属性。本章关注其中的点记录。点的位置由坐标和坐标参考系（CRS）共同确定。

坐标表示位置，其他字段说明对象或事件。 同一坐标可能对应几家店，也可能对应同一地点在不同时刻的多次签到。坐标因此不能充当记录的唯一编号；一条记录还可能包含地点名称、事件类型、时间和数据来源。

读点图前，先辨认观察单位

地点：一家店、一座车站或一个设施入口。
事件：某地某时的一次投诉、事故或签到。
移动位置：一条轨迹中的某个采样点，或一次行程的起点、终点。

路灯可以直接记录为点；居民小区或城市在小比例尺地图上也常用代表点表示，此时点只概括位置，不能表示对象的完整范围。点位可能来自设备定位，也可能来自地址地理编码（geocoding），即将地址转换为坐标。坐标精度和覆盖范围取决于数据的采集方式。

从真实记录读点

下面选取纽约中城六座 Citi Bike 站点，将同一批记录放在表和地图中。字母 A—F 仅用于两边对照；每行对应一座站点。capacity 是站点的车位容量，不代表当时可借出的自行车数量。

六座 Citi Bike 站点的属性与几何
记录	站点名称	车位	`geometry`（便于阅读的表示）
A	W 35 St & 8 Ave	45	`POINT (-73.992805 40.752762)`
B	E 40 St & 5 Ave	93	`POINT (-73.981632 40.752062)`
C	W 42 St & 8 Ave	41	`POINT (-73.990985 40.757570)`
D	W 44 St & 5 Ave	68	`POINT (-73.980144 40.755003)`
E	W 49 St & 8 Ave	67	`POINT (-73.987882 40.762272)`
F	W 51 St & 6 Ave	81	`POINT (-73.980420 40.760660)`

纽约中城的六座 Citi Bike 站点（A—F 与左表对应）

geometry 是 sfc_POINT 几何列，表中的 POINT (x y) 是便于阅读的文本表示。这里 x = 经度、y = 纬度，CRS 为 EPSG:4326。点的几何维度是 0，没有长度或面积；定位仍需两个坐标值。站点资料来自 Citi Bike 官方站点信息（2026 年 10 月 14 日快照），背景道路来自 OpenStreetMap。

点分布与统计单元

原始点记录可以逐条绘制。若把同一批点汇总到行政区或网格，统计结果还会受到单元划分方式的影响。可变面元问题（Modifiable Areal Unit Problem, MAUP）包括两种常见情形：尺度效应指单元变大或变小时，局部差异可能被平滑或放大；分区效应指单元数量与大小相近，但边界位置改变后，各单元的计数、均值或相关性仍可能改变。

同一批模拟点：面元大小与边界改变格内计数

布鲁克林一块 $6\times6$ km 的研究窗口内有 380 座 Citi Bike 站点。下图把同一批点汇总到不同网格，统一以每平方千米的站点数表示。

同一批 380 座站点：500 m 与 1.5 km 网格下的站点密度

尺度改变时，500 m 网格有 10 个零站点单元；合并为 1.5 km 网格后，16 个单元都包含站点，局部空白随之消失。再保持每格面积为 $3\text{ km}^2$，分别采用 $2\times1.5$ km 与 $1.5\times2$ km 的矩形网格，比较边界划法：

同一批 380 座站点：单元面积相同、网格方向不同时的站点密度

两组图使用同一批站点与同一研究范围，只改变统计单元。为演示分区的影响，暂以超过 12 站/km²标记高密格：横向划分有 3 格达到阈值，纵向划分仅有 1 格。12 站/km² 是本图的示例阈值，不能直接推广为站点密度标准。站点数据来自 Citi Bike 官方站点信息（2026 年 10 月 14 日快照）。

选举中的分区效应（模拟）

九个相邻小格各投出 10 张有效选票，A、B 代表两方。图中票数均按 A∶B 书写；细线区分小格，粗线围出选区。每区由三个小格组成，合计后票数较多的一方获得该区一席。下面只改变分区方向，选票总数始终为 A 48、B 42。

同一批模拟选票：按行或按列划区，区内合计票数决定席位

按行划区：三行各为 A 16∶B 14，A 获得 3 席。

按列划区：前两列各为 A 14∶B 16，第三列为 A 20∶B 10；A 获得 1 席，B 获得 2 席。

左图最下行虽有两个小格由 B 领先，合计后仍由 A 胜出。小格底色不能直接当作选区胜负。同一批选票汇总到不同选区，席位分配便可能改变；解读按区汇总的结果时，需要说明分区方案。

展开核对选区票数

以左图最下行为例：A 获得 $4+4+8=16$ 票，B 获得 $6+6+2=14$ 票，因此这一选区由 A 获得一席。换成右图第一列，A 为 $5+5+4=14$ 票，B 为 $5+5+6=16$ 票，因此这一选区由 B 获得一席。

1.1 常见的空间点记录

餐厅是地点，在餐厅发生的一次签到是事件。二者即使坐标相同，记录单位也不同。以下按记录对象与形成方式介绍四类常见点数据。

1. 地点与设施点（POI）

兴趣点（Point of Interest, POI）通常表示一家商店、一座车站、一个公园入口等可定位的地点。名称、类别、营业信息或服务评分可以作为属性。若对象本身有面积，地图上的点通常是入口或代表位置，并不等于它的完整范围。

POI 可用于观察设施分布、服务覆盖和不同功能的邻近关系。平台收录规则与更新频率会影响点数。Inside Airbnb 的房源点位经过位置偏移，不能当作精确门牌位置。

下面四幅图截取伦敦、纽约、巴黎和东京各一块约 $1.6\times1.6$ km 的街区，按相同规则绘出 OpenStreetMap 中标注的咖啡馆、图书馆与博物馆。

伦敦 · 特拉法加广场周边

纽约 · 曼哈顿中城

巴黎 · 卢浮宫周边

东京 · 上野周边

地图数据：© OpenStreetMap contributors，提取时间为 2026 年 10 月。

2. 行为与事件点

事件记录通常说明何时、何地发生了什么。一次纽约市 311 非紧急市政服务请求、交通事故或社交媒体签到都可以作为点记录；时间信息有时精确到分钟，有时只有日期。位置也可能由地址或路口转换而来，因此不一定是事件发生的精确坐标。

同一地点可以反复产生记录。按时段绘图可观察记录数量与分布的变化；平台用户构成、上报意愿和数据缺漏会影响观察到的空间分布。

事件点的位置与发生时刻（模拟示意）

图中的街区和事件均为模拟资料。八个时间层覆盖一天，相邻时间层间隔三小时；高度表示发生时刻，不表示事件强度。

3. 轨迹采样点与起讫点

轨迹由同一移动对象按时间排序的位置采样点组成。两个采样点之间的实际路径通常没有被直接记录；采样间隔、定位误差和信号缺失都会影响路线判断。起讫记录（OD）则成对保存一趟行程的起点（Origin）与终点（Destination），可能另有出发时间、到达时间或交通方式。

OD 可以由完整轨迹简化，也可以来自租借、进出站等系统记录。地图上连接起终点的直线只表示两地之间存在行程，不代表实际经过的道路。

两条模拟轨迹的采样点与同一起讫关系

4. 从其他空间对象得到的代表点

有些点由行政区、建筑轮廓等面对象计算而来，用于摆放标签、表示区域位置或连接区域间的 OD 关系。选择代表点时，应说明它概括的是形状、人口还是某个出入口。

表示方式	常见用途	解释时注意
几何重心 `st_centroid()`	概括形状的几何位置	凹形或带孔的面，其重心可能落在面外。
面内点 `st_point_on_surface()`	在面上放置符号或标签	点位于面上，但未必接近居民或活动集中的位置。
人口加权质心	表示人口分布的重心	依赖人口及其空间分布数据，也可能落在区域边界外。

延伸阅读 · 这些点数据从哪里来？

地点与事件：OpenStreetMap 的公开要素可通过 osmdata 与 Overpass 查询；NYC 311 数据页提供市政服务请求记录。
行程：Citi Bike 公开逐次骑行的起终点站记录；NYC TLC 的现行出租车公开数据主要使用出租车分区标识起终地，不能当作连续 GPS 轨迹。
平台数据：社交媒体签到由用户发布，商业 POI 由平台收录；数据开放程度和空间覆盖因平台而异。

本章的两个案例采用不同的记录单位：专题研究 1 先按房源 ID 汇总跨期报价，再比较六边形区域；专题研究 2 以社交媒体发帖记录估计活动分布。判断记录是否重复，需要结合 ID、时间和研究对象，不能仅凭坐标相同就删除。

2.0 点数据实战 1

2.0 专题研究 1：曼哈顿 Airbnb 每夜报价

专题研究 1：从房源点到空间关联

分析导航

曼哈顿 Airbnb 空间探索流程

本例从曼哈顿主岛的三期公开记录出发，检查有效每夜报价，按房源 ID 求均价并筛选均价的 P5–P95，再汇总为六边形中位报价。随后定义邻居与权重，用 Moran’s I 检验整体关联，用 LISA 定位局部高低报价组合。

研究课题

2.1 研究问题与 ESDA 框架

Airbnb 的每夜报价在城市内部如何分布，相邻地区的报价是否相似？本例以曼哈顿主岛为研究区，将整套住宅、私人房间等挂牌房源汇总到六边形，描述区域报价差异，再检验整体与局部空间关联。房型和入住条件与报价差异的关系留待后续分组研究。

数据来自 Inside Airbnb 的三个公开发布批次。price 记录每夜报价，单位为美元／夜。研究关注这些批次中可观测房源的挂牌报价及其空间差异。

本例的研究问题

在 2026 年 3、6、9 月三个发布批次中，曼哈顿主岛的每夜报价如何分布？汇总到六边形后，中位报价是否呈现相邻关联？显著的局部高低值组合出现在哪里？

ESDA 如何组织这一研究？

探索性空间数据分析（exploratory spatial data analysis, ESDA）将地图、分布图和空间统计结合起来，用于描述空间差异、识别集聚与离群，并形成后续研究假设。Anselin、Syabri 与 Kho（2006）对这类方法及其联系作了系统介绍。本例按下列框架展开：

分析环节	要回答的问题	本例的实现
研究对象与数据检查	一条记录代表什么，价格字段是否有效？	主岛边界、房源 ID、价格缺失与分布检查。
可视化与空间聚合	差异在哪里，采用什么统计单元？	按 ID 求均价后筛选 P5–P95；归入 0.10 km² 六边形，至少 5 个不同 ID 才计算格内中位数。
空间权重	哪些格算作“周围”？	有效六边形的 Queen 邻接与行标准化权重 $W$。
全局关联	相近报价是否倾向于相邻？	全局 Moran’s I、散点图与随机置换检验。
局部关联	哪些格与周围相似，哪些形成反差？	局部条件置换检验；以原始 p 值为主，另比较 BH、BY 等调整。
反思与再分析	哪些判断依赖样本或尺度选择？	检查价格变量、筛选范围、网格尺度、ID 门槛与空间权重的影响。

点图首先呈现房源的位置与报价，Moran 分析进一步检验六边形中位报价的空间关联。房源数量与区域中位报价是两个变量；房源密集的地方未必报价较高。

分析单位

原始表的一行是一条房源—批次记录；汇总后的一个点是一个房源 ID，即使几个 ID 共用公开坐标也分别保留；空间统计的一条观测是一个六边形。后续每一步均标明所用的统计单位。

先认识三个方法概念

概念	含义	在本例中如何使用
空间权重 $W$	明确单元间的邻接关系及权重。	将“周围”落实为可计算的邻居集合。
全局 Moran’s I	概括整个研究区中属性值与其空间滞后的关联。	判断中位报价相近的六边形是否倾向于相邻。
LISA 与局部 Moran	局部空间关联指标（local indicators of spatial association, LISA）是一类指标；本例采用其中的局部 Moran。	区分显著的同向集聚与空间离群格。

方法延展 · 空间权重、Moran 散点图与 LISA

1. 把“邻近”写成空间权重。 Rook 邻接要求两个面共享边，Queen 邻接允许共享边或顶点；距离阈值与最近邻规则则按位置间的距离确定邻居。下图用方格说明两种接壤规则。

方格上的 Queen 与 Rook 邻接；浅蓝色表示中心格的邻居。

规则六边形密铺中，Queen 与 Rook 通常得到相同邻接；本例的有效格也已核对为相同。距离阈值或最近邻规则可提供不同的“周围”定义。完整说明见 GeoDa 接壤权重教程。

本例采用行标准化：若第 $i$ 格有 $k_i$ 个邻居，则各邻居的权重为 $1/k_i$，其他格的权重为零。令 $z_i$ 为第 $i$ 格标准化后的中位报价，空间滞后（spatial lag）为

\[ (Wz)_i=\sum_j w_{ij}z_j. \]

因此，$Wz$ 可理解为邻居标准化报价的加权平均。本例先检查 365 个达到 ID 门槛的完整六边形，去掉 3 个没有 Queen 邻居的格，再对其余 362 格重建邻接与行标准化权重。地图的海岸遮罩只用于显示。

2. 全局 Moran’s I 与散点图。 一般形式为

\[ I=\frac{n}{S_0}\frac{\sum_i\sum_j w_{ij}z_i z_j}{\sum_i z_i^2}, \qquad S_0=\sum_i\sum_j w_{ij}. \]

其中 $n$ 是统计计算使用的有效样本量，$S_0$ 是全部权重之和。正向关联表示相似数值倾向于相邻，负向关联表示相邻单元更常呈现反差。随机重排属性值时，期望为 $-1/(n-1)$，在样本量较大时接近零；统计量的可达范围取决于 $W$。

本例在剔除孤立格后，重新计算 362 格报价的均值与标准差。所有权重行和均为 1，因而 $n=S_0=362$；以 $z$ 为横轴、$Wz$ 为纵轴，通过原点拟合的斜率等于本设定下的 Moran’s I。全局置换检验固定格网位置与 $W$，随机重排报价，比较观测关联是否超出这种随机排列下常见的程度。参见 Moran（1950）的经典文献与 GeoDa 全局空间自相关教程。

3. 从四象限走向局部检验。 Anselin（1995）提出了 LISA 框架。下表中的高低均相对于本次有效格的整体均值；纵轴对应周围格的加权平均。

组合	该格	周围格	含义
High–High（HH）	高	高	相对高值的同向集聚。
Low–Low（LL）	低	低	相对低值的同向集聚。
High–Low（HL）	高	低	高值格与周围形成反差。
Low–High（LH）	低	高	低值格与周围形成反差。

象限先描述方向，局部检验再判断是否达到显著性门槛。本例采用局部条件置换的双侧模拟 $p$ 值，主图显示原始 $p\leq0.05$ 的局部类别；另以相同的原始 $p$ 值比较 BH、BY、Holm 与 Bonferroni 调整，检查显著范围对方法的敏感性。方法说明见 GeoDa 局部空间自相关教程；R 的具体实现见 spdep 的局部置换文档。

延伸阅读 · Airbnb、租差与城市住房研究

空间探索可以衔接更具体的城市研究问题：高报价集聚区与旅游设施是否重合？房型构成是否解释部分地区差异？平台住宿供给与长期租赁、社区变化有怎样的关系？

纽约与租差。 Wachsmuth 与 Weisler（2018）的 Airbnb and the rent gap: Gentrification through the sharing economy 将纽约 Airbnb 的地理差异与租差、绅士化联系起来。可阅读期刊原文入口或作者公开稿，了解如何把空间格局转化为城市研究问题。租差指现有土地使用获得的地租与潜在更高地租之间的差额。
住房价格与租金。 Barron、Kung 与 Proserpio（2021）的 The Effect of Home-Sharing on House Prices and Rents: Evidence from Airbnb 研究平台住宿与住房市场的关系。论文入口可用于比较：检验价格关联，与识别 Airbnb 对租金或住房价格的影响，需要哪些不同资料和研究设计？

本例完成的是报价格局的 ESDA。上述机制可作为后续研究方向，需要补充住房、房型、设施及时间变化等资料再检验。

研究区域与数据

2.2 研究区与三期数据

纽约市有五个行政区。本例仅取曼哈顿行政区中面积最大的连续多边形，即曼哈顿主岛；罗斯福岛等离岛不在分析范围内。定位附图显示主岛在纽约市的位置。

曼哈顿主岛研究区及其在纽约市的位置；右侧定位图标出五个行政区。

认识研究地点

纽约市（New York City, NYC）由 Manhattan、Brooklyn、Queens、The Bronx 和 Staten Island 五个行政区（boroughs）组成，位置见上图右侧。曼哈顿主岛与曼哈顿行政区的范围不同：行政区还包括 Marble Hill 和若干小岛。五区关系可查阅 NYC311 官方说明，行政区范围见纽约市行政法典。

主岛位于 Hudson River 与 East River 之间，岛体狭长；中央公园、街网和沿岸水域形成清楚的地理参照。研究限定在主岛，既便于观察岛内报价差异，也便于明确分析边界。

城市延展 · 从纽约影像到曼哈顿街景

先从国际空间站视角辨认曼哈顿及周边水域，再进入时代广场与中央公园，观察商业街道、公共绿地和高层建筑的空间差异。以下四张真实照片摄于 2025—2026 年，分别提供城市区域与街道尺度的背景。

纽约日景｜2026-02-13。国际空间站视角。NASA/Chris Williams，iss074e0315625。

2025年国际空间站拍摄的纽约与新泽西夜景，曼哈顿灯光带、中央公园暗区及河流形成对照。

纽约夜景｜2025-09-11。城市灯光与水域。NASA，iss073e0658285。

时代广场｜2025-08-31。娱乐商业区的屏幕与人流。Nielsoncaetanosalmeron，CC BY 4.0。

中央公园与中城｜2025-05-19。从大都会艺术博物馆望出。Boaventuravinicius，CC BY 4.0。

原图与说明：NASA 日景 · NASA 夜景 · 时代广场 · 中央公园。NASA 影像按其教育与信息用途指引使用；后两张来自 Wikimedia Commons，按 CC BY 4.0署名。网页仅调整图片尺寸和编码压缩，保留原始画面比例。

Airbnb 与 Inside Airbnb 提供什么资料？

Airbnb 是连接房东与旅客的住宿平台，房源包含整套房源、私人房间等类型。Inside Airbnb 是独立整理平台公开信息的项目，提供房源、日历与评论等数据下载。一次发布是当时可见信息的快照：同一房源可以在不同批次出现，房源上下架及信息更新也会改变样本。

本例使用房源位置与公开报价。Inside Airbnb 的数据说明指出，公开坐标经过匿名化，与实际地址可能相差约 0–150 m。公开坐标不等于房源身份：不同 ID 即使坐标相同，仍作为不同房源保留；同一 ID 跨批次的报价才汇总。

Inside Airbnb 在 2026 年 1—9 月公开的纽约市详细房源文件中，可用于本例的是 3 月 16 日、6 月 14 日和 9 月 14 日三个发布批次。原始文件见 Inside Airbnb 数据目录，再发布数据时须核对其数据使用政策。

字段与记录	本例的读取方式	用途
`id` 与发布批次	同一批次中每个房源 ID 至多一条记录；同一房源可跨批次出现。	检查重复记录及批次构成。
`longitude`、`latitude`	房源的公开经纬度。	限定主岛、确定每个 ID 的代表位置及空间聚合。
`price`	本次公开报价折算的每夜金额。	价格检查、筛选与汇总。
`price_quote_checkin_date`、`price_quote_checkout_date`	入住日与退房日；两者相减得到报价夜数。	辅助理解每次报价的入住日期和时长。
`minimum_nights`、`room_type`	最低入住限制及房型。	理解报价条件与开展分组分析。

发布日期、入住日期与每夜报价

发布批次记录数据的发布时间；报价所对应的入住、退房日期由 price_quote_checkin_date 与 price_quote_checkout_date 给出，两者之差为报价夜数。实际入住日期可能跨到下一年。

minimum_nights 是房源的最低入住限制；本次报价夜数应由报价所列的入住、退房日期计算。price 已是每夜报价，无需再除以入住夜数。

各次报价的入住日期、时长、房型和折扣可能不同，可用这些字段进一步分组比较。三期汇总概括的是三个快照中的挂牌报价，不能解释为连续九个月的成交记录。

下列代码优先读取本机缓存；缓存不存在时从官方地址下载。地图边界优先读取本例保存的曼哈顿主岛边界，缺失时从纽约市边界服务获取。

先检查价格字段。 下表将三期原表与曼哈顿主岛记录分别按价格状态计数。“价格缺失”包括空白值；“解析失败”指非空价格无法转换为数值。

三期 Airbnb 价格字段状态
统计范围	记录数	价格缺失	解析失败	负价	零价	正价
纽约市三期原表	96,156	33,800	0	0	0	62,356
曼哈顿主岛	43,589	15,853	0	0	0	27,736

这三期数据中没有负价或零价。主岛内 43,589 条记录中，15,853 条价格缺失，其余 27,736 条为有效正价记录。下方地图与下一节的原价、对数分布图均使用这 27,736 条记录。

三张地图分别展示各发布批次的有效正价记录，共用地图范围、Fisher 七级断点与色盘。分类在合并三期报价的 log10 尺度上计算，图例以美元／夜表示。

正文的分析模块重建样本、网格、权重与统计量，图件读取同一研究流程导出的浅深色图片。可下载专题研究 1 分析脚本，在第5章目录下按模块运行；当前图件与输入核对说明记录图源及数据指纹。研究区定位的制图实现见定位脚本。运行时需保留本章的数据和缓存目录。

2026 年三期曼哈顿主岛的 27,736 条有效正价记录；三幅地图共用 Fisher 七级断点。

查看并复制数据读取、研究区与价格检查 R 代码

模块 1：依赖声明

# 模块 1：声明本专题依赖；缺少的包只需在 Console 安装一次。
# install.packages(c("dplyr", "readr", "sf", "spdep", "ggplot2", "classInt"))
library(dplyr)
library(readr)
library(sf)
library(spdep)
library(ggplot2)
library(classInt)

模块 2：三期文件读取与字段解析

# 模块 2：按发布批次读取原始表，解析每夜报价与报价日期。
# 缓存不存在时下载；此处不按报价夜数筛选。
data_dir <- "figures/manhattan_2026_multisnapshot_prototype/data"
release_dates <- c("2026-03-16", "2026-06-14", "2026-09-14")
source_base <- "https://data.insideairbnb.com/united-states/ny/new-york-city"

read_release <- function(release_date) {
  file_path <- file.path(data_dir, paste0("nyc_", release_date, ".csv.gz"))
  if (!file.exists(file_path)) {
    dir.create(data_dir, recursive = TRUE, showWarnings = FALSE)
    source_url <- paste0(source_base, "/", release_date, "/data/listings.csv.gz")
    download.file(source_url, file_path, mode = "wb")
  }
  quote_fields <- c("id", "longitude", "latitude", "price", "room_type",
                    "minimum_nights", "price_quote_checkin_date",
                    "price_quote_checkout_date")
  read_csv(file_path,
           col_select = all_of(quote_fields),
           col_types = cols(.default = col_character(),
                            longitude = col_double(), latitude = col_double()),
           show_col_types = FALSE, progress = FALSE) %>%
    mutate(snapshot = release_date,
           price_usd = parse_number(price,
                                    locale = locale(grouping_mark = ",")),
           quote_checkin = as.Date(price_quote_checkin_date),
           quote_checkout = as.Date(price_quote_checkout_date),
           quote_nights = as.integer(quote_checkout - quote_checkin))
}

listing_rows <- bind_rows(lapply(release_dates, read_release))

模块 3：研究区边界

# 模块 3：读取主岛边界；缺失时从官方曼哈顿多边形取最大陆地主体。
boundary_path <- file.path(data_dir, "manhattan_main_island_original_boundary.gpkg")
if (file.exists(boundary_path)) {
  main_island <- st_read(boundary_path, quiet = TRUE)
} else {
  boundary_url <- paste0(
    "https://services5.arcgis.com/GfwWNkhOj9bNBqoJ/arcgis/rest/services/",
    "NYC_Borough_Boundary/FeatureServer/0/query?where=1=1&outFields=*&outSR=4326&f=pgeojson"
  )
  main_island <- st_read(boundary_url, quiet = TRUE) %>%
    filter(BoroName == "Manhattan") %>%
    st_make_valid() %>%
    st_cast("POLYGON") %>%
    mutate(area_ft2 = as.numeric(st_area(st_transform(geometry, 2263)))) %>%
    slice_max(area_ft2, n = 1, with_ties = FALSE) %>%
    select(BoroName)
  st_write(main_island, boundary_path, quiet = TRUE)
}
main_island <- st_transform(main_island, 4326)

模块 4：坐标与主岛范围

# 模块 4：先检查坐标，再确定主岛内记录；保留报价缺失以便统计。
# 先检查原表，再用有效坐标确定主岛内的记录；此时尚未筛掉价格缺失值。
valid_location <- is.finite(listing_rows$longitude) &
  is.finite(listing_rows$latitude)
located_rows <- listing_rows[valid_location, ]
on_main_island <- lengths(st_intersects(
  st_as_sf(located_rows, coords = c("longitude", "latitude"), crs = 4326),
  main_island
)) > 0L
main_island_rows <- located_rows[on_main_island, ]

模块 5：报价状态检查

# 模块 5：分别报告原表与主岛记录的价格缺失、解析和正负状态。
count_price_status <- function(rows, scope) {
  price_missing <- is.na(rows$price) | trimws(rows$price) == ""
  price_valid <- is.finite(rows$price_usd)
  tibble::tibble(
    scope = scope,
    records = nrow(rows),
    missing = sum(price_missing),
    parse_failed = sum(!price_missing & !price_valid),
    negative = sum(price_valid & rows$price_usd < 0),
    zero = sum(price_valid & rows$price_usd == 0),
    positive = sum(price_valid & rows$price_usd > 0)
  )
}
price_eda <- bind_rows(
  count_price_status(listing_rows, "纽约市三期原表"),
  count_price_status(main_island_rows, "曼哈顿主岛")
)
price_eda %>%
  mutate(across(-scope, ~ format(.x, big.mark = ",", trim = TRUE))) %>%
  knitr::kable(
    format = "html",
    caption = "三期 Airbnb 价格字段状态",
    col.names = c("统计范围", "记录数", "价格缺失", "解析失败",
                  "负价", "零价", "正价"),
    align = c("l", rep("r", 6))
  )

模块 6：有效正价

# 模块 6：保留可用且大于零的每夜报价，供分布检查与 ID 汇总。
positive_rows <- main_island_rows %>%
  filter(is.finite(price_usd), price_usd > 0)

制图之后，按什么顺序读图？

读图先确认统计对象，再识别模式，最后结合检验形成判断。后续各图均可沿用下列顺序。

顺序	重点查看	在本例中的检查
先认对象	研究范围、记录单位与时间。	主岛边界；房源—批次记录、房源 ID 或六边形；三个发布批次。
再读编码	色阶、单位、参照值与未纳入区域。	美元／夜；颜色对应的价格区间；少于 5 个房源 ID 的格另列。
观察模式	高低值在哪里，是否相邻，哪里有反差。	沿岛体比较南北及东西差异，并关注成片区域与孤立点。
检验判断	视觉模式是否得到统计支持。	先用 Moran’s I 概括整体关联，再用局部检验定位显著组合。

本图的基本解读。 三期点图均可见报价的空间差异：北段较多较低报价点，中南段混合分布着更多较高报价点。点位较密的地方存在遮挡，后续用六边形聚合后再比较区域模式。三个发布批次中的房源构成与实际入住日期不同，跨期价格变化还需在相同报价条件下另行比较。

数据预处理

2.3 价格分布、房源均价与筛选

预处理按以下顺序执行。 前一节完成范围与价格检查，本节接着形成每个 ID 的代表报价，再做分位筛选；后续才进入网格和空间检验。

顺序与对应位置	处理内容	输出及统计单位
1 · 模块 2–4	读取三期快照，保留有限经纬度并限定曼哈顿主岛。	43,589 条房源—批次记录。
2 · 模块 5–6	检查缺失与解析状态，只保留有限且大于零的每夜报价。	27,736 条有效正价记录；据此绘制原价与对数分布。
3 · 模块 7	检查 ID 非空、同一 ID 在同一批次中唯一；按 ID 求已观测报价的算术均值，位置取最新有效正价记录。	11,853 个房源 ID，各有一个均价和代表位置。
4 · 模块 8–9	在上述 ID 均价上计算 P5、P95，按包含两端的区间筛选并核对数量。	10,669 个 ID，均价约为 63.96–1,000 美元／夜。
5 · 2.4 节	将每个 ID 唯一归入完整六边形；每格至少 5 个不同 ID 才计算格内中位数。	653 格中有 365 个达标格。
6 · 2.5–2.6 节	检查邻居，排除 3 个孤立格，在其余格上重建权重并重新标准化，再运行 Moran 与 LISA。	两项检验共用 362 个参检格及同一套权重。

缺失检查针对分析所需的坐标、报价和 ID；其他字段缺失不构成整行删除的理由。同一 ID 跨批次出现是预期的重复观测，由均价汇总；同一 ID、同一批次若重复，代码会报错，须回查来源。不同 ID 共用坐标时仍分别保留。

先按 ID 求均价，再筛 P5–P95。 这一顺序使每个房源只向分位数计算贡献一个值。若先在逐条报价中截尾，部分房源的观测会先被剔除，随后得到的均价和入选 ID 都可能改变。

先看分布。 下图比较主岛内 27,736 条有效正价记录的原价与 log10 分布。原价图保留金额差异；对数图压缩高价长尾，便于观察主体分布，相同横向距离表示相同的价格倍数变化。log10 用于分布展示，均价与 P5–P95 均在原始美元／夜尺度上计算。

两幅直方图的柱高均表示相应区间内的记录数，核密度曲线（kernel density estimate, KDE）概括平滑后的分布，底部观测短线（rug）标示每条有效正价记录的位置。KDE 分别在原价与 log10 尺度上估计，再乘以样本量与各自组距，使曲线与频数轴对应。两图的组距位于不同尺度，应比较分布形状，而不直接比较对应柱高或曲线高度。

两图的中位数参考线都表示 237.75 美元／夜，在对数图中放置于 log10(237.75) 的位置。水平网格线便于读取频数；这些参照线与 KDE、rug 均不参与价格筛选。

27,736 条有效正价记录的原价与对数分布、KDE 及 rug；两图中位数参考线均为 237.75 美元／夜，此时尚未按 ID 汇总或筛选 P5–P95。

按房源 ID 求均价。 对 27,736 条有效正价记录按 id 分组，计算各房源已观测报价的算术平均，得到 11,853 个房源均价。每个 ID 无论出现一、二还是三期，都使用同一个 mean() 规则；未观测批次不补零，不同 ID 即使坐标相同也分别保留。均价可平滑部分已观测报价波动，仍会受到入住条件与季节差异的影响。

在房源均价层筛选 P5–P95。 对 11,853 个均价计算第 5 与第 95 百分位，得到 63.96–1,000 美元／夜的保留区间。分位点采用 R 默认的 type = 7，筛选直接使用 quantile() 返回的完整精度值，并包含上下界。均价低于 P5 和高于 P95 的房源各有 592 个，最终保留 10,669 个房源 ID。每个 ID 只贡献一个均价，不会因出现批次较多而在分位数计算中重复计权。

处理阶段	统计单位	保留数量
主岛内原始记录	房源—批次记录	43,589
有效正价记录	房源—批次记录	27,736
按 ID 求均价	不同房源 ID	11,853
均价 P5–P95 筛选后	不同房源 ID	10,669

11,853 个房源均价上的 P5–P95：保留 10,669 个房源 ID；阈值按原始美元均价计算。

地图展示 P5–P95 保留的 10,669 个房源 ID。每个点代表一个 ID 的已观测均价，位置取该 ID 最新一条有效正价记录的公开坐标。坐标经过匿名化，不能解释为精确地址；不同 ID 共用坐标时，点位可能重叠。

右侧直方图展示全部 11,853 个 ID 均价。P5–P95 内的部分沿用地图的 Fisher 七级颜色；低于 P5 和高于 P95 的两端各有 592 个 ID，以灰色表示。Fisher 断点仍由筛选后的均价在 log10 尺度上计算，图例以美元／夜表示。

直方图的 KDE 与 rug 均使用全部 11,853 个均价。KDE 在 log10 尺度上估计，并按样本量与组距缩放到频数轴；rug 标出每个 ID 均价的位置。这里的 KDE 描述一维价格分布，其高度不表示原美元尺度上的概率密度；专题研究 2 将学习位置空间中的核密度。

地图展示 10,669 个保留 ID；直方图展示全部 11,853 个 ID 均价，彩色部分与地图共用 Fisher 七级色阶，灰色两端各为 592 个剔除 ID。

查看并复制同 ID 均价与 P5–P95 筛选 R 代码

模块 7：ID 检查与代表报价

# 模块 7：按房源 ID 汇总全部有效正价，计算原始美元每夜报价的均值。
# 单期、两期和三期均使用同一规则；未出现的批次不补零。
stopifnot(!anyNA(positive_rows$id),
          !any(positive_rows$id == ""),
          !anyDuplicated(positive_rows[c("id", "snapshot")]))
listing_mean_all <- positive_rows %>%
  arrange(id, snapshot) %>%
  group_by(id) %>%
  summarise(longitude = dplyr::last(longitude),
            latitude = dplyr::last(latitude),
            listing_mean_usd = mean(price_usd),
            n_releases = n_distinct(snapshot),
            n_records = n(), .groups = "drop")

模块 8：均价分位筛选

# 模块 8：每个 ID 贡献一个均价；在均价层计算 P5–P95，包含上下界。
# 分位点使用完整精度；显示时的四舍五入不参与筛选。
price_bounds <- quantile(listing_mean_all$listing_mean_usd,
                         probs = c(0.05, 0.95), type = 7)
price_lower <- unname(price_bounds[1])
price_upper <- unname(price_bounds[2])
listing_clean <- listing_mean_all %>%
  filter(listing_mean_usd >= price_lower,
         listing_mean_usd <= price_upper)
stopifnot(nrow(positive_rows) == 27736L,
          nrow(listing_mean_all) == 11853L,
          nrow(listing_clean) == 10669L,
          !anyDuplicated(listing_clean$id),
          all(listing_mean_all$n_records == listing_mean_all$n_releases))

模块 9：阶段数量核对

# 模块 9：按处理阶段报告数量，并明确记录与房源 ID 的单位变化。
data.frame(
  stage = c("主岛内原始记录", "有效正价记录", "按 ID 求均价",
            "均价 P5–P95 筛选后"),
  unit = c("房源—批次记录", "房源—批次记录", "不同房源 ID", "不同房源 ID"),
  count = c(nrow(main_island_rows), nrow(positive_rows),
            nrow(listing_mean_all), nrow(listing_clean))
) %>% knitr::kable(col.names = c("处理阶段", "统计单位", "保留数量"),
                   format.args = list(big.mark = ","))

读图顺序：先看分布，再看筛选后的空间模式。

辨认横轴与样本。 第一组图比较 27,736 条有效正价记录的原价与对数分布；均价筛选图以 11,853 个房源 ID 为单位。对数轴中相同的横向距离表示相同的价格倍数变化。
按图辨认标记。 初步原价／对数图用 KDE 显示平滑形状、用 rug 标示记录位置，并标出同一个美元中位数。均价筛选图以柱体与 P5–P95 参照展示保留范围；最终地图旁的 KDE 与 rug 基于全部 11,853 个 ID 均价。
核对筛选范围，再回到地图。 地图展示保留的 10,669 个 ID。直方图中保留区间的彩色部分与地图共用 Fisher 断点；灰色两端对应被剔除的 ID，各有 592 个。

本图的发现。 27,736 条有效正价记录呈右偏长尾，中位数为 237.75 美元／夜，少量高价记录将右尾拉长。按 ID 求均价并筛选 P5–P95 后，样本包含 10,669 个房源，其已观测均价处于 63.96–1,000 美元／夜。后续将结合六边形聚合，继续比较报价的空间差异。

初步分布图的一条观测是房源—批次记录，均价筛选图与最终点图的一条观测是一个房源 ID。同 ID 跨已观测批次取均值会平滑该房源的部分报价变化；它不等于消除了季节变化，也不代表连续九个月的成交均价。

网格化与聚合

2.4 六边形网格与中位报价

逐点设色存在点位叠压，难以直接比较区域的典型报价。本节将 2.3 节保留的 10,669 个房源 ID 归入固定面积网格：每个 ID 携带一个已观测均价和最新有效正价记录的公开坐标，同坐标的不同 ID 分别参与汇总。

本例采用面积 $0.10\text{ km}^2$、边长约 196 m 的完整正六边形，并要求每格至少包含 5 个不同 ID，再计算格内房源均价的中位数。面积、网格起点与 ID 门槛都是需要说明的分析设定，可通过改变设定比较结果。

距离与面积计算使用 EPSG:2263，其单位为美国测量英尺，按 $1\text{ m}=3937/1200\text{ ft}$ 换算。st_make_grid() 的六边形 cellsize 是对边距离，本例约为 339.8 m。统计使用与主岛相交的完整格；地图才将格网裁切到主岛边界，并在 EPSG:4326 下显示。

与主岛相交的 653 格中，365 格至少包含 5 个不同 ID；其余 288 格不计算中位报价，不以 0 补值。n_coordinates 记录不同公开坐标数，帮助检查空间覆盖，但不作为纳入门槛。listing_members 保留每个入选 ID 的报价与格编号，可用于后续样本支持诊断。

下图左侧显示这 365 格的中位报价，右侧直方图、KDE 与 rug 展示同一批格的数值分布。Fisher 七级仅用于设色，后续 Moran 与 LISA 使用未分级、未取对数的美元中位报价。

0.10 km² 完整六边形的中位报价：365 格至少包含5个不同房源ID；288格不足门槛，价格不补零。地图与直方图共用Fisher七级色阶。

查看并复制六边形网格分析 R 代码

library(sf)
library(dplyr)
library(spdep)

# 1. 网格构建：使用上一节的入选房源与原始主岛边界。
stopifnot(nrow(listing_clean) == 10669L,
          !anyDuplicated(listing_clean$id),
          all(is.finite(listing_clean$listing_mean_usd)),
          all(listing_clean$listing_mean_usd > 0))
main_island_2263 <- st_transform(main_island, 2263)
target_area_m2 <- 100000
feet_per_metre <- 3937 / 1200
cellsize_ft <- sqrt(2 * target_area_m2 / sqrt(3)) * feet_per_metre
hex_side_m <- sqrt(2 * target_area_m2 / (3 * sqrt(3)))

# 沿用主岛原始外包框的默认起点，先编号，再保留相交的完整格。
grid_origin_ft <- unname(st_bbox(main_island_2263)[c("xmin", "ymin")])
stopifnot(max(abs(grid_origin_ft -
  c(978886.7706350009, 194226.33300118599))) < 1e-6)
hex_grid <- st_make_grid(main_island_2263,
                         cellsize = cellsize_ft,
                         offset = grid_origin_ft,
                         square = FALSE) %>%
  st_sf() %>%
  mutate(hex_id = row_number()) %>%
  st_filter(main_island_2263, .predicate = st_intersects)
stopifnot(nrow(hex_grid) == 653L, all(st_is_valid(hex_grid)))

current_esda_dir <- paste0(
  "figures/manhattan_2026_multisnapshot_prototype/output/",
  "current_esda_20261026"
)

# 2. 空间归属：每个 ID 只归入一个格，不按坐标删除房源。
listing_points_2263 <- st_as_sf(listing_clean,
  coords = c("longitude", "latitude"), crs = 4326,
  remove = FALSE) %>% st_transform(2263)
hex_hits <- st_intersects(listing_points_2263, hex_grid)
within_hits <- st_within(listing_points_2263, hex_grid)

# 当前数据中，每个点恰好落入一个完整格；出现零归属或多归属时先检查。
stopifnot(all(lengths(hex_hits) == 1L),
          all(lengths(within_hits) == 1L),
          identical(vapply(hex_hits, `[[`, integer(1), 1),
                    vapply(within_hits, `[[`, integer(1), 1)))
listing_members <- listing_clean %>%
  mutate(hex_id = hex_grid$hex_id[
    vapply(within_hits, `[[`, integer(1), 1)])
stopifnot(nrow(listing_members) == 10669L,
          !anyDuplicated(listing_members$id),
          !anyNA(listing_members$hex_id))

# 3. 先数不同 ID，再对达到门槛的格计算中位数。
hex_support <- listing_members %>%
  group_by(hex_id) %>%
  summarise(n_ids = n_distinct(id),
            n_coordinates = n_distinct(longitude, latitude),
            .groups = "drop")
hex_medians <- listing_members %>%
  semi_join(filter(hex_support, n_ids >= 5L), by = "hex_id") %>%
  group_by(hex_id) %>%
  summarise(median_usd = median(listing_mean_usd), .groups = "drop")
hex_values <- left_join(hex_support, hex_medians, by = "hex_id")
hex_all <- hex_grid %>%
  left_join(hex_values, by = "hex_id") %>%
  mutate(n_ids = coalesce(n_ids, 0L),
         n_coordinates = coalesce(n_coordinates, 0L),
         eligible = n_ids >= 5L)
hex_eligible <- hex_all %>% filter(eligible) %>% arrange(hex_id)
stopifnot(nrow(hex_all) == 653L, nrow(hex_eligible) == 365L,
          sum(hex_all$n_ids) == 10669L,
          all(is.na(hex_all$median_usd[!hex_all$eligible])),
          all(is.finite(hex_eligible$median_usd)))

制图核心。 使用 365 格的连续中位报价计算 Fisher 七级，裁切只用于显示。下列代码省略底图和比例尺，完整浅深色样式见本节图源。

# 4. Fisher七级与显示裁切；将figure_mode改为"dark"可查看深色配色。
figure_mode <- "light"
price_colours <- list(
  light = c("#471063", "#453882", "#375A8C", "#2A788E",
            "#1F958B", "#2CB17E", "#62CB5F"),
  dark = c("#AB9AF5", "#83A9F5", "#63C5E8", "#60DBCB",
           "#A7E38F", "#F4D773", "#FFAB70"))
fisher_breaks <- classInt::classIntervals(hex_eligible$median_usd,
  n = 7, style = "fisher", largeN = Inf)$brks
price_labels <- paste0(round(head(fisher_breaks, -1)), "–",
                       round(tail(fisher_breaks, -1)))
hex_plot <- hex_all %>%
  mutate(price_class = cut(median_usd, breaks = fisher_breaks,
    include.lowest = TRUE, right = TRUE, labels = price_labels)) %>%
  st_intersection(st_union(main_island_2263)) %>% st_transform(4326)
ggplot2::ggplot(hex_plot) +
  ggplot2::geom_sf(ggplot2::aes(fill = price_class),
    colour = NA) +
  ggplot2::scale_fill_manual(values = price_colours[[figure_mode]],
    na.value = if (figure_mode == "light") "#D9E1E8" else "#354A60",
    drop = FALSE, name = "美元／夜") +
  ggplot2::coord_sf(datum = NA) + ggplot2::theme_void() +
  ggplot2::theme(
    plot.background = ggplot2::element_rect(colour = NA,
      fill = if (figure_mode == "light") "#F8FBFE" else "#101D2D"),
    text = ggplot2::element_text(
      colour = if (figure_mode == "light") "#1D3448" else "#F0F6FB"))

读图顺序：样本支持 → 代表报价 → 相邻差异。

先看支持范围。 有报价色阶的格包含至少 5 个不同 ID；灰色格不足门槛，不能读作零报价或低报价。
再比较报价。 颜色与直方图表示格内各 ID 均价的中位数。北段较多格处于低报价等级，中南段有较多高报价组合。
观察相邻格。 相近色阶的连续分布和局部反差提供空间线索，是否超出随机排列下的常见程度还需检验。

网格将房源层面的均价转换成区域层面的代表值。每个 ID 贡献一次；同坐标的不同 ID 仍分别参与计数和中位数计算。Fisher 颜色仅用于读图，空间检验直接使用连续报价。

全局空间自相关

2.5 全局 Moran’s I

全局 Moran’s I 检查相似报价是否倾向于相邻。本例先对 365 个达标完整格构建 Queen 邻接，其中 3 格没有邻居。将这 3 格排除后，在其余 362 格上重建权重，并用这批格的均值与标准差重新计算标准化报价。365 格仍是达到房源数量门槛的上游样本，362 格才是空间检验样本。

## Warning in poly2nb(hex_eligible, queen = TRUE, row.names = as.character(hex_eligible$hex_id)): some observations have no neighbours;
## if this seems unexpected, try increasing the snap argument.

## Warning in poly2nb(hex_eligible, queen = TRUE, row.names = as.character(hex_eligible$hex_id)): neighbour object has 5 sub-graphs;
## if this sub-graph count seems unexpected, try increasing the snap argument.

## Warning in poly2nb(analysis_hex, queen = TRUE, row.names = as.character(analysis_hex$hex_id)): neighbour object has 2 sub-graphs;
## if this sub-graph count seems unexpected, try increasing the snap argument.

Queen 把共边或共顶点的格定义为邻居；行标准化权重 W 的各行和为 1，对角线为 0。当前参检图有 2 个连通部分，没有孤立格。构建邻接时关于多个子图的提示反映这一拓扑结构，可结合邻居列表核对。

$I=0.621$，单侧置换 $p=0.0001$。 在当前报价口径、0.10 km² 网格和 Queen 行标准化权重下，相近中位报价的格倾向于相邻，表现为全局正空间自相关。$I$ 衡量空间关联，其数值不表示已解释的价格差异比例。

同一批362个有Queen邻居的完整格：横轴为标准化中位报价z，纵轴为邻居的行加权均值Wz；金色过原点线斜率为0.620857。颜色仅表示四象限组合，不表示局部显著性。

查看并复制全局空间自相关 R 代码

# 4. 邻居与孤立格：先识别无邻居格，再重建参检格的W。
neighbors_365 <- poly2nb(hex_eligible, queen = TRUE,
                        row.names = as.character(hex_eligible$hex_id))
queen_isolates <- card(neighbors_365) == 0L
isolated_hex_ids <- hex_eligible$hex_id[queen_isolates]
stopifnot(sum(queen_isolates) == 3L,
          identical(as.integer(isolated_hex_ids), c(1042L, 1389L, 1596L)))
analysis_hex <- hex_eligible[!queen_isolates, ]
neighbors <- poly2nb(analysis_hex, queen = TRUE,
                    row.names = as.character(analysis_hex$hex_id))
w <- nb2listw(neighbors, style = "W", zero.policy = FALSE)
stopifnot(nrow(analysis_hex) == 362L, all(card(neighbors) > 0L),
          n.comp.nb(neighbors)$nc == 2L,
          identical(as.character(analysis_hex$hex_id), attr(w, "region.id")))

# 5. 标准化与全局检验：使用连续的原美元中位报价。
hex_price <- analysis_hex$median_usd
z <- as.numeric(scale(hex_price))
wz <- lag.listw(w, z, zero.policy = FALSE)
analysis_n <- length(z)
weight_sum <- Szero(w)
global_i <- analysis_n / weight_sum * sum(z * wz) / sum(z^2)

# 行标准化且没有孤立格时n=S0；过原点拟合线斜率等于Moran I。
origin_fit <- lm(wz ~ 0 + z)
stopifnot(abs(weight_sum - analysis_n) < 1e-12,
          abs(unname(coef(origin_fit)[1]) - global_i) < 1e-12,
          abs(global_i - 0.62085739252026317) < 1e-10)

RNGkind("Mersenne-Twister", "Inversion", "Rejection")
set.seed(20261024)
global_permutation <- moran.mc(hex_price, w, nsim = 9999,
  alternative = "greater", zero.policy = FALSE)
stopifnot(abs(unname(global_permutation$statistic) - global_i) < 1e-12,
          abs(global_permutation$p.value - 0.0001) < 1e-12)

制图核心。 横轴为同一批362格的标准化报价，纵轴为邻域加权均值；四象限仅表示正负组合。过原点线使用本节计算出的全局斜率。

# 6. 沿用z与wz；纵轴不再单独标准化。
figure_mode <- "light"
quad_colours <- list(
  light = c(HH = "#B83F4D", LH = "#087D69", LL = "#2867A4", HL = "#7B4AA7"),
  dark = c(HH = "#FF9A94", LH = "#65D7BF", LL = "#8ABCF6", HL = "#C4A4F4"))
scatter_values <- data.frame(hex_id = analysis_hex$hex_id, z = z, wz = wz) %>%
  mutate(quadrant = case_when(z >= 0 & wz >= 0 ~ "HH",
    z < 0 & wz >= 0 ~ "LH", z < 0 & wz < 0 ~ "LL", TRUE ~ "HL"))
ggplot2::ggplot(scatter_values, ggplot2::aes(z, wz, colour = quadrant)) +
  ggplot2::geom_hline(yintercept = 0, colour = "grey60", linetype = 2) +
  ggplot2::geom_vline(xintercept = 0, colour = "grey60", linetype = 2) +
  ggplot2::geom_point(size = 2, alpha = 0.75) +
  ggplot2::geom_abline(intercept = 0, slope = global_i,
    colour = if (figure_mode == "light") "#956207" else "#F6D380") +
  ggplot2::scale_colour_manual(values = quad_colours[[figure_mode]]) +
  ggplot2::labs(x = "标准化中位报价 z", y = "邻域加权均值 Wz",
    colour = "象限组合") + ggplot2::theme_minimal() +
  ggplot2::theme(
    plot.background = ggplot2::element_rect(colour = NA,
      fill = if (figure_mode == "light") "#F8FBFE" else "#101D2D"),
    text = ggplot2::element_text(
      colour = if (figure_mode == "light") "#1D3448" else "#F0F6FB"),
    axis.text = ggplot2::element_text(
      colour = if (figure_mode == "light") "#4D6778" else "#B8CEDD"))

Moran 散点图的读图顺序

顺序	怎样读	本图中的含义
1 · 横轴	读取该格的标准化报价。	$z=0$ 对应362个参检格的中位报价均值；$z>0$ 为高于这一均值。
2 · 纵轴	读取邻居标准化报价的行加权平均。	$Wz>0$ 表示邻域平均报价高于同一参照，不要求每个邻居都偏高。
3 · 象限	对照本格与邻域的高低。	HH122格、LH56格、LL165格、HL19格；这些是正负组合，还没有筛局部显著性。
4 · 全局线与检验	查看整体方向，再看随机置换的证据。	过原点线向右上倾斜，斜率为0.620857；9,999次单侧置换的$p=0.0001$。

本图的过原点拟合方程为

\[ \widehat{Wz}=0.6209z. \]

当某格的中位报价高于样本均值1个标准差时，直线上对应的邻域加权平均约高出同一均值0.621个标准差。散点仍可偏离直线，整体相似性与局部反差可以同时存在。

进阶 · 为什么斜率等于 Moran’s I

纵轴$Wz$由已经标准化的邻居值计算，没有再单独标准化。全局统计量为

\[ I=\frac{n}{S_0}\frac{z^\mathsf{T}Wz}{z^\mathsf{T}z}. \]

本例排除了无邻居格，在同一362格上重新标准化，且每行权重和为1，因此$n=S_0=362$，上式等于过原点拟合线的斜率。更换权重规范、样本或孤立格处理后，需要重新核对这一关系。

图中没有普通独立观测 OLS 的置信带。空间关联的显著性来自固定图结构、重排报价的置换检验；$I$不作为方差解释比例。参见 GeoDa全局空间自相关教程。

局部空间自相关

2.6 LISA：局部关联在哪里

全局 Moran’s I 概括整体关联，局部 Moran 检验逐格考察本格与邻域的组合。本节沿用相同的 362 个参检格、Queen 权重、$z$ 与 $Wz$。高低以这批格的整体均值为参照，与 Fisher 地图上的颜色等级分别解读。

每格进行 9,999次双侧条件置换：固定该格的报价，从其他格抽取邻居值，同一次模拟的邻域内不重复抽取同一个格。本例默认用原始模拟$p\leq0.05$识别局部类别，同时保留 BH、BY、Holm 与 Bonferroni 调整作比较。

## Warning in parallel_setup(iseed): no cluster in ClusterOption, parallel set to
## no

原始$p\leq0.05$下，112格通过局部检验：HH44格、LH2格、LL64格、HL2格；其余 250格未通过。地图中的不足5ID 格及3个无 Queen 邻居格未参加本轮检验，统一作为背景灰显示，不能当作“局部不显著”或零报价。

同一批362格的局部Moran地图与Moran散点：原始双侧模拟p不超过0.05时设色，共112格；HH44、LH2、LL64、HL2。250个参检格未通过阈值；另291个未检验格只作背景灰显示。

p值处理	通过阈值	HH	LH	LL	HL	未通过
none	112	44	2	64	2	250
BH	52	24	0	28	0	310
BY	11	5	0	6	0	351
holm	0	0	0	0	0	362
bonferroni	0	0	0	0	0	362

BH 调整后保留 52格（HH24、LL28），BY 调整后保留 11格（HH5、LL6）；原始p下的4个 LH/HL 未在这两种调整后保留。局部数量与范围对检验口径敏感，不能把原始112格称为经多重检验确认的112格。

BH 针对错误发现比例的期望，其保证需要独立或特定正相关条件；BY 允许更一般的依赖结构。Holm 和 Bonferroni 控制整组至少一次误报的概率。本例9,999次双侧模拟的最小p为0.0002，高于$0.05/362$，因此后两种方法的0格还受到模拟分辨率限制，不能据此断言没有局部空间关联。

查看并复制局部空间自相关 R 代码

# 6. 局部条件置换：沿用同一z和W，固定单核与随机种子。
set.coresOption(1L)
set.mcOption(FALSE)
set.ClusterOption(NULL)
RNGkind("Mersenne-Twister", "Inversion", "Rejection")
set.seed(20261025)
local_permutation <- localmoran_perm(z, w, nsim = 9999,
  alternative = "two.sided", iseed = 20261025,
  no_repeat_in_row = TRUE, mlvar = TRUE, sample_Ei = TRUE,
  adjust.x = FALSE, zero.policy = FALSE)

# 只提取双侧模拟秩p；不用正态近似列或folded尾概率替代。
p_raw <- as.numeric(local_permutation[, "Pr(z != E(Ii)) Sim"])
quadrant <- case_when(
  z > 0 & wz > 0 ~ "HH",
  z < 0 & wz > 0 ~ "LH",
  z < 0 & wz < 0 ~ "LL",
  z > 0 & wz < 0 ~ "HL",
  TRUE ~ "Axis")
stopifnot(length(p_raw) == 362L, all(is.finite(p_raw)),
          !any(quadrant == "Axis"))
analysis_hex$lisa_cluster <- ifelse(p_raw <= 0.05, quadrant, "NS")
stopifnot(sum(analysis_hex$lisa_cluster == "HH") == 44L,
          sum(analysis_hex$lisa_cluster == "LH") == 2L,
          sum(analysis_hex$lisa_cluster == "LL") == 64L,
          sum(analysis_hex$lisa_cluster == "HL") == 2L,
          sum(analysis_hex$lisa_cluster == "NS") == 250L)

# 7. 分类比较：每种方法一次性处理全部362个p值，不重新模拟。
adjustment_methods <- c("none", "BH", "BY", "holm", "bonferroni")
p_adjusted <- vapply(adjustment_methods,
  function(method) p.adjust(p_raw, method = method), numeric(length(p_raw)))
lisa_class_by_method <- apply(p_adjusted, 2,
  function(p) ifelse(p <= 0.05, quadrant, "NS"))
lisa_method_counts <- bind_rows(lapply(adjustment_methods, function(method) {
  classes <- lisa_class_by_method[, method]
  data.frame(method = method, significant = sum(classes != "NS"),
             HH = sum(classes == "HH"), LH = sum(classes == "LH"),
             LL = sum(classes == "LL"), HL = sum(classes == "HL"),
             NS = sum(classes == "NS"))
}))
stopifnot(identical(lisa_method_counts$significant, c(112L, 52L, 11L, 0L, 0L)))
lisa_method_counts %>%
  knitr::kable(col.names = c("p值处理", "通过阈值", "HH", "LH", "LL", "HL", "未通过"))

制图核心。 按hex_id把当前 LISA 类别连回全部完整格；未参加检验的291格保留背景灰，参检但未通过的250格标为 NS。散点与地图使用同一组类别。

# 8. 当前分类以原始p为例；切换校正方法只替换lisa_cluster。
figure_mode <- "light"
lisa_colours <- list(
  light = c(HH = "#B83F4D", LH = "#087D69", LL = "#2867A4",
            HL = "#7B4AA7", NS = "#959EAA", Untested = "#D9E1E8"),
  dark = c(HH = "#FF9A94", LH = "#65D7BF", LL = "#8ABCF6",
           HL = "#C4A4F4", NS = "#96A7B8", Untested = "#354A60"))
lisa_lookup <- analysis_hex %>% st_drop_geometry() %>%
  select(hex_id, lisa_cluster)
lisa_plot <- hex_all %>% left_join(lisa_lookup, by = "hex_id") %>%
  mutate(lisa_cluster = coalesce(lisa_cluster, "Untested")) %>%
  st_intersection(st_union(main_island_2263)) %>% st_transform(4326)
ggplot2::ggplot(lisa_plot) +
  ggplot2::geom_sf(ggplot2::aes(fill = lisa_cluster), colour = NA) +
  ggplot2::scale_fill_manual(values = lisa_colours[[figure_mode]],
    breaks = c("HH", "LH", "LL", "HL", "NS"),
    drop = FALSE, name = "局部类别") +
  ggplot2::coord_sf(datum = NA) + ggplot2::theme_void() +
  ggplot2::theme(
    plot.background = ggplot2::element_rect(colour = NA,
      fill = if (figure_mode == "light") "#F8FBFE" else "#101D2D"),
    text = ggplot2::element_text(
      colour = if (figure_mode == "light") "#1D3448" else "#F0F6FB"))

LISA 地图的读图顺序

区分未检验与未通过。 背景灰含288个不足5ID 格及3个无邻居格；“NS”只指362格中完成局部检验但未过门槛的250格。
核对高低参照。 HH 是本格中位报价与邻域加权平均都高于362格均值；LL 是两者都低于该均值，不要求每个邻居同为高或低。
定位同类组合与局部反差。 原始p下 HH 主要位于主岛中南段，LL 主要位于北段；LH 与 HL 各2格表示本格和邻域的高低相异。
比较检验口径。 对照原始p、BH 和 BY，观察局部证据保留多少，再讨论位置和范围。

这里的“热点／冷点”指报价的高—高／低—低局部组合。它们不表示人流热度、入住率或出租量；本节也未另行进行 Getis–Ord 热点检验。解释实际成交、房型与社区因素的作用，需要额外数据与研究设计。

北部报价较低的网格与中南段报价较高的网格提供了可继续研究的空间线索。社区历史、住房与社会经济资料可用于提出假设，但当前 LISA 本身没有识别价格差异的原因。

进阶 · 北部社区、社会经济不利与低报价片区

从社区历史认识北部。 Harlem 是非裔美国人的重要居住与文化中心，20 世纪初的人口迁入及 Harlem Renaissance 塑造了其文化地位；此后的就业变化、建筑修复与新建持续改变着社区，见美国国家公园管理局的 Harlem 历史介绍。Washington Heights 与 Inwood 则具有鲜明的移民社区背景。纽约市规划局《The Newest New Yorkers 2026》采用 2019—2023 年 ACS 数据，估计两地外出生人口分别约占 44% 和 41%，见报告第三章。这些资料展示了北部社区的历史与人口差异。

社会经济不利怎样进入分析？ 社会经济不利（socioeconomic deprivation）涉及收入、就业、教育与住房条件等方面的资源不足。可先用贫困率和住房负担等具体指标描述社区，再检查它们与低报价片区的关系。

NYU Furman Center 社区网页中标为 2024 年的指标提供了以下背景：

社区统计区	贫困率	严重租金负担家庭占租房家庭比例
Central Harlem，MN10	26.7%	29.3%
Washington Heights/Inwood，MN12	19.2%	30.2%

同一来源给出的全市贫困率为 18.0%；严重租金负担指租金支出超过家庭收入的 50%。北部各社区的社会经济条件存在差异，低住宿报价与居民较重的住房负担也可能同时出现：前者衡量公开住宿价格，后者还取决于居民收入。

由此可以提出的假设是： 在统一入住时长，并考虑入住日期、房型、可住人数与区位后，社区社会经济条件与六边形报价之间是否仍有关联？检验时需将社区资料与六边形进行空间匹配，核对统计年份与样本覆盖。这里的社区资料按 Community District 或报告中的社区单元汇总，未与本例 LL 格逐格叠加；它们提供解释线索，当前 LISA 识别的对象仍是报价的局部关联。

为什么整体显著，多数格却未通过局部检验？

全局 Moran’s I 把全部362个参检格的相邻关系汇总为一个统计量；局部检验分别考察各格，门槛与多重检验处理也会影响识别范围。本例默认原始p下有250格未通过，BH 后为310格、BY 后为351格。未通过表示证据未达到所设门槛，不能据此认定没有空间关联；另291格没有参加检验，不应并入“NS”。参见 GeoDa局部空间自相关教程。

3.0 专题研究 1：回顾与反思

本例完成了获取数据、预处理、空间聚合、全局 Moran’s I 与 LISA 局部分析的基本 ESDA 流程。回顾分为三个层次：处理后的一条观测代表什么，当前结果说明什么，以及哪些研究选择需要进一步检验。

3.1 研究回顾：问题与回答

基本流程	已完成的操作	得到的分析对象或证据
数据获取与初探	读取 3、6、9 月三个快照，限定主岛，检查报价状态并绘制原价／对数分布。	研究范围内的房源—批次记录。
数据预处理	保留有效正价，核查 ID 与批次唯一性；按 ID 求均价，再筛均价的 P5–P95。	每个挂牌房源 ID 一个代表报价与位置。
六边形密铺与聚合	每个 ID 唯一归格；至少含 5 个不同 ID 的格计算中位报价。	固定面积的区域观测及各格样本量。
全局 Moran’s I	排除无邻居格，重建 Queen 权重和标准化值，绘制散点图并作全局置换。	整体空间关联的方向、量级及检验结果。
LISA	沿用同一批参检格与权重，作局部条件置换，比较原始与调整后的 p 值分类。	局部集聚、空间离群及其位置。

记录数为何逐步变化？ 主岛 43,589 条原始记录 → 27,736 条有效正价记录 → 11,853 个 ID 均价 → P5–P95 后 10,669 个 ID → 365 个至少含 5 个不同 ID 的达标六边形 → 362 个有 Queen 邻居的参检格。

这条流程包含筛选，也包含分析单位的转换。按 ID 汇总将同一挂牌房源的跨期观测合成一个代表报价；六边形聚合再将房源报价概括为区域中位数。计算比例或均值时，应说明分子、分母及其单位，不能将记录数、房源数与六边形数当作同一总量替换。坐标相同的不同 ID 在整个流程中仍分别保留。

研究问题	当前结果	怎样解释？
报价如何分布？	有效报价右偏、长尾；六边形图中，北段较多格报价偏低，中南段有更多报价较高的格。	描述的是当前样本的挂牌报价；房型、入住条件和缺失分布仍可能影响地区差异。
相近报价是否倾向于相邻？	362 格的 $I=0.621$，9,999 次单侧置换的 $p=0.0001$。	当前尺度与权重下呈全局正空间自相关；I 不表示价格差异的解释比例。
局部集聚在哪里？	原始 $p\leq0.05$ 下，HH 44 格主要在中南段，LL 64 格主要在北段；LH、HL 各 2 格。	HH／LL 是本格与邻域相对于整体均值的同向组合；LH／HL 提示局部反差。
局部判断依赖检验口径吗？	原始 p 下有 112 个显著格，BH 后为 52 格，BY 后为 11 格。	显著范围会随多重检验处理变化；仅凭主图不能断言局部边界稳定。

全局结果概括整体关联，局部结果将高低值组合落实到位置。二者共同为房型、区位与社区条件研究提供线索；价格成因仍需补充资料和解释设计。下文先回顾空间聚合，再讨论影响这些判断的研究设置。

3.2 知识回顾：从点到区域

为什么把点汇总到网格？

原始点适合查看位置与点间差异。本例关注相邻局部区域的典型报价是否相似，因此将每个房源 ID 的代表点归入固定面积的格，以格内中位价比较区域。点数据也可直接使用距离或 $k$ 最近邻权重，方法见 GeoDa 点距离权重教程。

六边形分箱 （hexagonal binning / hexbin）：以无空隙、不重叠的六边形划分平面，将每个点唯一归格，再计算数量、均值或中位数。密铺规定单元怎样覆盖空间，分箱规定点怎样归组和汇总。

420 个相同的模拟点及其属性；方格与六边形完整单元面积相同，颜色表示点属性或格内中位数。

同一批模拟点：原始点、等面积方格与六边形分箱

读图。 同一批模拟点分别汇总到等面积方格与六边形，整体高低值位置相近，边界附近的混合和局部细节有所不同。聚合得到区域代表值，也压缩了格内差异。

为什么选六边形？方格可以吗？

比较项	正六边形	正方形
完整密铺内部的邻居	6 个共边邻居，中心距离相等。	4 个共边邻居；加入角邻居后为 8 个，角邻居距离是边邻居的 $\sqrt{2}$ 倍。
形状与数据衔接	同面积时更紧凑，适合比较六个方向的邻域，仍有固定方向。	便于与栅格像元、行列索引衔接，也可旋转。
本例需要核对	少样本格退出后，邻居可能不足 6 个；实测 Queen 与 Rook 列表相同。	Queen 与 Rook 通常不同，需随几何重新定义。

几何解释见 Esri：Why hexagons?。两种形状均需检查面积、起点和边缘。

等面积比较的参数。 当 $A=0.10\text{ km}^2$ 时，方格边长 $\sqrt{A}\approx316.23$ m；正六边形边长 $\sqrt{2A/(3\sqrt{3})}\approx196.19$ m，对边距离约 339.81 m。st_make_grid() 的六边形 cellsize 指对边距离；只切换 square 而保留原 cellsize 会改变面积。长度须换成投影单位，本例 EPSG:2263 使用美国测量英尺。见 st_make_grid() 文档。

选学进阶 · Uber H3 的层级网格与索引压缩

Uber 的 H3 是以六边形为主的全球层级空间索引（hierarchical geospatial index）。单元有唯一编号，可查询父格、子格和邻居；分辨率 0–15，数值越大，单元越小。

按密度调整大小。 可先在同一细层级编码所有点，再设置点数阈值：过密区域继续细分，稀疏区域沿父索引合并。最终每点只归入一个组，各组重新计算统计量。规则需另行设定；H3 已答讨论与自适应网格交互示例展示了实现思路。

索引压缩。 compactCells() 将齐全的同层子格集合递归换成父格；它读取索引集合，不读取房源点数。下图展示这一操作。

同一组加州 H3 细格的压缩前后对照。左图统一细分辨率，右图以父格表示完整子格集合；颜色表示 H3 分辨率。

加州 H3 索引：统一细分辨率与压缩后的混合分辨率表示

加州的 10,652 个分辨率 6 索引压缩为 914 个混合层级索引，数量减少约 91.4%；展开回分辨率 6 后，索引集合完全恢复。内部容易凑齐完整子格组，边缘较多细格保留。颜色表示分辨率，格的大小由子格集合是否齐全决定，与人口或房源密度无关。

图按 Uber 原文与 H3 官方索引说明的思路独立计算；行政轮廓采用 U.S. Census Bureau 制图边界。

使用边界。 父子关系在索引上精确，原生几何仅近似嵌套，多层级绘图可出现间隙或重叠，见边界交互示例。每个分辨率还包含 12 个五边形，同级面积并非处处相同；尺度表中的分辨率 9 平均面积约 0.105 km²，只可作为尺度参照。

应用阅读。 Uber 原文介绍出行供需；Volkswagen／NVIDIA 案例以每批约 2,500 条位置记录为计算约束，逐级细分过密组。若用混合层级格比较报价，还需记录实际覆盖与点数，重新定义邻居。本例仍保留固定面积的 sf 平面六边形。

3.3 研究设置的不足与改进

数据检查可以发现空值、解析失败和重复批次记录，但不会自动决定研究范围、价格变量或空间尺度。本例采用的一组设定使分析得以实施，仍存在其他合理选择。敏感性分析通过有计划地改变这些设定，检查主要发现的适用范围及发生变化的位置。

环节	当前设定	可以怎样比较？
范围与时点	曼哈顿主岛，3、6、9 月三个快照。	分批次或按房型比较，检查样本构成、缺失与覆盖差异。
价格与代表值	ID 算术均价；格内中位数；原美元统计变量。	分别比较 ID 层代表值、格层代表值和统计变量变换，避免同时改动。
价格筛选	在 ID 均价层保留 P5–P95。	对比不截尾、P1–P99，并记录被排除 ID 及有效覆盖。
空间单元与支持	0.10 km² 六边形，至少 5 个不同 ID。	分开比较面积、起点、等面积形状与 ID 数量门槛。
邻居与权重	Queen、行标准化，排除孤立格后重建。	比较有依据的距离或最近邻规则，检查连通部分和长连线。
显示与统计判断	Fisher 七级设色，局部原始 $p\leq0.05$。	分级数只改显示；p 调整方法和显著性阈值分别比较。

样本与报价口径

当前使用平台每夜挂牌报价。实际入住日期、时长、房型、可住人数和折扣仍有差异，三个批次还存在上下架及价格缺失。可先按批次、房型检查缺失比例和报价分布，再比较相近住宿产品。跨期均价概括已观测时点，没有消除季节变化，也不能代替成交价或入住率。

按房源 ID 汇总后，每个 ID 对分位筛选和格内中位数各贡献一个代表报价；不同 ID 即使位于同一公开坐标也分别保留。代表位置取该 ID 最新一条有效正价记录的公开坐标，坐标跨期变化时，归格可能受此规则影响。未出现的批次保留为未观测状态。房源 ID 标识公开挂牌房源，不能据此证明是不同建筑、经过核实的独立物理住处或实际成交。

代表值与代表位置可以分别检查。 当前存在两次不同层级的汇总：ID 层概括跨期报价，格层概括区域内的房源报价。比较 ID 均值与中位数时，从模块 7 开始重算分位筛选及下游结果；本例每个 ID 至多三条报价，只有一条或两条时，R 默认中位数与均值相同。比较格内中位数与均值时，则固定入选 ID 和位置，从 2.4 节的格内汇总代码开始重算，以检查高价房源对区域代表值的影响。

位置规则可先比较最新与最早有效正价记录的坐标：固定 ID 报价与分位入选集合，只更换代表位置，统计多少 ID 改格，再重做格内汇总、邻居和检验。经纬度均值可能落在未观测的位置，采用前须另作说明。

价格筛选与空间覆盖

P5–P95 聚焦房源均价的主体范围，也改变格内组成及达到 ID 门槛的区域。当前在 11,853 个 ID 均价中保留 10,669 个，显示范围为 63.96–1,000 美元／夜，两端各排除 592 个 ID。区间外报价仍可能有效；分位阈值应与研究对象相联系。

价格范围的对照应同时检查报价与覆盖。 可固定 11,853 个 ID 的代表均价与位置，比较不截尾、P1–P99 和当前 P5–P95。每次从模块 8 重算分位边界与入选 ID，再按同一网格重新汇总，检查达标格、孤立格与空间检验。

比较层次	具体做法	回答什么？
各方案的完整结果	每个价格范围都运行完整流程，报告入选 ID、达标格、参检格、I 与局部类别。	采用这一筛选范围后，实际得到怎样的覆盖和结果？
固定网格下的共同格	另取各方案共同参检的格，在相同顺序下重新检查邻居、构建 W 并计算各方案的 z 与 Wz。	限定在同一空间覆盖后，报价差异仍会怎样影响结果？

共同格分析不能直接截取旧 W 或沿用旧 z：邻居减少会改变权重，参检集合改变也会改变标准化参照。若共同集合出现孤立格，还须按统一规则处理并重新核对。比较局部类别时，分别记录“退出检验”与“仍参检但由 HH／LL 变为 NS”，避免把样本退出读成集聚消失。

上述价格范围的系统对照尚未执行。当前 365 个达标格、362 个参检格和 I 只对应主方案，不能沿用到新筛选样本。

空间聚合与尺度选择：当前模式是否稳定？

本例采用 0.10 km² 六边形、每格至少 5 个不同挂牌房源 ID。同坐标不同 ID 分别计数，同一 ID 跨批次只计一次；该门槛不能解释为五处不同物理住处或五次成交。这项选择同时影响局部细节、格内数据量及有效覆盖，需要结合敏感性分析讨论。

15 分钟城市：生活圈与计算单元

15 分钟城市强调居民通过约 15 分钟步行或骑行获得日常服务。Moreno 等（2021）从密度、邻近、多样性与数字化讨论这一规划框架。它为本例提供邻里尺度的观察视角。

Bruno 等（2024）的作者公开稿方法部分采用边长 200 m的六边形，再用 OSRM 计算到设施的路网时间。其计算单元与本例约 196 m 的边长接近，为细尺度选择提供了研究实现参照。

本例的 0.10 km² 承载报价汇总；15 分钟服务范围取决于路网、出行方式和设施位置。 时间目标不能直接推导固定面积，面积是否适合当前数据仍需检验。

先看一次聚合隐藏了什么。 当前两个真实格的中位报价相近，内部报价却很不同：

检查项	格 662	格 381
中位报价（美元／夜）	232.74	239.45
不同房源 ID 数 n	46	22
中间一半房源均价所在区间 Q1–Q3（美元／夜）	216.48–259.74	86.74–466.98

IQR 是上述区间的宽度 $Q_3-Q_1$，两格分别为 43.26 与 380.24 美元。它衡量格内离散程度；n 表示汇总所依据的不同房源数量，两者共同补充中位价图。

查看当前网格的诊断图 · 两格分布与全岛 n／IQR

662格与381格的中位报价分别为232.74和239.45美元，房源ID数为46和22，四分位距为43.26和380.24美元。

两个中位报价相近的真实六边形：房源均价与箱线图

两图纵轴一致，每点代表一个房源 ID 的已观测报价均值。662、381 两格作为固定教学例，说明相近的区域中位数可能对应差异很大的格内分布；它们不能代表全岛所有格，也不能据此给出中位数的精度或置信区间。

当前365个达标格；左图是不同房源ID数，右图是房源均价的IQR，分别以五级Fisher分级，不足5个ID的格以固定中性色显示。

当前0.10 km²网格的样本支持：不同房源ID数与格内报价IQR（分别采用五级Fisher）

365 个达标格的房源 ID 数为 5–299，中位数 20；格内房源均价 IQR 的中位数为 148.83 美元／夜。两图分别以五级 Fisher 显示数量与 IQR，再与中位报价、LISA对读。统计基于每个 ID 的代表均价，分位点用 type = 7；房源数与 IQR 不能直接解释为中位数的精度，至少 5 个 ID 只规定纳入资格。定义见 R IQR 文档，重抽样推断可选读 Bootstrap 教程。

当前数量与 IQR 的计算代码

# 当前样本支持：从同一份ID归格表计算数量和格内报价离散程度。
# 诊断使用全部365个达标格，区别于Moran/LISA的362个参检格。
hex_support_current <- listing_members %>%
  group_by(hex_id) %>%
  summarise(n_id = n_distinct(id),
            median_usd = median(listing_mean_usd),
            q1 = quantile(listing_mean_usd, 0.25, type = 7),
            q3 = quantile(listing_mean_usd, 0.75, type = 7),
            iqr_usd = IQR(listing_mean_usd, type = 7), .groups = "drop") %>%
  filter(n_id >= 5L)
stopifnot(nrow(hex_support_current) == 365L,
          min(hex_support_current$n_id) == 5L,
          median(hex_support_current$n_id) == 20,
          max(hex_support_current$n_id) == 299L,
          abs(median(hex_support_current$iqr_usd) - 148.83) < 1e-8)

再改变空间划分，检验 MAUP。 上述图检查当前固定网格；改变单元面积或边界后重算，才是在比较尺度效应与分区效应。可回看导论及 NCRM 的 MAUP 教学案例。

起步方案：先做三组有区分度的检查（尚未运行）。

改变什么？	怎样设置？	重点比较什么？
面积：尺度效应	0.05、0.10、0.20 km²；固定价格口径、网格方向与基准起点，并记录参数。	区域中位价、有效覆盖及 HH／LL 片区的位置与连片程度。
起点：分区效应	固定 0.10 km²，先从格心坐标核对 x、y 方向各自的重复间距，再分别沿一个方向平移该方向间距的一半；扩展边缘候选格以覆盖整个主岛。	同一批点重新归格后，局部边界和中位价如何变化，是否有漏归点。
房源 ID 门槛：样本支持	固定网格，比较至少 3、5、10 个不同房源 ID。	哪些区域失去覆盖，统计量与局部类型是否随之改变。

方格对照可沿用等面积参数另设一组。房源 ID 数量门槛改变的是样本纳入，需与改变空间单元的 MAUP 检查分开记录。

网格几何改变后，怎样比较位置？ 面积、起点或形状一变，相同 hex_id 就不再代表同一地面单元。应先并列显示各方案，再在主岛的共同可检验范围内比较 HH／LL 的地面覆盖、交叠比例与主要片区位置。为比较范围而叠加得到的碎片，只用于地图对照，不作为新的独立观测重新计算 Moran’s I。

怎样组织一次敏感性分析？

预先列出候选。 各组均以主方案为基准，其余设置保持主方案。可依次比较价格范围、ID 门槛和面积，每组只改一项；比较门槛时固定 0.10 km²，比较面积时固定至少 5 个 ID。
从受影响的模块开始重算。 每组参数对应一个 scenario_id，分别保存输入、参数和结果。上游样本或网格改变后，不复用旧的格属性、W、z 或置换结果。主方案中的数量断言须改为该方案的实际核对，保留每个 ID 唯一归格、单位和观测顺序检查。
同时报告完整覆盖与可比范围。 固定几何时，补充共同格分析；几何改变时，比较地面范围。对照图沿用同一底图和 LISA 色盘；同一价格变量尽量使用固定参考断点，并完整显示超出参考范围的值。
记录改变了什么。 除 I 和显著格数外，还检查哪些区域退出、哪些局部类别转换、哪些片区在各方案中持续出现。保留无显著格或覆盖不足的候选，不只呈现符合预期的地图。

每个方案可用下表整理结果。价格、网格和邻居方案仍是后续研究设计，当前未据此声称结果稳健。

报告项目	至少记录什么？
参数与输入	`scenario_id`、改变项、分位点、入选 ID 数及代表值／位置规则。
空间覆盖与支持	达标格、参检格、未纳入 ID；格内 ID 数及 IQR；孤立格、连通部分与邻居数分布。
全局关联	I、n 与权重和；置换次数、方向、种子和 p 值；适用时另列共同格结果。
局部关联	HH／LL／LH／HL／NS 数量及检验口径；另列未检验格、类别转换与片区交叠。

敏感性分析用于说明发现依赖哪些条件。参数应依据研究尺度、报价含义、位置误差和样本支持选择，不以 I 更大、p 更小或显著格更多来决定主方案。

对北部 LL 的解释还缺哪些资料？

北部 LL 提供了局部低报价集聚的线索。可先按房型重新执行房源 ID 汇总、价格筛选和网格聚合，检查相似住宿产品是否仍有南北差异。不同房型的覆盖可能不同，需同时报告各组有效格；同一公开坐标下的不同 ID 继续分别保留。

北部社区资料提供历史、人口及住房负担背景。若要检验社区条件能否解释报价差异，需要匹配边界与年份，再纳入房型、入住日期和区位信息。当前图件完成的是空间分异描述。

3.4 地理学方法拓展与阅读

聚合之后，怎样定义“周围”？

权重方式	表达什么关系？	检查重点
接壤邻接	Queen 包括共边或共顶点；Rook 要求共边。	实际邻居列表、边缘单元。
距离阈值	指定距离内纳为邻居。	距离依据、单位、无邻居格。
k 最近邻	每格选最近的 k 个其他格。	稀疏区长连线与关系是否对称。

交互演示：从邻接到空间滞后。 切换网格形状与邻居规则，再依次查看邻接、权重和取值。点击计算区可展开完整矩阵或逐项计算；图中报价为教学设定。

“看取值”计算的是报价 $x$ 的空间滞后 $Wx$；将输入换成标准化报价 $z$，即可对应 Moran 散点图的纵轴 $Wz$。完整矩阵同时展示边缘格，便于比较各行邻居数与权重。

365 个达标格的 Queen 图有 5 个连通部分，其中 3 格孤立；去掉这 3 格后，362 个参检格仍分为 2 个连通部分，分别为 357 格和 5 格。连通部分由邻接路径划分，与 HH／LL 类型不同。

邻居方案可先比较格心距离 500／1,000 m，或 k=4／6／8；距离须换成 EPSG:2263 的美国测量英尺。k 最近邻还需说明保留有向关系，还是按“任一方选中”或“互相选中”对称化。当前规则六边形的 Queen 与 Rook 已核对为相同，因此切换二者只能核对实现，距离或 k 最近邻才会改变本例的邻居定义。

比较权重时，固定房源、网格和格属性，从 2.5 节的邻居构建代码开始重算，再重算标准化值与全局／局部检验。若不同规则产生不同参检集合，先报告各自完整结果，再补充满足各规则的共同格对照；同时检查孤立格、连通部分、邻居数和稀疏区长连线。

邻接图如何变成矩阵？

设有四个依次相邻的分析单元 $1\leftrightarrow2\leftrightarrow3\leftrightarrow4$，各行、各列均按 1、2、3、4 排列。这个教学例子演示权重计算。

二元邻接矩阵 $B$ 记录邻居关系，相邻为 1、不相邻为 0，对角线为 0：

\[ B=\begin{pmatrix} 0&1&0&0\\ 1&0&1&0\\ 0&1&0&1\\ 0&0&1&0 \end{pmatrix}. \]

第 $i$ 行有 $k_i=\sum_j b_{ij}$ 个邻居。对有邻居的行，行标准化令 $w_{ij}=b_{ij}/k_i$，得到

\[ W=\begin{pmatrix} 0&1&0&0\\ 1/2&0&1/2&0\\ 0&1/2&0&1/2\\ 0&0&1&0 \end{pmatrix}. \]

假定四格报价为 $x=(100,200,300,400)^\mathsf{T}$ 美元／夜，则 $Wx=(200,200,300,300)^\mathsf{T}$。例如第 2 格周围是第 1、3 格，两者的平均为 $(100+300)/2=200$。将 $x$ 换成标准化报价 $z$，得到的就是 Moran 散点图纵轴 $Wz$。

这个矩阵还说明：邻居关系可以对称，行标准化后的数值权重却未必对称。第 1 格只有一个邻居，给第 2 格的权重为 1；第 2 格有两个邻居，给第 1 格的权重为 $1/2$。因此，检查邻接与检查权重数值要分开进行。矩阵行列的 ID 顺序必须与价格向量一致。

选读 Paula Moraga 的空间邻近矩阵教材；定义和实现见 GeoDa 接壤权重教程及 spdep 的 nb2listw() 文档。

zero.policy = TRUE 可约定无邻居格的空间滞后为零，但这不是观测到的邻域均价。本例先排除 3 个孤立格，再重建 362 格的权重，并用 zero.policy = FALSE 检查每格都有邻居。全局置换在全部参检格间重排；若改为连通部分内部重排，原假设也随之改变。见 spdep 连通部分说明。

数值变换与显示分级会怎样影响结果？

固定当前 362 格、观测顺序和 Queen 权重，原美元中位报价 $x$ 与其线性标准化值 $z=(x-\bar{x})/s_x$ 得到相同的 $I=0.621$。标准化统一平移和缩放数值，不改变这一统计量。

对格中位报价 $x$ 取对数，会改变格间差距，因而可能改变 I 和局部类别。可固定当前格集合与 W，从统计输入开始比较 $x$ 与 $\log(x)$，分别重算 z、Wz 及置换检验，观察主要片区在绝对金额与比例差异两种视角下是否持续出现。若改为先对每条报价取对数再按 ID 汇总，则须从房源代表值开始重算，它与“对格中位数取对数”回答的问题不同。上述变量对照尚未执行。

Fisher 分级与仅用于绘图的对数色标只改变显示；检验继续输入原美元中位报价时，统计结果不受设色等级影响。价格范围筛选则改变纳入的 ID 和空间覆盖，见价格方案比较。

局部检验与区域解释各有什么边界？

362 个参检格使用同一批局部双侧模拟 $p$ 值。原始 $p\leq0.05$ 得到 112 格；对全部 362 个 $p$ 值一起作 BH 调整后为 52 格，BY 调整后为 11 格。局部数量和边界随检验口径变化，不能仅凭显著格较多选择方法。BH 关注被判为显著结果中错误发现比例的期望，即错误发现率（false discovery rate, FDR），保证依赖有效 $p$ 值及相应依赖条件；BY 允许更一般的依赖结构。方法见 p.adjust()。

只改变 p 的处理时，无需重跑置换。 固定参检格、变量与 W 后，可复用保存的全部原始 p 值，分别比较调整方法及 0.05／0.01 门槛；调整时纳入该方案所有参检格，不能只调整某个象限或已显著的格。此时散点位置和全局线保持不变，不对显著子集重新拟合。变量、样本或 W 一旦改变，则须重新进行相应置换。

生态谬误。 即使将来发现报价较低的网格集中在贫困率较高的社区，也只能先描述区域关联；判断低收入居民或房东是否提供更低报价，需要相应个体资料。可阅读 Robinson（1950）。

解释模型之前，还可检查两件事。 一是按格心距离查看空间相关图，同时列出距离区间的格对数量，观察空间关联随距离的变化；南北方向的趋势需另用纬度—报价关系或方向性检查。二是考虑公开位置约 0–150 m 的匿名偏移，按明确规则做位置扰动，记录归格和覆盖变化。两项尚未运行；扰动是情景比较，不能还原真实位置。若后续拟合模型，剩余关联应使用适合残差的检验，如 lm.morantest()。

拓展阅读：从概念到研究实施

15 分钟城市的概念。 Moreno 等（2021），Smart Cities，可读出版社全文第 5–6 节：时间预算怎样联系邻近、多样性与出行方式？
细尺度单元与可达性。 Bruno 等（2024），Nature Cities及作者公开稿：先看图 1，再读方法 3.2–3.3，区分网格、路网时间及人口加权。
从研究到规划。 C40：Green and Thriving Neighbourhoods：邻近规划怎样落实到街区服务与公共空间？
空间统计的后续问题。 Spatial Data Science讨论趋势、遗漏变量与自相关；GeoDa 教程连接全局、局部指标与地图。
Airbnb 与城市变化。 Wachsmuth 与 Weisler（2018）联系纽约租差与绅士化；Barron、Kung 与 Proserpio（2021）帮助理解研究住房价格和租金影响还需要哪些数据与识别设计。

4.0 专题研究 2

4.0 专题研究 2: 社交媒体点

专题研究 2: 上海六区微博记录空间分析

分析导航

微博记录热点与交通设施的空间对比

从记录分布到高密度区，再到设施距离比较

本节选取黄浦、徐汇、长宁、静安、杨浦与闵行六区，使用微博记录与公共交通 POI，依次完成范围筛选、核密度估计、高密度区提取和设施邻近性比较。

点击下方卡片，快速跳转至相应模块：

课题与目标

4.1 研究课题与分析目标

示例课题

上海六区微博记录热点与公共交通设施邻近性比较

公共交通设施提供出行接入点，地理标记微博留下部分平台活动的位置线索。本课题结合两类点数据，比较六区微博记录高密度区内外的设施距离，建立从原始点、KDE密度表面到AOI内外比较的分析流程。

累计热点是否来自日常活动、短期事件或重复发帖，以及相近的设施距离是否对应相近的出行条件，需要时间分段、路网与人群资料进一步检验。

三个递进问题

热点在哪里？ 2020–2025年汇总的微博记录，在各区主要集中于哪些位置？
设施有多近？ 各区高密度区内外的记录位置，到最近公交站、地铁站的距离分布有何差异？
比较是否稳定？ 改变带宽、密度阈值或时间窗口后，热点位置和距离差异是否保持？

前两问对应本节的主要流程，第三问用于后续验证。AOI内设施距离更短、更长，或差异不明确，均按实际结果报告。

从记录集聚到城市活力

记录集聚

描述样本记录在空间上的集中程度。一条记录对应一次被纳入样本的发布或签到行为，同一用户可以留下多条记录。

城市活力

涉及活动的强度、持续性、多样性及时间变化，不同研究的定义和指标有所区别（Liu et al., 2023）。讨论活力还需结合独立用户、活动类型、人口或出行资料。

累计热点的持续性需另行检验。 六年汇总图叠加了日常活动与短时事件。按月、季节或年份重复分析，才能判断哪些高值区持续出现。

研究延伸：怎样从设施邻近性走向可达性与空间公平？

研究层次	可以进一步提出的问题	需要补充的资料
平台活动与时空热点	热点是否在工作日、周末、昼夜或不同年份反复出现？活动类型是否一致？	可比较的时间窗口、匿名用户、活动类型及独立资料。
交通可达性	从活动地点出发，步行多久可到站，乘车后能到达哪些就业、医疗或教育机会？	路网、站点入口、时刻表、换乘规则与目的地。
空间公平	哪些人群或地区的出行需求较高，却获得较少的服务机会？	居民分布、出行需求、年龄或收入等人群资料，以及明确的公平评价标准。

设施距离短，可以减少接近站点的空间阻隔；实际使用仍受入口、班次、票价和出行能力等条件影响（Cooke et al., 2022）。空间公平进一步追问，交通服务的机会与负担如何分配，哪些人群仍受到使用条件的限制（Verlinghieri & Schwanen, 2020）。这需要把服务供给与人群需求放在一起比较；仅凭微博点或站点数量，无法判断谁获得了更充分的服务。

已有研究利用地理标记社交媒体，按时段、活动类型和人口特征分析城市活动的时空差异，可作为本例分时分析的设计参考（Niu & Silva, 2023）。本例的匿名样本没有可核实的人口属性，不能直接从发帖内容推断用户的年龄或收入；若需比较人群，应另补充合适资料。

上海购物中心的一项研究整合抖音、小红书、微博、大众点评和百度指数，从关注、活动、体验三个层面评价商业空间活力（Zhang et al., 2025）。它为多源指标设计提供了一个本地范例：网络关注、位置记录和使用体验应分别测量，再考察它们是否一致。

活动与站点相邻也可能有多种解释：交通条件吸引活动，商业与就业场所同时吸引两者，或特定地点的发帖习惯造成记录集中。检验这些解释，可逐步加入用地、就业、场所类型与分时段数据；判断交通变化的作用，则需要设施变化前后的比较及合适的对照。

分析范围：2020–2025年上海微博4%随机样本，共215,968条记录，全部72个月均有记录；地图汇总整个时期。实际出行可达性与居民服务公平性留待补充资料后研究。

研究区域与数据

4.2 研究区域与数据获取

本例选择黄浦、徐汇、长宁、静安、杨浦五个中心城区，并保留闵行作为原案例的对照。各区面积和场所结构不同，适合比较高密度区的位置、形态及其交通设施邻近性。

数据	本例的作用	读取与解释口径
微博 4% 随机样本	识别记录分布与高密度区。	2020–2025 年，215,968 条；WGS84坐标，保留样本内匿名用户编号。
公交站、地铁站 POI	计算到最近设施的直线距离。	按文件声明读取坐标；设施搜索保留上海全市范围。
2025 年行政区划边界	确定空间筛选与展示范围。	使用更新的上海16区边界，专题2与旧练习边界文件分开保存。

距离计算以本例所用交通设施 POI 为参照。现有文件未注明设施采集日期与版本，要素是否对应独立车站、入口或其他设施位置仍需核实；POI 条数不能直接作为独立车站数，也不能代表各发帖年份的设施状态。

样本怎样得到？ 数据提供文件共 5,399,200 条记录。先对全部记录进行无放回简单随机抽样，固定种子为 20261016，抽出精确的4%；抽样前不按年份或地区筛选，也不去重。抽样后再按行政边界筛选研究范围。

课程样本只保留匿名用户编号、发布时间、经度、纬度、年份和月份。匿名编号仅用于识别本样本中的重复用户，不与旧版编号连接。样本与抽样说明位于 Example_data 数据仓库；下方代码使用本地 data/ 文件，便于离线复现。

4% 是记录的比例。随机抽样可减少课堂计算量，仍保留原始资料的平台覆盖与发帖选择特征。它不代表抽取了4%的居民或4%的独立访客。不同年份的记录量也有差异，跨年比较需进一步核对采集与平台使用条件（Petutschnig et al., 2021）。

运行前：核对文件与目录

工作目录设为第5章目录；代码开关位于图形下方，按编号运行，后续模块沿用前一步对象。

data/weibo_sample_4pct.csv: 2020–2025年微博4%随机样本
data/shanghai_boundary_2025.geojson: 上海16区行政边界
data/公交站.geojson: 全市公交设施位置
data/地铁站.geojson: 全市地铁设施位置
figures/weibo-multidistrict/amap_corrected_basemap.R: 高德瓦片坐标匹配辅助函数

运行顺序：依赖与参数 → 数据读取 → 底图与视野 → 地图版式 → 研究区总览。课程已使用 tidyverse，其中的读取、整理和绘图包无需重复加载；本节另加载 sf 与 patchwork，其他包通过 包名::函数名() 调用。首次配置时在 Console 补装缺少的包。

字段先规范：用 rename() 将来源表头转为英文 snake_case，后续使用 user_id、posted_at、longitude、latitude、year 和 month。原始CSV保留来源格式，中文用于图题、图例与表头。

计算与显示分别设置坐标：WGS84经纬度记录位置；KDE、面积与距离使用UTM 51N米制坐标，面积再换算为平方千米。地图采用Web Mercator显示，使真北朝上、经纬网横平竖直；显示投影不参与分析计算。

底图与坐标：高德瓦片怎样匹配？

底图默认使用高德地图，提供道路、水系与地名背景。将 basemap_provider <- "amap" 改为 "esri"，可切换到保留的Esri浅灰底图；两个来源的缓存分别保存。首次绘图需要联网，下载后的瓦片保存在 cache/map_tiles/。

高德瓦片采用GCJ-02坐标。辅助函数先对底图像素作近似坐标匹配，再将底图与矢量图层放入共同的显示投影。分析数据仍按WGS84读取、在UTM下计算。底图只提供地理背景，不用于测量距离或判定行政归属。

上海全市行政边界与六个研究区。红色半透明高亮研究范围；真北朝上，经纬网横平竖直，1–6编号与后续分区图一致。

底图来源：底图：高德地图（style 7）；仅对底图作坐标匹配重采样

1 · 依赖包、路径与共同参数

# 模块 1 · 依赖包、路径与共同参数

# tidyverse 已提供 readr、dplyr、tidyr 和 ggplot2，无需再逐一加载。
# 首次使用时在 Console 安装；课程环境已有 tidyverse 时，只补装其余包。
# install.packages("tidyverse")
# install.packages(c("sf", "patchwork", "MASS", "terra", "ggspatial",
#                    "ragg", "isoband", "rosm", "png", "digest"))
library(tidyverse)  # 数据读取、整理与 ggplot2 绘图
library(sf)         # 空间对象、投影转换与距离计算
library(patchwork)  # 多面板排版

# 其余包使用“包名::函数名()”调用。

# 工作目录设为第5章。
chapter_dir <- "."
data_dir <- file.path(chapter_dir, "data")
image_dir <- file.path(chapter_dir, "images")
result_dir <- file.path(chapter_dir, "figures", "weibo-multidistrict", "results")
tile_cache_dir <- file.path(chapter_dir, "cache", "map_tiles")
dir.create(result_dir, recursive = TRUE, showWarnings = FALSE)
dir.create(tile_cache_dir, recursive = TRUE, showWarnings = FALSE)

# 区名顺序在地图、表格和距离图中保持一致。保留闵行作为原案例参照。
district_order <- c("黄浦区", "徐汇区", "长宁区", "静安区", "杨浦区", "闵行区")
# 分析参数：米制坐标用于密度、面积和距离计算。
target_crs <- 32651       # UTM 51N，坐标单位为米；WGS84经纬度的单位是度。
display_crs <- 3857L     # 仅用于地图显示：真北朝上，经线竖直、纬线水平。
cell_m <- 75             # 评价网格边长；自动带宽较小的城区也需保留足够绘图细节。
bandwidth_multiplier <- 1 # 自动带宽倍率；敏感性试验可分别设为0.5、1或2。
aoi_quantile <- 0.95     # 每个区分别按区内密度网格的第95百分位数提取AOI。

# 显示参数：字体、配色与底图来源。
font_family <- "Microsoft YaHei"
ink <- "#253e49"
point_colour <- "#245d73"
study_colour <- "#C65B61"
boundary_colour <- "#496672"
aoi_colour <- "#00FF00"
# 密度图采用 Inferno 与透明度叠加，AOI 使用绿色虚线；这些设置只影响显示。
# 修改这一个选项即可切换底图，两个来源的缓存分别保留。
basemap_provider <- "amap"       # 可选 "amap" 或 "esri"

2 · 读取数据与规范字段

# 模块 2 · 读取数据与规范字段

# 1) 读取微博；匿名用户编号按字符保存，字段统一为snake_case。
source_columns <- c(
  user_id = "匿名用户ID",
  posted_at = "发布时间",
  longitude = "经度",
  latitude = "纬度",
  year = "年份",
  month = "月份"
)
weibo_raw <- read_csv(
  file.path(data_dir, "weibo_sample_4pct.csv"),
  col_types = "ccddii",
  show_col_types = FALSE
) %>%
  rename(all_of(source_columns))

# 2) 行政区界：投影转换、几何修复与研究区排序。
shanghai_districts <- st_read(
  file.path(data_dir, "shanghai_boundary_2025.geojson"),
  quiet = TRUE
) %>%
  st_transform(target_crs) %>%
  st_make_valid()

study_districts <- shanghai_districts %>%
  filter(district_name %in% district_order) %>%
  mutate(district_name = factor(district_name, levels = district_order)) %>%
  arrange(district_name)

# 3) 交通 POI：保留全市范围，供后续最近设施匹配使用。
# 保留文件声明的CRS，再进行转换；不靠st_set_crs()掩盖来源坐标不一致的问题。
poi_bus <- st_read(file.path(data_dir, "公交站.geojson"), quiet = TRUE)
poi_subway <- st_read(file.path(data_dir, "地铁站.geojson"), quiet = TRUE)
bus_proj <- st_transform(poi_bus, target_crs)
subway_proj <- st_transform(poi_subway, target_crs)

3 · 底图、导出与分区视野（辅助函数，可直接复用）

# 模块 3 · 底图、导出与分区视野

# 1) 输出设备与显示对象。
# PNG使用ragg，PDF使用Cairo。
export_board <- function(plot, name, width = 11.2, height = 16.2) {
  ggsave(
    file.path(image_dir, paste0(name, ".png")),
    plot,
    device = ragg::agg_png,
    width = width,
    height = height,
    units = "in",
    dpi = 300,
    bg = "white"
  )
  ggsave(
    file.path(result_dir, paste0(name, ".pdf")),
    plot,
    device = grDevices::cairo_pdf,
    width = width,
    height = height,
    units = "in",
    family = font_family,
    bg = "white"
  )
}

show_course_plot <- function(plot) {
  if (interactive() || isTRUE(getOption("knitr.in.progress"))) print(plot)
}

# 分析继续使用UTM 51N；显示几何另外转换为Web Mercator，不修改分析对象。
# Web Mercator中的真北朝向正上方；面积和距离仍由UTM分析几何计算。
shanghai_display <- st_transform(shanghai_districts, display_crs)
study_display <- st_transform(study_districts, display_crs) %>%
  mutate(district_id = factor(seq_len(n()), levels = seq_along(district_order)))

# 2) 底图：读取匹配函数，配置来源与缓存。
# 高德瓦片先匹配GCJ-02与WGS84；仅调整底图，分析几何保持不变。
# 首次绘制需联网，后续复用cache/map_tiles中的瓦片。
source(
  file.path(chapter_dir, "figures", "weibo-multidistrict", "amap_corrected_basemap.R"),
  encoding = "UTF-8"
)
esri_tile_url <- paste0(
  "https://services.arcgisonline.com/ArcGIS/rest/services/Canvas/",
  "World_Light_Gray_Base/MapServer/tile/${z}/${y}/${x}"
)
esri_credit <- paste0(
  "Tiles: Esri, HERE, Garmin, (c) OpenStreetMap contributors, ",
  "and the GIS user community"
)
amap_credit <- "底图：高德地图（style 7）；仅对底图作坐标匹配重采样"
basemap_credit <- if (basemap_provider == "amap") amap_credit else esri_credit
rosm::register_tile_source(
  weibo_light_gray = rosm::as.tile_source(
    esri_tile_url,
    extension = "jpg",
    max_zoom = 13,
    attribution = esri_credit
  )
)
basemap_layer <- function(window, zoom, alpha = 0.4, output_width = 1600L) {
  if (basemap_provider == "amap") {
    background <- amap_corrected_raster(
      window = window,
      zoom = zoom,
      cache_dir = tile_cache_dir,
      target_crs = display_crs,
      output_width = output_width
    )
    return(amap_corrected_layer(background, alpha = alpha))
  }
  ggspatial::annotation_map_tile(
    type = "weibo_light_gray",
    zoom = min(zoom, 13L),
    alpha = alpha,
    cachedir = tile_cache_dir,
    quiet = TRUE
  )
}

# 瓦片层级随显示视野选择，避免小区沿用低层级时把地名放得过大。
# 这是底图细节参数，与KDE的75米网格和自动带宽无关。
choose_tile_zoom <- function(window, target_pixels = 1600) {
  span <- max(
    window["xmax"] - window["xmin"],
    window["ymax"] - window["ymin"]
  )
  zoom <- round(log2(40075016.69 * target_pixels / (256 * span)))
  as.integer(max(10, min(16, zoom)))
}

# 3) 视野与经纬网：同一区在不同地图中保持共同范围。
# 用独立sf线几何画经纬网，避免网格被研究区填色或KDE面盖住。
map_graticule <- function(window) {
  extent <- st_as_sfc(st_bbox(
    setNames(
      unname(window[c("xmin", "ymin", "xmax", "ymax")]),
      c("xmin", "ymin", "xmax", "ymax")
    ),
    crs = st_crs(display_crs)
  ))
  st_graticule(extent, crs = st_crs(display_crs), datum = st_crs(4326))
}

# 每区独立缩放以便看清内部结构；同一区在点图/KDE/AOI中使用相同视野。
map_windows <- lapply(seq_len(nrow(study_display)), function(i) {
  box <- st_bbox(study_display[i, ])
  centre_x <- mean(box[c("xmin", "xmax")])
  centre_y <- mean(box[c("ymin", "ymax")])
  half_side <- 0.60 * max(
    box["xmax"] - box["xmin"],
    box["ymax"] - box["ymin"]
  )
  c(
    xmin = unname(centre_x - half_side),
    xmax = unname(centre_x + half_side),
    ymin = unname(centre_y - half_side),
    ymax = unname(centre_y + half_side)
  )
})
names(map_windows) <- district_order

4 · 地图主题、经纬网与制图标记（辅助函数，可直接复用）

# 模块 4 · 地图主题、经纬网与制图标记 ----------------------
# 沿用上一步的底图与视野；map_base() 和 map_references() 供后续分区图共用。

# 1) 共同主题：图例置右，地图边框与文字采用统一样式。
map_theme <- theme_void(base_family = font_family, base_size = 12) +
  theme(
    plot.title = element_text(
      colour = ink,
      face = "bold",
      size = 14,
      margin = margin(b = 4)
    ),
    plot.subtitle = element_text(
      colour = "#526574",
      size = 11,
      margin = margin(b = 5)
    ),
    panel.grid = element_blank(),
    panel.border = element_rect(fill = NA, colour = "black", linewidth = 1.0),
    plot.margin = margin(7, 7, 7, 7),
    plot.background = element_rect(fill = "white", colour = NA),
    panel.background = element_rect(fill = "white", colour = NA),
    legend.position = "right",
    legend.title = element_text(size = 11, face = "bold", colour = ink),
    legend.text = element_text(size = 10, colour = ink),
    legend.background = element_rect(fill = "white", colour = NA),
    axis.text = element_text(colour = ink, size = 8),
    axis.ticks = element_line(colour = "#71818c", linewidth = 0.25),
    axis.ticks.length = grid::unit(1.4, "mm"),
    plot.caption = element_blank()
  )

# 2) 图层底座：先铺底图，再叠加上海各区边界。
map_base <- function(region_name, index, subtitle = NULL, tile_alpha = 0.20) {
  ggplot() +
    basemap_layer(
      window = map_windows[[region_name]],
      zoom = choose_tile_zoom(map_windows[[region_name]]),
      alpha = tile_alpha
    ) +
    geom_sf(
      data = shanghai_display,
      fill = NA,
      colour = "#75828a",
      linewidth = 0.35
    ) +
    labs(title = paste(index, region_name), subtitle = subtitle) +
    map_theme
}

# 3) 定位标记与坐标：经纬网、比例尺、真北针和区内裁剪。
map_references <- function(region_name) {
  window <- map_windows[[region_name]]
  list(
    geom_sf(
      data = map_graticule(window),
      colour = "#6f7e89", linewidth = 0.24,
      linetype = "dashed", alpha = 0.48,
      inherit.aes = FALSE
    ),
    ggspatial::annotation_scale(
      location = "bl",
      width_hint = 0.27,
      pad_x = grid::unit(0.12, "cm"),
      pad_y = grid::unit(0.12, "cm"),
      height = grid::unit(0.08, "cm"),
      line_width = 0.35,
      text_cex = 0.65,
      text_family = font_family,
      text_col = ink
    ),
    # 显示投影已与真北对齐；仍明确让指北针读取真北方向。
    ggspatial::annotation_north_arrow(
      location = "br",
      which_north = "true",
      height = grid::unit(0.85, "cm"),
      width = grid::unit(0.85, "cm"),
      pad_x = grid::unit(0.16, "cm"),
      pad_y = grid::unit(0.18, "cm"),
      style = ggspatial::north_arrow_fancy_orienteering(
        line_width = 0.6,
        line_col = ink,
        fill = c("white", ink),
        text_col = ink,
        text_family = font_family,
        text_size = 7
      )
    ),
    # 坐标设置放在所有空间图层之后，防止后续geom_sf的默认坐标覆盖视野。
    coord_sf(
      crs = display_crs,
      default_crs = st_crs(display_crs),
      datum = st_crs(4326),
      label_axes = list(bottom = "E", left = "N"),
      xlim = window[c("xmin", "xmax")],
      ylim = window[c("ymin", "ymax")],
      expand = FALSE,
      # 分区视野外仍有全市区界；开启裁剪，避免线条越过框线与图例。
      clip = "on"
    )
  )
}

5 · 绘制研究区总览

# 模块 5 · 绘制研究区总览 ----------------------
# 使用显示用区界与底图；输出 study_area_plot，并保存 PNG 和 PDF。

# 总览保留上海16区的完整市域范围；所有岛屿也保留在市界内。
# 区域高亮采用透明填色，便于同时辨认底图道路、水系和区界。
# 1) 准备全市视野与编号位置。
overview_window <- st_bbox(shanghai_display)
overview_display_window <- setNames(
  unname(overview_window[c("xmin", "xmax", "ymin", "ymax")]) +
    c(-4500, 4500, -4500, 4500),
  c("xmin", "xmax", "ymin", "ymax")
)
city_boundary <- st_union(shanghai_display)
overview_labels <- study_display
st_geometry(overview_labels) <- st_point_on_surface(st_geometry(study_display))
# 2) 叠加底图、区界与研究区；编号对应右侧图例。
study_area_plot <- ggplot() +
  basemap_layer(
    window = overview_display_window,
    zoom = choose_tile_zoom(overview_display_window, target_pixels = 1800),
    alpha = 0.16,
    output_width = 1800L
  ) +
  geom_sf(
    data = shanghai_display,
    fill = NA,
    colour = "#566771",
    linewidth = 0.48
  ) +
  geom_sf(
    data = city_boundary,
    aes(colour = "city_boundary"),
    fill = NA,
    linewidth = 0.9
  ) +
  geom_sf(
    data = study_display,
    aes(fill = "study_area"),
    colour = "#A43742",
    linewidth = 0.9,
    alpha = 0.40
  ) +
  scale_fill_manual(
    values = c(study_area = study_colour),
    labels = c(study_area = "研究区"),
    name = NULL,
    guide = guide_legend(order = 1, override.aes = list(alpha = 0.5))
  ) +
  scale_colour_manual(
    values = c(city_boundary = "#3c5260"),
    labels = c(city_boundary = "上海市界"),
    name = NULL,
    guide = guide_legend(order = 2)
  ) +
  geom_sf(
    data = map_graticule(overview_display_window),
    colour = "#6f7e89", linewidth = 0.26,
    linetype = "dashed", alpha = 0.50,
    inherit.aes = FALSE
  ) +
  # 使用编号对照右侧图例，避免中心城区的长区名相互遮挡。
  geom_sf(
    data = overview_labels,
    shape = 21, size = 6.8, fill = "white", colour = "#A43742",
    stroke = 0.45, show.legend = FALSE
  ) +
  geom_sf(
    data = overview_labels,
    aes(shape = district_id),
    colour = ink,
    size = 4.8
  ) +
  scale_shape_manual(
    values = setNames(
      as.character(seq_along(district_order)), seq_along(district_order)
    ),
    labels = paste0(seq_along(district_order), "：", district_order),
    name = "研究区编号",
    guide = guide_legend(
      order = 3,
      override.aes = list(shape = NA),
      theme = theme(legend.key.width = grid::unit(0, "pt"))
    )
  ) +
  coord_sf(
    crs = display_crs,
    default_crs = st_crs(display_crs),
    datum = st_crs(4326),
    label_axes = list(bottom = "E", left = "N"),
    xlim = overview_display_window[c("xmin", "xmax")],
    ylim = overview_display_window[c("ymin", "ymax")],
    expand = FALSE,
    clip = "off"
  ) +
  ggspatial::annotation_north_arrow(
    location = "br",
    which_north = "true",
    height = grid::unit(1.0, "cm"),
    width = grid::unit(1.0, "cm"),
    pad_x = grid::unit(0.3, "cm"),
    pad_y = grid::unit(0.3, "cm"),
    style = ggspatial::north_arrow_fancy_orienteering(
      line_width = 0.7,
      line_col = ink,
      fill = c("white", ink),
      text_col = ink,
      text_family = font_family,
      text_size = 10
    )
  ) +
  ggspatial::annotation_scale(
    location = "bl",
    width_hint = 0.2,
    text_family = font_family,
    text_cex = 0.9,
    pad_x = grid::unit(0.4, "cm"),
    pad_y = grid::unit(0.35, "cm")
  ) +
  map_theme +
  theme(
    plot.margin = margin(14, 14, 14, 14),
    axis.text = element_text(colour = ink, size = 10),
    axis.ticks = element_line(colour = ink, linewidth = 0.35),
    axis.ticks.length = grid::unit(1.6, "mm"),
    legend.box = "vertical",
    legend.box.just = "left",
    legend.title = element_text(colour = ink, size = 13, face = "bold"),
    legend.text = element_text(colour = ink, size = 12),
    legend.key.height = grid::unit(0.60, "cm"),
    legend.key.width = grid::unit(0.60, "cm"),
    legend.spacing.y = grid::unit(0.18, "cm")
  )
# 3) 导出文件，并在交互运行或 Knit 时显示图形。
export_board(study_area_plot, "weibo-study-area-2020-2025", width = 11.5, height = 9.2)
show_course_plot(study_area_plot)

这一步应得到什么？

study_districts 保存六个研究区，weibo_raw 保存完整4%样本，bus_proj 与 subway_proj 保存转换后的全市交通设施点。总览中的红色仅表示研究区域，不编码微博数量。市界完整展示上海16区；图内编号与右侧图例对应，后续地图和云雨图沿用同一顺序：1黄浦、2徐汇、3长宁、4静安、5杨浦、6闵行。

数据预处理

4.3 数据预处理与初步展示

本节把微博经纬度转成空间点，先判断是否落在上海市界内，再分配到六个研究区。坐标合法与属于研究范围是两项不同的检查；界外记录不自动视为错误记录。

检查时间与坐标：保留有效经纬度与2020–2025年记录。
创建点几何并转换坐标：st_as_sf() 按WGS84生成点对象，st_transform() 转到米制投影。
按区截取：用点与区界的空间关系归属记录，边界点按固定规则处理，避免重复计数。每区后续只使用本区内的微博记录。
画出所有研究记录：六幅点图按1–6编号排列，保留点的真实位置；副标题列出记录数和单位面积的累计样本记录数。

2020–2025年上海微博4%随机样本的区内原始点分布。编号与研究区总览一致；副标题密度为累计样本记录数除以区面积。各区独立缩放，面板比例尺不同。

6 · 检查记录并按行政区筛选

# 模块 6 · 按时间、坐标和行政区筛选

# 1) 检查时间与坐标，并创建米制空间点。
# 先保留完整4%样本，再在分析阶段按明确规则筛选；不重新抽样或按区补齐数量。
# 源数据未提供时区；Asia/Shanghai为教学解析约定，开展分时研究前须核实。
weibo_clean <- weibo_raw %>%
  mutate(
    posted_at = as.POSIXct(
      posted_at,
      format = "%Y-%m-%d %H:%M:%S",
      tz = "Asia/Shanghai"
    )
  ) %>%
  filter(
    is.finite(longitude),
    is.finite(latitude),
    between(longitude, -180, 180),
    between(latitude, -90, 90),
    !is.na(posted_at),
    year %in% 2020:2025
  )
# 先声明原始经纬度的WGS84坐标系，再转换为米制投影；声明坐标系不会移动点。
weibo_points <- st_as_sf(
  weibo_clean,
  coords = c("longitude", "latitude"),
  crs = 4326,
  remove = FALSE
) %>%
  st_transform(target_crs)

# 2) 判断市内范围，再将记录唯一分配到研究区。
# 坐标合法不等于落在上海。行政边界判断另行进行，并记录范围外数量。
city_union <- st_union(shanghai_districts)
inside_city <- lengths(st_intersects(weibo_points, city_union)) > 0L
weibo_shanghai <- weibo_points[inside_city, ]
# st_intersects包含边界点；接触两区的记录按固定区序唯一归属，避免重复计数。
district_hits <- st_intersects(weibo_shanghai, study_districts)
district_index <- vapply(
  district_hits,
  function(hit) {
    if (length(hit)) hit[1] else NA_integer_
  },
  integer(1)
)

weibo_study <- weibo_shanghai[!is.na(district_index), ] %>%
  mutate(
    district_name = factor(
      district_order[district_index[!is.na(district_index)]],
      levels = district_order
    )
  )

# 3) 汇总各区面积、记录数与用户数。
district_summary <- study_districts %>%
  mutate(area_km2 = as.numeric(st_area(.)) / 1e6) %>%
  st_drop_geometry() %>%
  dplyr::select(district_name, district_code, area_km2) %>%
  left_join(
    st_drop_geometry(weibo_study) %>%
      group_by(district_name) %>%
      summarise(
        records = n(),
        users = n_distinct(user_id),
        .groups = "drop"
      ),
    by = "district_name"
  )

write_csv(district_summary, file.path(result_dir, "district_summary.csv"))

# 地图副标题使用累计样本记录数/区面积，不作总体推算。
map_display_summary <- district_summary %>%
  mutate(
    district_id = match(as.character(district_name), district_order),
    sample_record_density_km2 = records / area_km2
  ) %>%
  dplyr::select(
    district_id, district_name, district_code, area_km2,
    records, sample_record_density_km2
  )

weibo_study %>%
  st_drop_geometry() %>%
  count(district_name, year, month) %>%
  write_csv(file.path(result_dir, "district_month_counts.csv"))

7 · 绘制六区微博点分布

# 模块 7 · 绘制六区微博点分布 ----------------------
# 使用全部区内研究记录；各区沿用固定视野和编号。

# 逐区绘制全部研究点，再按统一顺序拼图；点的透明度只影响显示。
point_panels <- lapply(seq_along(district_order), function(i) {
  region_name <- district_order[i]
  points_i <- weibo_study %>% filter(district_name == region_name)
  map_base(
    region_name,
    i,
    paste0(
      format(nrow(points_i), big.mark = ","), " 条；样本密度 ",
      formatC(
        map_display_summary$sample_record_density_km2[i],
        format = "f", digits = 1, big.mark = ","
      ),
      " 条/km²"
    )
  ) +
    # aes()内的颜色映射生成图例；放大图例点不改变地图上的点大小。
    geom_sf(
      data = st_transform(points_i, display_crs),
      aes(colour = "微博记录（4%样本）"),
      alpha = 0.30,
      size = 0.16
    ) +
    scale_colour_manual(
      values = c("微博记录（4%样本）" = point_colour),
      name = NULL,
      guide = guide_legend(override.aes = list(size = 2.2, alpha = 1))
    ) +
    geom_sf(
      data = study_display[i, ],
      fill = NA,
      colour = "#253e49",
      linewidth = 0.90
    ) +
    map_references(region_name)
})
point_board <- (
  wrap_plots(point_panels, ncol = 2, guides = "collect")
) &
  theme(legend.position = "right", legend.margin = margin(0, 0, 0, 3))
export_board(point_board, "weibo-points-six-districts-2020-2025")
show_course_plot(point_board)

读图与检查

本次上海边界内有 215,403 条记录，其中六区共 112,569 条。各区数量与面积如下；用户数指样本内匿名代码数。

研究区	面积 km²	记录数	样本记录密度条/km²	样本内用户数
黄浦区	20.68	29900	1445.8	24902
徐汇区	53.48	22989	429.9	18587
长宁区	38.19	13057	341.9	10530
静安区	36.83	19688	534.6	16132
杨浦区	59.77	9160	153.3	7168
闵行区	376.22	17775	47.2	13296

副标题中的密度怎样理解？ 这是2020–2025年上海微博4%随机样本的累计记录数÷区面积，单位为条/km²，描述区内平均数量强度。它不代表全部微博的记录密度，本例不据此推算总体记录数。下一节KDE进一步观察区内不同位置的局部变化。

六区地图按两列三行、1–6编号排列，可点击图片放大查看。所有地图真北朝上，经纬网横平竖直。各面板比例尺不同，同一区在点图、KDE和AOI图中保持相同视野。点重叠会掩盖记录数量，局部色深也受透明度影响，不能只凭点图的颜色深浅判断密度。

核密度估计

4.4 核密度估计：从点到连续表面

核密度估计（Kernel Density Estimation, KDE）在每个记录位置放置一个平滑核，叠加各点的贡献，用于描述记录在空间上的相对集中程度。

本例按“区内点 → 自动带宽 KDE → AOI”逐区分析：

输入记录：每次只使用已经截取的本区微博点，不纳入邻区记录，也不设置区外缓冲区。
自动带宽：采用 MASS::kde2d() 的默认经验规则，由各区点坐标分别估计分析投影X、Y方向的带宽，大致对应东西、南北方向。代码将估计值显式传入 h，bandwidth_multiplier = 1 为基准；实际平滑尺度因区而异。
共同网格与单位：统一使用75米网格；将概率密度乘以本区输入记录数，再换算为样本记录/km²。像元值表示平滑强度，不能当作该像元内的精确记录数。

地图采用Inferno色带，由低到高叠加等值面，密度透明度为0.1–0.8。色带与透明度只影响显示，AOI仍按原密度计算。

跨区读图：各区实际带宽和色标不同，相同颜色不能直接比较跨区强度。先看区内位置与形态，再结合各自的色条数值判断密度。

带宽、网格与图面分层，分别控制什么？

带宽：平滑尺度

由各区数据自动估计。较大时相邻高值区可能合并，较小时保留更多局部起伏；高斯核带宽不能直接当作步行服务半径。

网格：评价位置

本例网格边长为75米，决定在哪些位置计算密度及表面细度。增加网格数量，不能代替选择合适的带宽。

分层：显示级数

bins = 20 控制绘制多少级密度等值面。它只用于绘图，不参与自动带宽选择，也不决定AOI阈值。

原理与延伸阅读：密度值、核函数与空间平滑

对等权记录，二维 KDE 可写为：

\[ \hat f(x,y)=\frac{1}{n h_x h_y}\sum_{i=1}^{n} K\!\left(\frac{x-x_i}{h_x},\frac{y-y_i}{h_y}\right). \]

其中 $n$ 是记录数，$h_x$、$h_y$ 是分析投影X、Y方向的平滑尺度。高斯核的贡献随距离逐渐衰减。

MASS::kde2d() 返回概率密度估计。坐标以米表示时，概率密度单位为每平方米；本例再乘以记录数和 $10^6$，得到样本记录/km²。不同软件的核函数和带宽约定可能不同，参数值不能直接照搬。

本例只用区内点估计密度，再按区界裁剪，未做边界质量补偿。高斯核的一部分会延伸到区外，因此区内栅格的积分不必等于输入记录数。区界附近的估计也会受到未纳入邻区记录的影响。

ArcGIS Pro：How Kernel Density works：核函数叠加的图解；其具体实现与 R 函数需分别核对。
eks：Kernel density estimates for tidy and geospatial data：带空间参考的 KDE 操作示例。

$六区KDE：仅使用各区内记录，统一75米网格。副标题的KDE bandwidth列出分析投影X、Y方向的实际带宽，单位为米；各图区内自动带宽和色标不同。$

六区KDE：仅使用各区内记录，统一75米网格。副标题的KDE bandwidth列出分析投影X、Y方向的实际带宽，单位为米；各图区内自动带宽和色标不同。

8 · 逐区估计带宽并计算 KDE

# 模块 8 · 逐区估计带宽并计算 KDE

# 每区只使用上一步裁取的本区记录，先自动估计带宽，再计算密度表面。
# 本例不增加区外缓冲记录；边界影响留在回顾与反思中讨论。
# 1) 定义单区计算：估计带宽 → 建评价网格 → 计算密度 → 按区界掩膜。
make_density_surface <- function(boundary, district_points) {
  xy <- st_coordinates(district_points)

  # MASS分别估计X、Y方向带宽；kde2d内部将传入的h除以4作为高斯核尺度。
  auto_h <- c(
    MASS::bandwidth.nrd(xy[, 1]),
    MASS::bandwidth.nrd(xy[, 2])
  )
  bandwidth_h <- auto_h * bandwidth_multiplier
  effective_h_m <- bandwidth_h / 4

  box <- st_bbox(boundary)
  x_edges <- seq(
    floor(box["xmin"] / cell_m) * cell_m,
    ceiling(box["xmax"] / cell_m) * cell_m,
    by = cell_m
  )
  y_edges <- seq(
    floor(box["ymin"] / cell_m) * cell_m,
    ceiling(box["ymax"] / cell_m) * cell_m,
    by = cell_m
  )
  x_centres <- head(x_edges, -1) + cell_m / 2
  y_centres <- head(y_edges, -1) + cell_m / 2

  # 将自动带宽乘以指定倍率后显式传入h；n只控制评价网格数量。
  # 固定网格不等于固定带宽，各区的平滑尺度仍由本区记录估计。
  density_fit <- MASS::kde2d(
    xy[, 1],
    xy[, 2],
    h = bandwidth_h,
    n = c(length(x_centres), length(y_centres)),
    lims = c(range(x_centres), range(y_centres))
  )
  grid_data <- expand.grid(x = density_fit$x, y = density_fit$y)
  # kde2d返回概率密度。乘本区实际输入点数和1e6，转为“样本记录/km²”；
  # 这里不乘25推算总体，更不能解释为居民人数密度。
  grid_data$intensity <- as.vector(density_fit$z) * nrow(xy) * 1e6
  surface <- terra::rast(grid_data, type = "xyz", crs = "EPSG:32651")
  # xyz中的x/y是像元中心，避免把评价点误当成像元边界的半格偏移。
  # touches=FALSE按像元中心判断是否在区内，区外设置为NA而非0。
  surface <- terra::mask(surface, terra::vect(boundary), touches = FALSE)
  # 绘图保留完整矩形矩阵供色带插值，再将色带裁至区界。
  # 此处不把掩膜外NA补为0，也不增加区外输入点。
  density_field <- list(
    x = density_fit$x,
    y = density_fit$y,
    z = density_fit$z * nrow(xy) * 1e6
  )
  list(
    raster = surface,
    source_n = nrow(xy),
    hx_m = effective_h_m[1],
    hy_m = effective_h_m[2],
    field = density_field,
    grid = terra::as.data.frame(surface, xy = TRUE, na.rm = TRUE)
  )
}

# 2) 对六区分别调用函数；各区带宽取决于本区数据。
density_surfaces <- lapply(seq_along(district_order), function(i) {
  district_points <- weibo_study %>% filter(district_name == district_order[i])
  surface <- make_density_surface(study_districts[i, ], district_points)
  terra::writeRaster(
    surface$raster,
    file.path(result_dir, paste0("kde_", study_districts$district_code[i], ".tif")),
    overwrite = TRUE,
    datatype = "FLT4S",
    gdal = "COMPRESS=LZW"
  )
  surface
})
names(density_surfaces) <- district_order
# 3) 保存实际高斯核带宽及密度栅格，供读图和结果核对。
bandwidth_summary <- tibble(
  district_name = district_order,
  district_code = study_districts$district_code,
  records = vapply(density_surfaces, function(x) x$source_n, integer(1)),
  hx_m = vapply(density_surfaces, function(x) x$hx_m, numeric(1)),
  hy_m = vapply(density_surfaces, function(x) x$hy_m, numeric(1)),
  grid_m = cell_m
)
write_csv(bandwidth_summary, file.path(result_dir, "bandwidth_summary.csv"))
map_display_summary <- map_display_summary %>%
  left_join(
    bandwidth_summary %>% dplyr::select(district_name, hx_m, hy_m, grid_m),
    by = "district_name"
  )
district_maxima <- vapply(
  density_surfaces,
  function(x) max(x$grid$intensity),
  numeric(1)
)

9 · 绘制密度等值面

9.1 等值面与色标辅助函数，可直接复用

# 模块 9a · 定义密度等值面与色标辅助函数
# 将上一步的密度表面转换为显示图层；不在绘图时重新估计密度。

make_density_bands <- function(density_field, boundary, bins = 20L) {
  # 1) 整理矩阵和区界，保留原密度单位。
  x <- as.double(density_field$x)
  y <- as.double(density_field$y)
  z <- density_field$z
  storage.mode(z) <- "double"
  boundary_geometry <- sf::st_make_valid(sf::st_geometry(boundary))
  boundary_geometry <- sf::st_union(boundary_geometry)
  crs <- sf::st_crs(boundary)
  empty_result <- function() {
    sf::st_sf(
      level_low = numeric(),
      level_high = numeric(),
      level_mid = numeric(),
      geometry = sf::st_sfc(crs = crs)
    )
  }
  display_max <- density_field$display_max
  if (is.null(display_max)) display_max <- max(z)
  if (all(sf::st_is_empty(boundary_geometry)) || max(z) == 0 || display_max == 0) {
    return(empty_result())
  }

  # 2) 按区内最大值构造显示层级；不改变密度数值或AOI阈值。
  breaks <- seq(0, display_max, length.out = as.integer(bins) + 1L)
  level_low <- head(breaks, -1L)
  level_high <- tail(breaks, -1L)
  level_mid <- (level_low + level_high) / 2
  # 嵌套面由低到高绘制，累积透明度。
  level_high[] <- Inf
  # kde2d的z按[x,y]存储；isoband按[行y,列x]读取，必须转置。
  bands <- isoband::isobands(
    x = x,
    y = y,
    z = t(z),
    levels_low = level_low,
    levels_high = level_high
  )
  band_geometry <- sf::st_sfc(isoband::iso_to_sfg(bands), crs = crs)
  keep <- !sf::st_is_empty(band_geometry)
  if (!any(keep)) return(empty_result())
  band_data <- data.frame(
    level_low = level_low[keep],
    level_high = level_high[keep],
    level_mid = level_mid[keep]
  )
  band_geometry <- sf::st_make_valid(band_geometry[keep])

  # 3) 裁剪显示面，并移除空几何；不改分析栅格。
  # 先生成完整等值面，再严格裁至区界，避免NA掩膜造成边缘色带缺口。
  clipped <- sf::st_intersection(band_geometry, boundary_geometry)
  if (!length(clipped)) return(empty_result())
  source_index <- attr(clipped, "idx")[, 1L]
  result <- sf::st_sf(band_data[source_index, , drop = FALSE], geometry = clipped)
  result <- sf::st_make_valid(result)
  if (any(sf::st_geometry_type(result) == "GEOMETRYCOLLECTION")) {
    result <- sf::st_collection_extract(result, "POLYGON", warn = FALSE)
  }
  keep <- sf::st_geometry_type(result) %in% c("POLYGON", "MULTIPOLYGON") &
    !sf::st_is_empty(result)
  result <- result[keep, , drop = FALSE]
  if (!nrow(result)) return(empty_result())
  result <- result[order(result$level_low), , drop = FALSE]
  rownames(result) <- NULL
  result
}

# 定义色标：密度 alpha 随显示层级变化，底图 alpha 另行设置。
# 每区保留独立数值色标，便于查看本区结构。
# 相同颜色在不同区可对应不同强度，跨区比较须读取各自图例数值。
heat_scale <- function(maximum) {
  list(
    scale_fill_viridis_c(
      option = "inferno",
      limits = c(0, maximum),
      oob = scales::squish,
      name = "样本记录\n/ km²",
      breaks = scales::breaks_extended(n = 3),
      guide = guide_colourbar(
        barwidth = grid::unit(0.22, "cm"),
        barheight = grid::unit(2.5, "cm"),
        title.position = "top"
      )
    ),
    scale_alpha_continuous(
      range = c(0.1, 0.8),
      limits = c(0, maximum),
      oob = scales::squish,
      guide = "none"
    ),
    theme(
      legend.position = "right",
      legend.title = element_text(size = 9.0, face = "bold"),
      legend.text = element_text(size = 8.5),
      legend.margin = margin(0, 0, 0, 1)
    )
  )
}

9.2 生成六区显示图层并绘图

# 模块 9b · 绘制六区密度等值面

density_bands <- lapply(seq_along(density_surfaces), function(i) {
  density_field <- density_surfaces[[i]]$field
  density_field$display_max <- district_maxima[i]
  make_density_bands(density_field, study_districts[i, ], bins = 20)
})
# 等值带仅在显示时重投影；原KDE数值、UTM评价点和AOI阈值均不重算。
density_bands_display <- lapply(density_bands, st_transform, crs = display_crs)
kde_panels <- lapply(seq_along(district_order), function(i) {
  region_name <- district_order[i]
  map_base(
    region_name,
    i,
    sprintf(
      "KDE bandwidth：X %.0f m / Y %.0f m",
      bandwidth_summary$hx_m[i],
      bandwidth_summary$hy_m[i]
    ),
    tile_alpha = 0.12
  ) +
    geom_sf(
      data = density_bands_display[[i]],
      aes(fill = level_mid, alpha = level_mid),
      colour = NA
    ) +
    heat_scale(district_maxima[i]) +
    geom_sf(
      data = study_display[i, ],
      fill = NA,
      colour = "black",
      linewidth = 0.90
    ) +
    map_references(region_name)
})
kde_board <- wrap_plots(kde_panels, ncol = 2, guides = "keep")
export_board(kde_board, "weibo-kde-six-districts-2020-2025")
show_course_plot(kde_board)

本次自动带宽：地图副标题与下表报告高斯核实际使用的X、Y方向带宽，单位为米。这两个方向属于UTM分析投影，与真北朝上的显示方向略有差别；默认经验规则不表示已经找到最优带宽。

研究区	记录数	X方向带宽 m	Y方向带宽 m	网格 m
黄浦区	29900	134.4	170.9	75
徐汇区	22989	159.9	352.0	75
长宁区	13057	466.2	168.8	75
静安区	19688	146.1	298.8	75
杨浦区	9160	282.7	339.9	75
闵行区	17775	635.3	1049.6	75

这张密度图怎么读？

先查看各区高值部分的位置和连片情况，再结合本区带宽及图例读取密度。较小带宽通常保留更多局部起伏，较大带宽可能连接邻近高值区；跨区形态差异也会受到平滑尺度的影响。地图汇总2020–2025年的记录，高值区的活动类型与持续性仍需进一步核对。

接下来的任务是为高值区域建立可重复的分析边界，以便比较区域内外的记录位置。该边界由数据与参数共同决定。

下方为模拟点平滑演示：画布坐标与带宽使用像素，颜色来自未归一化的核叠加值，用于观察平滑形态。它不使用本例的真实微博数据，滑块数值不能直接套入 R 的空间分析。

边界提取与可视化

4.5 高密度区提取与 AOI 边界

本节从刚才计算的同一份KDE表面提取高值分析区，称为 AOI（Area of Interest）。绘图时不再重新估计另一份密度。

限定评价范围：只使用像元中心落在相应区内的有效网格。
逐区取阈值：分别计算每个区内密度的第95百分位数，提取达到该阈值的高值网格。
转面并裁剪：合并相邻高值网格，与区界相交，得到可用于记录分组的多边形。
统计连通部分：共享边或顶点、可以彼此连通的AOI合计为1个；面积按UTM投影计算所有AOI的合计面积。
叠加核对：AOI边界叠在原KDE上，检查位置、范围与边界是否对应。

Top 5% 是区内排名，六区的绝对阈值不同。 边界网格裁剪、密度并列等情况会使AOI面积不恰好等于区面积的5%。AOI随数据与参数变化，不代表固定城市功能区，也不等同于统计显著热点。

同一KDE表面上提取的AOI。绿色虚线表示各区内密度排名前5%的高值分析区；副标题列出连通AOI个数及合计面积，共享边或顶点的区域计为同一个。各图色标不同。

10 · 提取 AOI 并统计连通区域

10.1 按阈值提取 AOI，划分记录内外组

# 模块 10a · 提取 AOI 并按内外归组

# 阈值在各区的有效网格内分别计算，不能把六区合在一起取95%分位。
# 各区的绝对阈值不同，所得AOI用于区内高值区/其余地区的分组。
# 1) 逐区阈值化、转为多边形，再严格裁至本区边界。
aoi_list <- lapply(seq_along(district_order), function(i) {
  surface <- density_surfaces[[i]]$raster
  threshold <- unname(quantile(density_surfaces[[i]]$grid$intensity, aoi_quantile))
  high_cells <- terra::ifel(surface >= threshold, 1, NA)
  high_polygons <- terra::as.polygons(
    high_cells,
    dissolve = TRUE,
    na.rm = TRUE
  ) %>%
    st_as_sf() %>%
    st_make_valid()
  # 网格面还需与区界相交；边界格被裁切后，面积不保证恰好占全区5%。
  clipped <- st_intersection(
    st_geometry(high_polygons),
    st_geometry(study_districts[i, ])
  )
  st_sf(
    district_name = district_order[i],
    threshold = threshold,
    geometry = st_union(clipped)
  )
})
aoi_regions <- do.call(rbind, aoi_list)
st_write(
  st_transform(aoi_regions, 4326),
  file.path(result_dir, "aoi_six_districts.geojson"),
  delete_dsn = TRUE,
  quiet = TRUE
)

# 2) 将原始记录分为 AOI 内外，并汇总面积与记录覆盖率。
# 区内样本按自身所在区的AOI分类，不用其他区的AOI给本区记录分组。
aoi_hits <- st_intersects(weibo_study, aoi_regions)
own_aoi_index <- match(as.character(weibo_study$district_name), aoi_regions$district_name)
aoi_inside <- mapply(
  function(hit, own_index) own_index %in% hit,
  aoi_hits,
  own_aoi_index,
  USE.NAMES = FALSE
)
weibo_study$aoi_group <- factor(
  ifelse(aoi_inside, "AOI 内", "AOI 外"),
  levels = c("AOI 内", "AOI 外")
)
aoi_summary <- district_summary %>%
  left_join(st_drop_geometry(aoi_regions), by = "district_name") %>%
  mutate(
    aoi_km2 = as.numeric(st_area(aoi_regions)) / 1e6,
    aoi_area_percent = 100 * aoi_km2 / area_km2
  ) %>%
  left_join(
    st_drop_geometry(weibo_study) %>%
      group_by(district_name) %>%
      summarise(
        aoi_records = sum(aoi_group == "AOI 内"),
        aoi_record_percent = 100 * mean(aoi_group == "AOI 内"),
        .groups = "drop"
      ),
    by = "district_name"
  )
write_csv(aoi_summary, file.path(result_dir, "aoi_summary.csv"))

10.2 统计连通部分与地图指标

# 模块 10b · 统计连通 AOI 并汇总地图指标

# 按共边或共点关系统计连通AOI；孔洞不单独计数。
count_aoi_components <- function(aoi_geometry) {
  polygon_geometry <- st_make_valid(st_geometry(aoi_geometry))
  if (any(st_geometry_type(polygon_geometry) == "GEOMETRYCOLLECTION")) {
    polygon_geometry <- st_collection_extract(
      polygon_geometry, "POLYGON", warn = FALSE
    )
  }
  polygon_parts <- polygon_geometry %>%
    st_union() %>%
    st_cast("POLYGON", warn = FALSE)
  polygon_parts <- polygon_parts[!st_is_empty(polygon_parts)]
  if (!length(polygon_parts)) return(0L)

  neighbours <- st_intersects(polygon_parts)
  visited <- rep(FALSE, length(polygon_parts))
  component_count <- 0L
  for (start in seq_along(polygon_parts)) {
    if (visited[start]) next
    component_count <- component_count + 1L
    pending <- start
    visited[start] <- TRUE
    while (length(pending)) {
      current <- pending[1L]
      pending <- pending[-1L]
      new_parts <- neighbours[[current]][!visited[neighbours[[current]]]]
      visited[new_parts] <- TRUE
      pending <- c(pending, new_parts)
    }
  }
  component_count
}

aoi_map_metrics <- aoi_summary %>%
  dplyr::select(district_name, aoi_km2) %>%
  mutate(
    aoi_connected_components = vapply(
      seq_len(nrow(aoi_regions)),
      function(i) count_aoi_components(aoi_regions[i, ]),
      integer(1)
    )
  )
map_display_summary <- map_display_summary %>%
  left_join(aoi_map_metrics, by = "district_name")
write_csv(map_display_summary, file.path(result_dir, "map_display_summary.csv"))

11 · 叠加 AOI 与密度表面

# 模块 11 · 叠加 AOI 与密度表面 ----------------------
# 使用已提取的 AOI；显示投影转换不改变面积与内外分组。

# AOI与KDE色带使用同一显示投影；分析对象aoi_regions继续保留UTM坐标。
aoi_display <- st_transform(aoi_regions, display_crs)
aoi_panels <- lapply(seq_along(district_order), function(i) {
  region_name <- district_order[i]
  map_base(
    region_name,
    i,
    paste0(
      "连通AOI ", map_display_summary$aoi_connected_components[i], " 个；总面积 ",
      format(round(aoi_summary$aoi_km2[i], 2), nsmall = 2),
      " km²"
    ),
    tile_alpha = 0.16
  ) +
    geom_sf(
      data = density_bands_display[[i]],
      aes(fill = level_mid, alpha = level_mid),
      colour = NA
    ) +
    heat_scale(district_maxima[i]) +
    geom_sf(
      data = study_display[i, ],
      fill = NA,
      colour = "black",
      linewidth = 0.90
    ) +
    geom_sf(
      data = aoi_display[i, ],
      aes(colour = "AOI 边界"),
      fill = NA,
      linewidth = 0.75,
      linetype = "dashed"
    ) +
    scale_colour_manual(
      values = c("AOI 边界" = aoi_colour),
      labels = c("AOI 边界" = "AOI边界"),
      name = NULL,
      guide = guide_legend(
        order = 2,
        override.aes = list(linewidth = 0.8, alpha = 1)
      )
    ) +
    map_references(region_name)
})
# 色条与AOI线型均放在各面板右侧；独立色标不合并，避免混淆各区强度。
aoi_board <- wrap_plots(aoi_panels, ncol = 2, guides = "keep")
export_board(aoi_board, "weibo-aoi-six-districts-2020-2025")
show_course_plot(aoi_board)

提取后检查什么？

AOI 是否与密度表面的高值部分对应？边界附近是否超出研究范围？
连片与分散部分各在哪里？改变阈值后，面积、数量和连通性如何变化？
每组包含多少记录？是否由少数用户、重复位置或短期事件主导？

下表列出本次AOI的连通个数、阈值、合计面积和记录覆盖比例，便于对照地图检查。连通个数不等于独立功能区或活动中心的数量，它会受到网格、带宽和阈值影响。记录覆盖比例指样本记录的位置落在AOI内的比例，不能解释为人口比例。

研究区	连通AOI数	密度阈值记录/km²	AOI合计 km²	面积占比 %	记录覆盖 %
黄浦区	4	5012.6	1.03	5.01	34.02
徐汇区	4	1775.9	2.66	4.98	37.89
长宁区	5	937.3	1.91	5.00	32.27
静安区	3	2240.3	1.84	5.00	39.58
杨浦区	3	486.0	2.99	5.00	36.24
闵行区	4	163.1	18.81	5.00	38.34

AOI内外分组建立后，下一步比较各记录位置到最近交通设施的距离。

设施邻近性对比

4.6 到公交站、地铁站的距离比较

本节分别计算六区微博记录位置到最近公交站、地铁站的直线距离，在各区内比较AOI内外的分布。公交与地铁分开作图，同一设施图使用共同距离轴。

最近设施的搜索范围

微博记录按研究区分组，最近设施在上海全市公交或地铁POI中寻找，以保留区外的近邻站点。点到点的直线距离用于衡量设施邻近性；实际可达性还受入口、班次与出行条件等影响（Cooke et al., 2022）。

代码阅读顺序

统一米制投影，使微博点与设施点可以计算距离。
关联 AOI 分组，依据记录位置与多边形的关系，标记 AOI 内外。
匹配最近设施，用 st_nearest_feature() 找到对应站点，再以 by_element = TRUE 逐对计算距离，避免创建完整距离矩阵。
查看云雨图，结合半密度轮廓、细箱线与点阵，观察中间分布和长尾；红色菱形标出原始距离的均值。
比较均值差，在每个区内估计AOI外减AOI内的均值差，按用户账号计算聚类稳健标准误，并对12项比较统一作Holm校正。

图下代码依次分为距离计算、统计检验、分布准备、点阵整理、绘图函数、出图、保存结果七个模块。完整函数保留在同一开关内，按内部编号阅读；无需在运行各模块前重复安装或加载包。

完整记录与展示子样本，分别用于什么？

箱体、胡须、均值与检验使用全部有效记录，不因远距离或箱线图离群标记而删除记录。横轴用 pseudo_log(sigma = 50) 压缩长尾：0附近近似线性，较远距离逐渐压缩，标签仍以米表示。

规则点阵固定随机种子，每组最多随机展示260条记录，另补真实最小值与最大值。每点最多表示5条展示记录，不足5条时颜色较淡；点位对齐分箱中心，不能读取精确距离。各组总量看记录数标签和统计表，不能按点的多少比较。

半密度轮廓由全部记录在压缩视图中估计，在0处反射以避免负距离；各组轮廓归一化到相同高度，用于观察分布形状。

到最近公交POI的云雨图。箱体和均值保留完整距离，横轴压缩长尾；星号对应用户聚类稳健均值检验，经12项统一Holm校正。

到最近地铁POI的云雨图。均值差为AOI外减AOI内，方括号为逐项95%置信区间。规则点阵来自固定子样本，每点最多5条；不足5条颜色更淡，点位对齐分箱中心，不代表精确距离。

12 · 计算最近设施距离与分组统计

# 模块 12 · 计算最近设施距离与分组统计

# 最近设施可在研究区之外，因此保留全上海公交/地铁POI参与匹配。
# st_nearest_feature给出最近设施的行号；by_element逐条计算对应点对的距离。
# 两类几何已统一为UTM米制投影，因此st_distance返回的距离单位为米。
weibo_study$subway_m <- as.numeric(
  st_distance(
    weibo_study,
    subway_proj[st_nearest_feature(weibo_study, subway_proj), ],
    by_element = TRUE
  )
)
weibo_study$bus_m <- as.numeric(
  st_distance(
    weibo_study,
    bus_proj[st_nearest_feature(weibo_study, bus_proj), ],
    by_element = TRUE
  )
)
distance_long <- st_drop_geometry(weibo_study) %>%
  dplyr::select(district_name, user_id, aoi_group, bus_m, subway_m) %>%
  pivot_longer(
    c(bus_m, subway_m),
    names_to = "facility",
    values_to = "distance_m"
  ) %>%
  mutate(
    facility = factor(
      facility,
      levels = c("bus_m", "subway_m"),
      labels = c("公交", "地铁")
    )
  )
distance_summary <- distance_long %>%
  group_by(district_name, facility, aoi_group) %>%
  summarise(
    n = n(),
    users = n_distinct(user_id),
    median_m = median(distance_m),
    q25_m = quantile(distance_m, 0.25),
    q75_m = quantile(distance_m, 0.75),
    p95_m = quantile(distance_m, 0.95),
    max_m = max(distance_m),
    within_500m_percent = 100 * mean(distance_m <= 500),
    .groups = "drop"
  )
write_csv(distance_summary, file.path(result_dir, "distance_summary.csv"))

13 · 用户聚类稳健检验与 Holm 校正

本块的检验函数可直接复用。阅读时重点核对“AOI 外 − AOI 内”的均值差方向、按匿名账号处理重复记录的标准误，以及12项比较的 Holm 校正；矩阵计算可留作进阶阅读。

# 模块 13 · 用户聚类稳健检验与 Holm 校正

# 目标量为 AOI 外减 AOI 内的记录加权均值差，不把重复记录当作独立访客。
# CR1 聚类稳健标准误允许同一账号内的记录相关；区间为逐项95%置信区间。
# AOI 在此视为固定；该探索性比较未传播 AOI 估计或跨用户空间相关的不确定性。
# 1) 定义单项检验：估计差值、用户聚类标准误与逐项置信区间。
cluster_mean_test <- function(data) {
  model_data <- data.frame(
    distance_m = data$distance_m,
    group_out = as.integer(data$aoi_group == "AOI 外"),
    user_id = data$user_id
  )
  fit <- lm(distance_m ~ group_out, data = model_data)
  design <- model.matrix(fit)
  cluster_scores <- rowsum(
    design * residuals(fit), model_data$user_id, reorder = FALSE
  )
  n_records <- nrow(model_data)
  n_users <- n_distinct(model_data$user_id)

  bread <- solve(crossprod(design))
  correction <- n_users / (n_users - 1) *
    (n_records - 1) / (n_records - ncol(design))
  robust_vcov <- correction *
    bread %*% crossprod(cluster_scores) %*% bread
  estimate_m <- unname(coef(fit)[[2]])
  se_m <- sqrt(robust_vcov[2, 2])
  critical_t <- qt(0.975, df = n_users - 1)

  tibble(
    n_records = n_records, n_user_clusters = n_users,
    mean_inside_m = mean(data$distance_m[data$aoi_group == "AOI 内"]),
    mean_outside_m = mean(data$distance_m[data$aoi_group == "AOI 外"]),
    estimate_out_minus_in_m = estimate_m,
    se_user_cluster_m = se_m,
    ci_low_m = estimate_m - critical_t * se_m,
    ci_high_m = estimate_m + critical_t * se_m,
    t_statistic = estimate_m / se_m,
    df = n_users - 1,
    p_raw = 2 * pt(-abs(estimate_m / se_m), df = n_users - 1)
  )
}

# 2) 完成“六区×两类设施”的12项比较，再统一作Holm校正。
distance_tests <- distance_long %>%
  group_by(district_name, facility) %>%
  group_modify(~ cluster_mean_test(.x)) %>%
  ungroup() %>%
  mutate(
    p_holm = p.adjust(p_raw, method = "holm"),
    significance = case_when(
      p_holm < 0.001 ~ "***",
      p_holm < 0.01 ~ "**",
      p_holm < 0.05 ~ "*",
      TRUE ~ "ns"
    )
  )
write_csv(distance_tests, file.path(result_dir, "distance_cluster_tests.csv"))

14 · 准备云雨图统计与半密度轮廓

# 模块 14 · 准备云雨图统计与半密度轮廓

# pseudo_log 在 0 附近平滑过渡；这里只变换显示，不删远点或改变检验目标。
# 箱体和 1.5 IQR 胡须在原始米制距离上计算；密度在压缩视图中估计。
# 1) 为图面建立压缩距离坐标，保留原始 distance_m。
distance_transform <- scales::transform_pseudo_log(sigma = 50, base = 10)
distance_view <- distance_long %>%
  mutate(
    group_id = if_else(aoi_group == "AOI 内", 2, 1),
    distance_view = distance_transform$transform(distance_m),
    record_id = row_number()
  )
distance_colours <- c("AOI 内" = "#5B91A6", "AOI 外" = "#C58B70")

# 2) 先按原始米数计算箱线和均值，再将坐标转换到压缩视图。
box_summary <- distance_view %>%
  group_by(district_name, facility, aoi_group, group_id) %>%
  group_modify(~ {
    box_values <- boxplot.stats(.x$distance_m, do.conf = FALSE)$stats
    tibble(
      n_records = nrow(.x), n_users = n_distinct(.x$user_id),
      whisker_low_m = box_values[1], q25_m = box_values[2],
      median_m = box_values[3], q75_m = box_values[4],
      whisker_high_m = box_values[5], mean_m = mean(.x$distance_m)
    )
  }) %>%
  ungroup() %>%
  mutate(
    across(
      c(whisker_low_m, q25_m, median_m, q75_m, whisker_high_m, mean_m),
      distance_transform$transform, .names = "{.col}_view"
    ),
    count_label = paste0(
      "记录 ", format(n_records, big.mark = ",", trim = TRUE),
      " · 用户 ", format(n_users, big.mark = ",", trim = TRUE)
    )
  )

write_csv(
  box_summary %>%
    dplyr::select(
      district_name, facility, aoi_group, n_records, n_users,
      whisker_low_m, q25_m, median_m, q75_m, whisker_high_m, mean_m
    ),
  file.path(result_dir, "distance_box_summary.csv")
)

# 3) 反射边界估计半密度轮廓；全部有效记录参与。
# 将视图中的距离在 0 处反射后估计密度，消除不合理的负距离尾部。
# 所有有效记录参与估计；轮廓正值尾端自然收尖，曲线不代表新增观测。
half_density <- function(values) {
  bandwidth <- stats::bw.nrd0(values)
  density_values <- stats::density(
    c(-values, values), bw = 0.8 * bandwidth,
    from = 0, to = max(values) + 3 * bandwidth, n = 512
  )
  tibble(
    distance_view = c(0, density_values$x, tail(density_values$x, 1)),
    density_height = c(0, density_values$y / max(density_values$y), 0)
  )
}
density_polygons <- distance_view %>%
  group_by(district_name, facility, aoi_group, group_id) %>%
  group_modify(~ half_density(.x$distance_view)) %>%
  ungroup() %>%
  mutate(density_y = group_id + 0.045 + 0.29 * density_height)

15 · 整理规则点阵

# 模块 15 · 整理规则点阵

# 点阵用固定子样本；统计量仍使用全部记录。
# 同一分箱与堆叠规则下，每点最多显示5条，末点按余数调整透明度。
dot_binwidth <- 0.028
dot_records_per_point <- 5L
dot_stack_step <- 0.024

# 固定种子使每次展示同一批子样本；每组最多260条随机记录，另补极值。
set.seed(20261016)
display_points <- distance_view %>%
  group_by(district_name, facility, aoi_group, group_id) %>%
  group_modify(~ {
    selected_rows <- unique(c(
      sample.int(nrow(.x), min(260L, nrow(.x))),
      which.min(.x$distance_m), which.max(.x$distance_m)
    ))
    .x[selected_rows, ]
  }) %>%
  ungroup() %>%
  mutate(
    dot_bin = floor(distance_view / dot_binwidth),
    dot_x = (dot_bin + 0.5) * dot_binwidth
  ) %>%
  arrange(district_name, facility, aoi_group, dot_bin, distance_view) %>%
  group_by(district_name, facility, aoi_group, group_id, dot_bin) %>%
  mutate(
    dot_rank = row_number(),
    dot_stack = ceiling(dot_rank / dot_records_per_point)
  ) %>%
  ungroup()

# 在每个分箱内部汇总显示点；分箱总记录数与抽样集合严格一致。
dot_cells <- display_points %>%
  group_by(district_name, facility, aoi_group, group_id, dot_bin, dot_stack) %>%
  summarise(dot_x = first(dot_x), dot_records = n(), .groups = "drop") %>%
  mutate(
    dot_y = group_id - 0.13 - dot_stack_step * (dot_stack - 1),
    dot_alpha = dot_records / dot_records_per_point
  )

16 · 定义云雨图图层与坐标（辅助函数，可直接复用）

# 模块 16 · 定义云雨图图层与坐标 ----------------------
# 组合前两步的图层与 distance_tests；定义函数，此处尚不出图。

format_holm_p <- function(values) {
  ifelse(
    values < 0.001, "Holm P < 0.001",
    paste0("Holm P = ", formatC(values, format = "f", digits = 3))
  )
}

distance_plot <- function(facility_name) {
  # 1) 选择设施，并为完整分布和检验标记预留空间。
  filter_facility <- function(data) filter(data, facility == facility_name)
  box_data <- filter_facility(box_summary)
  cloud_data <- filter_facility(density_polygons)
  dot_data <- filter_facility(dot_cells)
  test_data <- filter_facility(distance_tests)

  # 六区共用完整范围，括号留在最远轮廓之外；不按显著性选择视野。
  maximum_view <- max(cloud_data$distance_view)
  bracket_x <- maximum_view * 1.035
  test_data <- test_data %>%
    mutate(
      bracket_x = bracket_x,
      result_label = paste0(
        "均值差（外 − 内） ", sprintf("%+.1f", estimate_out_minus_in_m),
        " m  [", sprintf("%.1f", ci_low_m), ", ",
        sprintf("%.1f", ci_high_m), "]；", format_holm_p(p_holm)
      )
    )
  break_m <- c(0, 50, 100, 250, 500, 1000, 2000, 5000, 10000, 20000)
  break_m <- break_m[
    distance_transform$transform(break_m) <= maximum_view * 1.02
  ]

  # 2) 云与雨分别使用完整分布轮廓和规则展示点阵。
  ggplot() +
    geom_polygon(
      data = cloud_data,
      aes(distance_view, density_y, group = aoi_group, fill = aoi_group),
      alpha = 0.56, colour = "#51646E", linewidth = 0.28
    ) +
    geom_point(
      data = dot_data,
      aes(dot_x, dot_y, fill = aoi_group, alpha = dot_alpha),
      shape = 21, size = 0.78, stroke = 0.07, colour = "white"
    ) +
    # 3) 叠加原始距离的箱线、均值与记录数标签。
    # 细胡须、白箱、中位数线均来自原始米制统计，不删除真实离群值。
    geom_segment(
      data = box_data,
      aes(x = whisker_low_m_view, xend = whisker_high_m_view,
          y = group_id, yend = group_id),
      linewidth = 0.38, colour = ink
    ) +
    geom_segment(
      data = box_data,
      aes(x = whisker_low_m_view, xend = whisker_low_m_view,
          y = group_id - 0.055, yend = group_id + 0.055),
      linewidth = 0.38, colour = ink
    ) +
    geom_segment(
      data = box_data,
      aes(x = whisker_high_m_view, xend = whisker_high_m_view,
          y = group_id - 0.055, yend = group_id + 0.055),
      linewidth = 0.38, colour = ink
    ) +
    geom_rect(
      data = box_data,
      aes(xmin = q25_m_view, xmax = q75_m_view,
          ymin = group_id - 0.075, ymax = group_id + 0.075),
      fill = "white", colour = ink, linewidth = 0.38
    ) +
    geom_segment(
      data = box_data,
      aes(x = median_m_view, xend = median_m_view,
          y = group_id - 0.075, yend = group_id + 0.075),
      colour = ink, linewidth = 0.55
    ) +
    geom_point(
      data = box_data, aes(mean_m_view, group_id),
      shape = 23, size = 2.25, fill = "#B9535F",
      colour = "white", stroke = 0.40
    ) +
    geom_text(
      data = box_data,
      aes(x = 0, y = group_id - 0.54, label = count_label),
      family = font_family, hjust = 0, size = 2.65, colour = "#647681"
    ) +
    # 4) 标注均值差、逐项95%区间与统一校正后的显著性。
    geom_segment(
      data = test_data,
      aes(x = bracket_x, xend = bracket_x, y = 1, yend = 2),
      linewidth = 0.40, colour = "#50616B"
    ) +
    geom_segment(
      data = test_data,
      aes(x = bracket_x - 0.025 * maximum_view, xend = bracket_x,
          y = 1, yend = 1),
      linewidth = 0.40, colour = "#50616B"
    ) +
    geom_segment(
      data = test_data,
      aes(x = bracket_x - 0.025 * maximum_view, xend = bracket_x,
          y = 2, yend = 2),
      linewidth = 0.40, colour = "#50616B"
    ) +
    geom_text(
      data = test_data,
      aes(x = bracket_x + 0.025 * maximum_view, y = 1.5, label = significance),
      size = 3.5, fontface = "bold", colour = "#A13C4C", hjust = 0
    ) +
    geom_text(
      data = test_data, aes(x = 0, y = 2.65, label = result_label),
      family = font_family, hjust = 0, size = 2.85, colour = ink
    ) +
    # 5) 分区面板共用坐标范围；横轴标签仍以米显示。
    facet_wrap(
      ~ district_name, ncol = 2, scales = "fixed",
      labeller = labeller(
        district_name = setNames(
          paste(seq_along(district_order), district_order), district_order
        )
      )
    ) +
    scale_fill_manual(values = distance_colours, guide = "none") +
    scale_alpha_identity() +
    scale_x_continuous(
      breaks = distance_transform$transform(break_m),
      labels = scales::label_number(big.mark = ",")(break_m),
      limits = c(0, maximum_view * 1.16), expand = expansion(mult = c(0, 0))
    ) +
    scale_y_continuous(
      breaks = c(1, 2), labels = c("AOI 外", "AOI 内"),
      limits = c(0.36, 2.80), expand = expansion(mult = c(0, 0))
    ) +
    labs(
      x = paste0("到最近", facility_name, " POI 的直线距离（米；压缩刻度）"),
      y = NULL
    ) +
    theme_classic(base_family = font_family, base_size = 11) +
    theme(
      strip.background = element_blank(),
      strip.text = element_text(face = "bold", colour = ink, size = 11.5),
      axis.text = element_text(colour = ink, size = 8.5),
      axis.text.y = element_text(face = "bold", size = 10),
      axis.title.x = element_text(colour = ink, size = 10.5, margin = margin(t = 8)),
      axis.line.y = element_blank(), axis.ticks.y = element_blank(),
      axis.line.x = element_line(colour = "#70838C", linewidth = 0.30),
      axis.ticks.x = element_line(colour = "#70838C", linewidth = 0.30),
      panel.grid.major.x = element_line(colour = "#EDF1F3", linewidth = 0.24),
      panel.spacing = grid::unit(7, "mm"),
      plot.background = element_rect(fill = "white", colour = NA),
      plot.margin = margin(10, 16, 10, 10)
    )
}

17 · 绘制公交与地铁距离图

# 模块 17 · 绘制公交与地铁距离图 ----------------------
# 调用上一步的 distance_plot()；分别展示并保存两类设施的距离图。

bus_distance_plot <- distance_plot("公交")
subway_distance_plot <- distance_plot("地铁")
export_board(
  bus_distance_plot, "weibo-bus-distance-six-districts-2020-2025",
  width = 13.6, height = 11.6
)
export_board(
  subway_distance_plot, "weibo-subway-distance-six-districts-2020-2025",
  width = 13.6, height = 11.6
)
show_course_plot(bus_distance_plot)
show_course_plot(subway_distance_plot)

18 · 保存分析结果

# 模块 18 · 保存分析结果

# RDS不直接序列化terra外部指针；密度面以GeoTIFF另存，便于独立核对。
# 1) 保存分析对象；栅格已在 KDE 模块另存。
saveRDS(
  list(
    districts = study_districts,
    weibo = weibo_study,
    aoi = aoi_regions,
    district_summary = district_summary,
    bandwidth_summary = bandwidth_summary,
    aoi_summary = aoi_summary,
    distance_summary = distance_summary
  ),
  file.path(result_dir, "analysis_results.rds")
)

本次结果：地铁方向一致，公交存在分区差异

地铁：六区AOI内记录到最近地铁POI的平均距离均较短，AOI外减内的均值差约为25.7–444.2米，12项统一Holm校正后均达到0.05阈值。
公交：黄浦、徐汇、静安的AOI内平均距离较短，差值约为8.5–25.6米；长宁、杨浦的方向相反，AOI外分别短约131.4米、28.2米，差异均达到校正后的0.05阈值。闵行未达到该阈值。

这些结果支持样本记录热点与地铁设施的邻近性存在一致方向的差异，公交关系则因区而异。现有比较没有设置“远近”的统一判定标准，也没有检验全部微博位置相对于城市其他位置是否更靠近站点，因此还不能概括为“微博都靠近公共交通”或“AOI可达性普遍更高”。

差异方向与量级：下表的均值差均为“AOI外−AOI内”，单位为米；正值表示AOI内平均更近，负值表示AOI外平均更近。$p$值需要结合差异米数及区间阅读，不能据其大小判断出行意义。

研究区	公交均值差 m	公交Holm p	地铁均值差 m	地铁Holm p
1 · 黄浦区	+25.6	<0.001	+40.1	<0.001
2 · 徐汇区	+8.5	<0.001	+98.7	<0.001
3 · 长宁区	-131.4	<0.001	+25.7	0.018
4 · 静安区	+14.5	<0.001	+228.9	<0.001
5 · 杨浦区	-28.2	<0.001	+61.7	<0.001
6 · 闵行区	-5.0	0.226	+444.2	<0.001

比较时按这个顺序阅读

先看样本量。 从图中的记录数、用户数标签及分组统计表判断组别规模，留意同一用户的重复记录。点阵只展示子样本，不用于比较总量。
再看中间分布。 比较中位数与四分位距，判断典型距离和离散程度；半密度轮廓辅助观察形状。
接着检查长尾。 查看差异集中在哪段距离，并核对统计表中的500米内比例。中位数较短，并不保证这一比例更高。
最后读均值差与区间。 红色菱形为原始距离均值；均值差为AOI外减AOI内，正值表示AOI内平均更近。结合差异米数、置信区间与Holm校正结果判断。

星号比较的是原始距离均值，不是中位数。 均值差为正表示AOI外记录平均更远，为负表示AOI外平均更近；白色箱体中线仍表示中位数，两者方向可能不同。*、**、*** 分别对应校正后 $p<0.05$、$p<0.01$、$p<0.001$，ns 表示未达到0.05阈值。

公交与地铁分别比较；更换带宽或AOI阈值后，需重新核对分组统计与差异方向。

查看六区分组统计：中位距离与500米内比例

研究区	设施	分组	记录数	中位距离 m	500米内 %
黄浦区	公交	AOI 内	10172	80.3	100.0
黄浦区	公交	AOI 外	19728	103.8	100.0
黄浦区	地铁	AOI 内	10172	278.1	94.7
黄浦区	地铁	AOI 外	19728	284.0	84.7
徐汇区	公交	AOI 内	8711	152.0	100.0
徐汇区	公交	AOI 外	14278	132.3	99.7
徐汇区	地铁	AOI 内	8711	238.5	88.8
徐汇区	地铁	AOI 外	14278	366.2	73.5
长宁区	公交	AOI 内	4213	143.5	71.5
长宁区	公交	AOI 外	8844	145.5	99.9
长宁区	地铁	AOI 内	4213	403.0	58.5
长宁区	地铁	AOI 外	8844	471.8	54.4
静安区	公交	AOI 内	7793	92.6	100.0
静安区	公交	AOI 外	11895	120.2	100.0
静安区	地铁	AOI 内	7793	123.9	95.7
静安区	地铁	AOI 外	11895	366.8	69.6
杨浦区	公交	AOI 内	3320	141.6	89.5
杨浦区	公交	AOI 外	5840	152.3	99.3
杨浦区	地铁	AOI 内	3320	322.3	77.1
杨浦区	地铁	AOI 外	5840	473.1	53.2
闵行区	公交	AOI 内	6815	131.4	91.7
闵行区	公交	AOI 外	10960	162.7	98.1
闵行区	地铁	AOI 内	6815	437.6	52.4
闵行区	地铁	AOI 外	10960	737.8	29.7

检验方法、完整结果与解释边界

比较目标是本样本内按记录加权的平均距离差，方向为“AOI外−AOI内”。同一账号的记录作为一个统计聚类，采用CR1聚类稳健标准误，以用户聚类数减1为自由度进行双侧t检验。六区的公交、地铁共12项比较统一采用Holm校正；图中的区间为逐项95%置信区间，未作同时区间校正。

研究区	设施	用户聚类数	均值差 m	95%区间 m	原始 p	Holm p	标记
黄浦区	公交	24902	25.6	[24.2, 27.0]	< 0.001	<0.001	***
黄浦区	地铁	24902	40.1	[34.8, 45.5]	< 0.001	<0.001	***
徐汇区	公交	18587	8.5	[5.2, 11.8]	< 0.001	<0.001	***
徐汇区	地铁	18587	98.7	[91.2, 106.3]	< 0.001	<0.001	***
长宁区	公交	10530	-131.4	[-142.8, -120.0]	< 0.001	<0.001	***
长宁区	地铁	10530	25.7	[6.4, 45.0]	0.00902	0.018
静安区	公交	16132	14.5	[11.6, 17.4]	< 0.001	<0.001	***
静安区	地铁	16132	228.9	[217.0, 240.8]	< 0.001	<0.001	***
杨浦区	公交	7168	-28.2	[-35.4, -21.1]	< 0.001	<0.001	***
杨浦区	地铁	7168	61.7	[37.5, 85.9]	< 0.001	<0.001	***
闵行区	公交	13296	-5.0	[-13.0, 3.1]	0.22617	0.226	ns
闵行区	地铁	13296	444.2	[410.8, 477.6]	< 0.001	<0.001	***

这项检验的解释范围

用户聚类处理同一账号内的相关性，不保证消除不同用户之间的空间相关。
当前检验将AOI边界视为固定，没有传播“由本样本估计AOI”的不确定性；后续可用独立时期验证，或在重采样中重新估计AOI。
记录加权的均值受重复发帖及真实长尾影响，需结合中位数、500米内比例及去重或分时段的敏感性分析判断。

本次属于探索性的组间比较。显著差异不等于交通设施造成热点集聚，也不能直接解释为居民可达性或服务公平性的差异。

5.0 专题研究 2：回顾与反思

5.1 研究回顾：问题与回答

研究问题	本例的回答	对应证据
微博记录主要集中在哪里？	各区识别出3–5个连通AOI；约5%的区内面积包含32.27%–39.58%的样本记录。	原始点、KDE表面、AOI范围及记录覆盖。
AOI内外的设施距离有何差异？	地铁六区均为AOI内平均更近。公交在黄浦、徐汇、静安为内近，长宁、杨浦为外近；闵行差异未达到校正后的0.05阈值。	距离分布、均值差与检验表。
改变参数或时段后，发现是否保持？	尚未完成这组比较，当前不能判断边界与距离差异的稳定性。	见下方设置与改进。

分析对象如何变化？ 全市4%样本215,968条 → 上海边界内215,403条 → 六区记录112,569条 → 六份KDE表面 → 各区高值AOI → 比较AOI内外记录到全市交通POI的距离。

同一区也要区分设施类型。 长宁AOI内记录到公交POI的平均距离比AOI外记录长约131.4米，到地铁POI则比AOI外记录短约25.7米。将两类设施合成“公共交通更近”的判断，会掩盖这种差异。本例可以讨论记录热点的设施邻近性；解释热点形成或出行便利，还需其他资料。

5.2 知识回顾：从记录到高密度区

原始点、KDE、AOI与设施距离分别回答什么？

原始点显示记录在哪里，KDE估计局部记录强度，AOI按明确规则选取高值范围，设施距离比较再描述这些范围内外的邻近性。每一步都增加了不同的分析设定。

区内“记录数÷面积”概括整体平均强度，KDE进一步表示局部差异；两者都是2020–2025年的累计样本记录强度。

为什么相同分位规则，记录覆盖比例却不同？

本例按每区有效网格密度的第95百分位数提取AOI。阈值规定密度值的排名位置，没有预先规定AOI中包含多少条记录。 下面选取六区中记录覆盖比例最低和最高的两区对照：

研究区	AOI面积占本区比例	AOI包含的样本记录比例
长宁	5.00%	32.27%
静安	5.00%	39.58%

两区的面积占比四舍五入后相同，记录覆盖却不同。这反映当前平滑与阈值设定下的集中程度，不能直接判断哪一区城市活力更高。分位数取值、并列密度、边缘格的纳入和几何裁剪，也使实际面积不必恰好等于区面积的5%；六区本次约为4.981%–5.006%。

密度高值、AOI与显著热点，各有判定依据

KDE颜色表示平滑强度，AOI是按阈值划出的区域。它们没有进行空间显著性检验。本例的95%分位阈值不对应显著性水平0.05，距离图中的星号也不检验AOI是否为显著热点。

哪些设置只改图，哪些设置会改变结果？

设置	在本例中改变什么？	比较时怎样处理？
KDE带宽	改变平滑尺度，进而影响高值位置、AOI及内外分组。	从密度估计开始重算。
75米评价网格	改变取值位置与离散近似，可能影响阈值和边界细节。	固定带宽作网格细度对照。
AOI分位阈值	改变从同一表面截取的范围及记录分组。	重新提取AOI，再比较距离。
色盘、透明度、显示投影	改变图面呈现；分析仍使用原UTM密度、面积与距离。	保留原数值，检查图例与显示范围。

各区使用相同的自动带宽规则，实际带宽仍不同：黄浦X／Y方向约134／171米，闵行约635／1,050米。独立色标和独立视野也会影响观感。比较相同尺度的局部强度时，可另做共同带宽、共同色标的对照；评价局部细节时，则需考虑各区的数据量和分布。

观测单位、权重与显著性怎样一起读？

本例比较原始米制距离的记录加权均值差，方向为“AOI外−AOI内”。正值表示AOI内平均更近，负值表示AOI外平均更近。中位数、均值和500米内比例分别描述典型位置、算术平均与阈值内覆盖，不能相互代替。

记录多的用户，在均值中仍有较大权重。 以单组内两名用户分别留下10条、1条记录为教学例，按记录平均时权重为10∶1；若先计算每名用户的组内平均，再将用户等权平均，权重才是1∶1。本例采用前者，用户聚类稳健标准误处理的是同一账号内记录的相关性，不会自动把估计量改为用户等权。

显著性与实际意义一起判断。 徐汇公交均值差约8.5米，Holm校正后仍显著；是否有实际出行意义，要结合步行路线、入口和服务条件判断。闵行公交差约−5.0米，逐项95%置信区间为[−13.0, 3.1]米，未达到校正后的0.05阈值；这也不能证明两组距离完全相同。

12项比较统一作Holm校正，图中的置信区间仍是逐项95%区间，不能解释为12个区间同时覆盖的保证。图面压缩长尾、调整颜色或抽取展示点，都不会替代这些统计判断。

5.3 研究设置的不足与改进

先把“可能不稳定”变成可以复现的对照

目前已完成六区主方案。以下均为尚未运行的敏感性分析设计，沿用六区研究范围，每次改变一项设置，分别判断空间范围与设施距离结论如何变化。

改变什么？	固定什么、怎样比较？	主要判断依据
带宽	各区当前X／Y实际带宽分别乘0.5、1、2；固定记录、75米网格、区界及设施。每份密度表面按本区P95规则重新取阈值。	AOI位置与连通块、记录归组变化，以及公交／地铁均值差的方向与米数。
AOI阈值	固定同一份KDE，对照P90、P95、P97.5；重新划界、分组和比较距离。	范围扩张或收缩后，记录／用户覆盖及距离差是否保持。
贡献权重	先固定主方案AOI和原距离，比较记录加权与“每名用户在各区、各AOI组内先求均值，再将用户等权平均”。	高频账号的记录贡献是否影响均值差；两种估计量分别代表什么。

带宽比较看位置，不能只看面积。 各方案仍取P95，AOI面积本来就接近本区的5%。可在相同评价范围内计算交并比（Jaccard index）：两个AOI的交集面积除以并集面积，越接近1，空间重合越多。再检查有多少条相同记录改变了内外分组，避免仅凭外观判断。

带宽试验可在模块1将 bandwidth_multiplier 设为 0.5、1 或 2，再从模块8顺序重跑KDE及后续分析。该倍率同时进入 kde2d(h = ...) 和带宽报告；effective_h_m 是报告实际高斯尺度的对象，不作为独立调参入口。

阈值比较看覆盖与效应。 同一表面的P90、P95和P97.5高值集合天然嵌套，较高阈值选出较小范围。这种包含关系来自规则；判断结论是否敏感，应结合记录归组、均值差和区间随阈值的变化。

权重比较先明确研究对象。 用户等权改变平均值的构成；用户聚类稳健标准误处理账号内相关性，两者作用不同。同一用户可在AOI内外都留有记录，分组汇总不自动构成配对检验。若另按用户贡献重新估计KDE，则属于另一项完整流程比较。

固定共同参照。 记录输入样本、时段、区界、评价网格、设施版本和基准参数；每项对照列出记录数及用户数。
重算受影响的环节。 改带宽需重跑KDE、AOI与分组；改阈值可复用原密度表面；只比较距离的贡献权重时保持AOI不变。
合看空间与效应。 每区同时报告边界重合、连通数量、记录归组改变率、用户覆盖，以及两类设施的均值差和区间。说明各方案的估计量与推断口径，参数选择不以星号增多为目标。

网格细度可作为数值近似检查：保持实际带宽和评价范围，另比50、75、100米网格。它改变评价位置与边界近似，不等同于改变核的平滑尺度。

样本、时间与行政边界分别限制什么？

样本与选择: 4%抽样针对已获取的记录，平台使用、地理标记与数据收录的选择仍然存在。改变抽样或贡献权重时，需分别记录数据组成与自动带宽的变化；将图例乘25不能恢复完整数据的热点形态。
时间匹配: 在时长和采集覆盖可比的时段重复分析，可检查累计AOI的持续性。先固定带宽、网格与分位规则，观察记录分布的变化；若再逐期自动选带宽，应另作一组工作流比较。只有一版设施时，各期距离表示相对于同一设施参照的变化，不应解释为发帖当年的交通条件。
区界与边缘: 逐区截点截断了邻区记录的平滑贡献。本例未做边界质量补偿，区内密度积分不必等于输入记录数。可保持基准带宽，在共同区内网格上对照是否纳入邻区缓冲记录；缓冲提供额外观测，边界修正调整核质量，两种处理应分别说明。

AOI由同一批记录估计，不确定性还缺哪一层？

当前推断将AOI视为固定。用户聚类处理账号内相关性，Holm校正处理12项比较的多重性；两者均没有传播重新估计AOI的变动，也未专门处理不同用户之间的空间相关。

进阶可用一个时段识别AOI，在未参与划界的另一个时段评价距离。若进行重抽样，应按研究目的处理用户或空间依赖，并在每轮重新估计KDE与AOI；固定边界后仅反复抽记录，描述的是边界给定条件下的变动。这些设计尚未实施。

5.4 概念拓展：从平台记录到城市问题

5.4.2 城市活力：记录多，还需要知道什么？

城市活力关注城市空间中人的活动及其组织状况，涉及强度、时间节律、活动类型和社会交往。本专题从活动记录的空间强度与集中程度切入，考察集中位置与公共交通的联系。

活力的观察方面	本专题提供什么？	进一步需要什么？
活动强度	记录数÷区面积，以及KDE局部强度。	相同时段的客流或现场活动核对。
参与范围	记录数与匿名账号数。	检查账号贡献，明确平台覆盖。
空间集中	AOI面积及记录覆盖比例。	结合场所用途解释集中位置。
时间与活动多样性	当前为多年累计位置记录。	分时分析、停留与活动类型观测。

例如，静安AOI约占区面积的5%，包含约39.58%的样本记录；长宁相应约为5%与32.27%。在本次设置下，这说明两区记录的空间集中程度不同。AOI面积约5%主要来自提取规则，不能据此给六区评定活力等级。

可以从六区AOI及AOI外选取调查地点，在相同季节和时段记录步行人数、停留与活动类型，检验微博指标覆盖了实际空间使用的哪些方面。概念与测量的核心阅读为Liu等（2023），见文末。

5.4.3 城市活动热区：相对哪里、按什么规则判定？

活动热区是在明确范围、时段、资料和规则下，活动记录相对集中的位置。本例从KDE表面按各区P95提取AOI，因此其含义是区内相对高值。

热区的解释须带上评价范围与观察尺度：共同带宽用于比较同一尺度，区内分位用于定位各区内部的集中位置。边界是否稳定按5.3的通用方案检验；与统计显著热点的区别见5.2。

5.4.4 空间可达性：给定条件下能获得哪些机会？

本专题已计算微博记录到最近公交或地铁POI的直线距离，并比较AOI内外记录的距离分布。这个指标描述设施邻近性，适合回答“记录集中位置是否更靠近站点”。当前结果中，六区地铁均为AOI内平均更近，公交差异的方向则因区而异。

空间可达性关心在特定交通方式、出发时刻和使用条件下，从一个地点能够获得哪些活动机会。到站距离是连接服务的一个环节，入口、候车、换乘及目的地分布共同影响能到达哪里。

沿本例继续研究	可以采用的指标	需补充的资料
记录位置是否靠近站点？	最近公交／地铁POI直线距离。	本例已完成。
给定时间内能到达哪些机会？	步行到站时间、指定时限内可达岗位或服务数。	入口、路网、班次、换乘及目的地。
可获得的服务是否足够？	同时考虑服务容量、需求与出行时间。	容量、需求人口与明确的服务范围。

Cooke等（2022）讨论距离之外的使用条件，见文末。

5.4.5 空间不平等与公平：差距有多大，怎样评价？

本例已经观察到不同区、不同设施类型及AOI内外之间的距离差异。空间不平等关注资源、机会和负担在地区或人群间怎样分布；若进一步研究交通服务，可比较哪些地区或人群获得的出行机会较少、付出的时间或费用较多。

公平评价需要明确标准。 基本服务门槛、相近需求获得相近机会、对行动受限人群提供支持，分别体现不同的评价依据。若研究居民服务差距，应将六区AOI与居民及需求分布结合，再比较机会与负担；微博发帖地点不直接代表居住地点。

Verlinghieri与Schwanen（2020）可作为公平评价的一条延伸阅读。设施位置变化怎样影响受益人口，可衔接专题1的十五分钟城市阅读。这些问题需要补充人口、路网和服务资料，尚未在本例实施。

5.5 地理学方法拓展与阅读

从区内比较继续提出什么问题？

当前比较的是AOI内外两组微博记录位置。若要追问微博位置是否比城市其他位置更靠近设施，可在预先确定的可比较范围设置背景位置，再计算相同的设施距离。

背景参照	对应的问题	设计时需要明确什么？
范围内均匀背景	微博位置相对于一般地理位置有何差异？	可进入空间或陆地范围、抽样规则与数量。
人口加权背景	微博位置相对于居民分布有何差异？	人口数据、年份及空间单元；居住位置与活动位置的区别。

这是一项尚未运行的研究拓展。商业、就业和用地可能同时吸引站点与活动；背景比较提供新的参照，解释设施变化的作用仍需匹配时期和适当对照。

核心阅读：每篇对应一个研究问题

城市活力测量了什么？ Liu等（2023），Urban vitality assessment at the neighborhood scale with geo-data: A review toward implementation，Journal of Geographical Sciences, 33, 1482–1504。论文 · 公开全文。先读第3–4节的维度、时间与测量讨论：本例指标覆盖了活力的哪个方面，还需要哪种独立观测？
随机抽样后还保留哪些选择？ Petutschnig等（2021），Evaluating the Representativeness of Socio-Demographic Variables over Time for Geo-Social Media Data，ISPRS International Journal of Geo-Information, 10(5), 323。论文。关注地理社交媒体的代表性与时间变化，思考对已获取记录抽取4%后，平台使用与数据收录的选择是否改变。
怎样从累计分布走向分时研究？ Niu与Silva（2023），Understanding temporal and spatial patterns of urban activities across demographic groups through geotagged social media data，Computers, Environment and Urban Systems, 100, 101934。论文。重点参考其分时与分组思路；人口分组需要另有人群资料，不能从本例匿名账号直接得到。
距离之外，什么影响服务使用？ Cooke等（2022），Proximity is not access: A capabilities approach to understanding non-motorized transport vulnerability in African cities，Frontiers in Sustainable Cities, 4, 811049。论文与公开全文。结合其研究设计，列出直线距离相同、实际出行条件仍可能不同的两个原因。

6.0 章节练习与测评

实战目标

本节先用闵行区咖啡POI练习DBSCAN密度聚类与聚类AOI，再自主选取点数据开展参数敏感性分析，最后通过六道自测回顾本章的方法。

任务 1：跟随示例运行

数据与研究问题

研究问题：在500米邻域、至少10点的规则下，闵行咖啡POI形成哪些密度簇？哪些点未被归簇，簇点的外轮廓怎样表达？

从课程数据页下载咖啡厅POI和上海市行政区界，放入第5章的data/目录。边界文件含上海各区，代码按区代码310112筛选闵行。

本地文件	用途
`data/咖啡厅.geojson`	咖啡POI位置；保留合法点，不按同名连锁店直接去重。
`data/shanghai_boundary.geojson`	上海各区边界；先检查CRS，再按区代码筛选闵行。

DBSCAN怎样把点连接成簇？

DBSCAN（Density-Based Spatial Clustering of Applications with Noise）按局部邻域的点数识别密度簇，无需预设簇数。

核心点：半径eps的邻域内至少有minPts个点，包括自身。
边界点：自身邻域点数不足，但位于某个核心点的邻域内，可归入该簇。
噪声点：既不是核心点，也没有归入核心点邻域；R结果中的簇号为0。

相邻核心点可以逐步连接成不规则簇，边界点不负责继续扩张。邻域半径约束点之间的连接距离，不限定整个簇的半径或直径。 噪声表示当前参数下未归簇，不等于错误数据。

DBSCAN本身不提供显著性检验的p值，识别出的密度簇不能直接解释为统计显著热点。

模拟点的DBSCAN分类，核心、边界和噪声用不同颜色区分；P、Q、R的同半径邻域及真实点数帮助核对分类。

DBSCAN点类别（模拟）：相同半径下，P满足核心点条件，Q为边界点，R为噪声点；邻域计数包含自身。

先看邻域，再调整参数

eps与minPts共同定义局部密度门槛：前者决定查多远，后者决定邻域中至少有多少点。计数包括被查看的点本身。

在二维平面上，半径翻倍，邻域圆面积变为4倍。点数门槛不变时，扩大邻域会让更多点有机会达到核心点条件。

调整哪一项？	固定条件	观察什么？
扩大`eps`	保持`minPts`不变。	更多点进入邻域，部分点成为核心点；新簇可能出现，原有簇也可能连接。
提高`minPts`	保持`eps`不变。	核心点的门槛变高；部分点转为边界或噪声，原簇也可能分裂。

调整参数或切换数据，点击点查看邻域。

观察重点：哪些点改变类别，哪些点簇连接或分离？簇号与颜色仅对应当前结果；凹包在后续 AOI 步骤生成。

方法或结果	本例中回答什么问题？	输出
KDE	哪些位置的记录强度较高？	连续密度表面；带宽控制平滑尺度。
DBSCAN	哪些点通过密集邻域连接在一起？	点的簇号及噪声标记。
聚类AOI	怎样表示已识别点簇的外轮廓？	由簇内点生成的面；凹包参数控制轮廓。

分步运行DBSCAN模板

工作目录设为第5章目录，代码块按顺序运行。首次绘制瓦片需要联网，后续复用本地缓存。DBSCAN与凹包均在UTM 51N米制坐标下计算，仅将展示副本转回经纬度。

eps_m与min_pts控制点的聚类，hull_ratio控制已识别点簇的外轮廓。修改聚类参数后，从步骤3起顺序重跑；只改变凹包参数时，重跑步骤6、7。

步骤 1：环境与参数

# 1. 加载分析与制图包；缺少的包先在Console安装。
# install.packages(c("dbscan", "sf", "tidyverse", "ggspatial", "rosm"))
library(dbscan)
library(sf)
library(tidyverse)
library(ggspatial)

# 2. 设置聚类和凹包参数。
eps_m <- 500
min_pts <- 10
hull_ratio <- 0.8

# 3. 注册底图源，并保存瓦片缓存。
esri_light_gray_url <- paste0(
  "https://server.arcgisonline.com/ArcGIS/rest/services/Canvas/",
  "World_Light_Gray_Base/MapServer/tile/${z}/${y}/${x}"
)
rosm::register_tile_source(
  esri_light_gray = rosm::as.tile_source(esri_light_gray_url, extension = "jpg")
)
dir.create(file.path("cache", "map_tiles"), recursive = TRUE, showWarnings = FALSE)

步骤 2：读入、裁取与投影

# 1. 读入区界，只保留区代码、名称与几何。
shanghai_boundary <- st_read("data/shanghai_boundary.geojson", quiet = TRUE) %>%
  select(district_code = "县代码", district_name = "县名") %>%
  st_transform(4326) %>%
  st_make_valid()
minhang_boundary <- shanghai_boundary %>% filter(district_code == 310112)

# 2. 读入咖啡POI，裁取闵行范围；保留原始行序供重复聚类使用。
poi_coffee <- st_read("data/咖啡厅.geojson", quiet = TRUE) %>%
  st_transform(4326)
minhang_cafe_input <- st_filter(
  poi_coffee, minhang_boundary, .predicate = st_intersects
)

# 3. 转为UTM 51N米制坐标，提取DBSCAN的输入矩阵。
cafe_for_clustering <- st_transform(minhang_cafe_input, 32651)
coords <- st_coordinates(cafe_for_clustering)[, c("X", "Y"), drop = FALSE]

步骤 3：聚类与点类别

# 1. 识别密度簇和核心点；minPts包含点本身。
db_res <- dbscan::dbscan(coords, eps = eps_m, minPts = min_pts)
core_point <- dbscan::is.corepoint(coords, eps = eps_m, minPts = min_pts)

# 2. 按原始行序写回结果，再按簇号整理显示顺序。
minhang_cafe <- minhang_cafe_input %>%
  mutate(
    cluster_id = db_res$cluster,
    point_type = case_when(
      core_point ~ "核心点",
      cluster_id > 0 ~ "边界点",
      TRUE ~ "噪声点"
    ),
    cluster_label = if_else(
      cluster_id == 0, "Noise", paste0("Cluster ", cluster_id)
    )
  ) %>%
  arrange(cluster_id)

# 簇号0表示噪声，其余编号用于绘图。
cluster_ids <- sort(unique(db_res$cluster[db_res$cluster > 0]))

步骤 4：分类结果汇总

下表汇总核心点、边界点与噪声点的数量。比较不同参数时，按空间位置和点成员核对，簇号只表示当前分组。

# 按点类别汇总数量与比例。
point_summary <- minhang_cafe %>%
  st_drop_geometry() %>%
  count(point_type, name = "n_points") %>%
  mutate(share = sprintf("%.1f%%", 100 * n_points / sum(n_points)))

knitr::kable(
  point_summary, format = "html",
  caption = paste0("点类别与数量（", length(cluster_ids), "个簇）"),
  col.names = c("点类别", "点数", "比例"), align = c("l", "r", "r")
)

点类别与数量（9个簇）
点类别	点数	比例
噪声点	381	66.7%
核心点	138	24.2%
边界点	52	9.1%

步骤 5：绘制点分类地图

$闵行区咖啡POI密度聚类地图，颜色区分当前参数下的簇，灰色十字表示未归簇点。$

闵行咖啡POI的DBSCAN密度簇；参数见图中副标题，噪声点单独显示。

查看并复制点图R代码

# 图例映射与点图
# 按簇号设置离散色盘；噪声单独显示。
cluster_cols <- setNames(
  grDevices::hcl.colors(length(cluster_ids), palette = "Dark 3"),
  sprintf("Cluster %s", cluster_ids)
)

# 设定噪声点颜色，并合并为统一的颜色映射表
plot_cols <- c("Noise" = "grey20", cluster_cols)

# 设定形状映射：噪声点为十字，聚类点为实心圆
plot_shapes <- c(
  "Noise" = 3,
  setNames(rep(16, length(cluster_ids)), names(cluster_cols))
)

# 设定大小映射：噪声点较小，聚类点相对更突出
plot_sizes <- c(
  "Noise" = 0.8,
  setNames(rep(2.5, length(cluster_ids)), names(cluster_cols))
)

# 5. 可视化：闵行区咖啡店 DBSCAN 聚类结果
# 使用本模板已注册的Esri浅灰底图；瓦片缓存可供后续复用。

ggplot() +
  # 5.1 底图图层
  annotation_map_tile(
    type = "esri_light_gray",
    zoom = 12,
    alpha = 0.40,
    cachedir = file.path("cache", "map_tiles"),
    quiet = TRUE
  ) +
  
  # 5.2 行政边界图层
  geom_sf(
    data = minhang_boundary,
    aes(linetype = "Minhang Boundary"),
    fill = NA,
    colour = "#496672",
    linewidth = 0.8,
    alpha = 0.8
  ) +
  
  # 5.3 点要素图层
  # 同时映射颜色、形状与大小，以保证图面与图例表达一致
  geom_sf(
    data = minhang_cafe,
    aes(
      color = cluster_label,
      shape = cluster_label,
      size = cluster_label
    ),
    alpha = 0.85
  ) +
  
  # 5.4 制图要素：比例尺与指北针
  annotation_scale(
    location = "bl",
    width_hint = 0.25,
    style = "ticks",
    text_cex = 0.6,
    text_col = "grey20"
  ) +
  annotation_north_arrow(
    location = "br",
    which_north = "true",
    height = grid::unit(0.8, "cm"),
    width = grid::unit(0.8, "cm"),
    pad_x = grid::unit(0.1, "cm"),
    pad_y = grid::unit(0.1, "cm"),
    style = north_arrow_orienteering(
      fill = c("white", "#253e49"),
      line_col = "#253e49",
      line_width = 0.6,
      text_col = "#253e49",
      text_size = 8
    )
  ) +
  
  # 5.5 图例与标度设置
  scale_linetype_manual(
    name = "Study Area",
    values = c("Minhang Boundary" = "solid")
  ) +
  scale_color_manual(
    name = "DBSCAN Result",
    values = plot_cols
  ) +
  scale_shape_manual(
    name = "DBSCAN Result",
    values = plot_shapes
  ) +
  scale_size_manual(
    name = "DBSCAN Result",
    values = plot_sizes
  ) +
  
  # 调整图例顺序，并覆盖边界图例的默认显示样式
  guides(
    linetype = guide_legend(
      override.aes = list(color = "#496672", linewidth = 0.8),
      order = 1
    ),
    color = guide_legend(order = 2),
    shape = guide_legend(order = 2),
    size = guide_legend(order = 2)
  ) +
  
  # 5.6 版式设置
  theme_bw() +
  theme(
    legend.position = "right",
    legend.title = element_text(size = 10, face = "bold"),
    legend.text = element_text(size = 9),
    plot.title = element_text(face = "bold", hjust = 0.5),
    plot.subtitle = element_text(color = "grey40", hjust = 0.5, size = 10),
    panel.grid.major = element_line(colour = "grey90", linetype = "dashed")
  ) +
  
  # 5.7 图题与附注
  labs(
    title = "Coffee Shop Clusters in Minhang District, Shanghai",
    subtitle = sprintf("DBSCAN: eps = %s m; minPts = %s", eps_m, min_pts),
    caption = paste(
      "Tiles: Esri, HERE, Garmin, (c) OpenStreetMap contributors,",
      "and the GIS user community"
    )
  )

读图顺序：先看簇的位置与点数，再看未归簇点的分布；连续核心点怎样连接，可回到上方交互模块查看。500米／10点基准下，571点中有190点归入9个簇，381点为噪声，约占66.7%。已归簇点主要分布在研究区北部与中部，南部和东南部较多点未归簇；这些位置没有达到当前密度门槛。交通枢纽、办公与居住功能可作为后续核查线索。

步骤 6：生成聚类AOI

按簇汇总点，再计算凹包，得到每个簇的外轮廓。该AOI表示算法选出的点集范围；识别实际商圈还需经营联系、客流或调查资料。

# 1. 提取已归簇点，在UTM 51N米制坐标下计算轮廓。
cluster_points_projected <- minhang_cafe %>%
  filter(cluster_id != 0) %>%
  st_transform(32651)

# 2. 按簇合并点，再计算凹包；hull_ratio只控制轮廓形态。
aoi_geometries_projected <- cluster_points_projected %>%
  group_by(cluster_label) %>%
  summarise(geometry = st_union(geometry), .groups = "drop") %>%
  st_concave_hull(ratio = hull_ratio, allow_holes = FALSE)

# 3. 只保留能够成面的簇；共址或共线点集可能只形成点或线。
aoi_polygons_projected <- aoi_geometries_projected %>%
  filter(st_geometry_type(.) %in% c("POLYGON", "MULTIPOLYGON"))

# 米制副本保留供分析，仅将地图显示副本转回经纬度。
# 当前参数无簇时，各对象均返回空sf，不会保留上次结果。
aoi_polygons <- aoi_polygons_projected %>%
  st_transform(4326)

每次计算都重新生成AOI。没有可成面的簇时，结果为空，不会沿用上次的面。

步骤 7：绘制聚类AOI地图

Zoom: 12

闵行区咖啡POI聚类凹包地图，各颜色多边形对应同一密度簇的点集外轮廓，噪声点仍单独显示。

同一DBSCAN结果的聚类AOI；聚类与凹包参数见图中副标题。

查看并复制AOI制图R代码

ggplot() +
  # 7.1 底图图层
  annotation_map_tile(
    type = "esri_light_gray",
    zoom = 12,
    alpha = 0.40,
    cachedir = file.path("cache", "map_tiles"),
    quiet = TRUE
  ) +

  # 7.2 行政边界图层
  geom_sf(
    data = minhang_boundary,
    aes(linetype = "Minhang Boundary"),
    fill = NA,
    colour = "#496672",
    linewidth = 0.8,
    alpha = 0.8
  ) +

  # 7.3 聚类AOI多边形图层
  # 使用填充色与边框色同时区分不同聚类簇
  geom_sf(
    data = aoi_polygons,
    aes(
      fill = cluster_label,
      color = cluster_label
    ),
    alpha = 0.35,    # 保持一定透明度，便于识别底图信息
    linewidth = 0.6  # AOI 边界线宽
  ) +

  # 7.4 噪声点图层
  # 将未被聚类识别的点保留在图中，用作空间对照
  geom_sf(
    data = minhang_cafe %>% filter(cluster_id == 0),
    aes(shape = "Noise"),
    color = "grey30",
    size = 0.8,
    alpha = 0.6
  ) +

  # 7.5 制图要素：比例尺与指北针
  annotation_scale(
    location = "bl",
    width_hint = 0.25,
    style = "ticks",
    text_cex = 0.6,
    text_col = "grey20"
  ) +
  annotation_north_arrow(
    location = "br",
    which_north = "true",
    height = grid::unit(0.8, "cm"),
    width = grid::unit(0.8, "cm"),
    pad_x = grid::unit(0.1, "cm"),
    pad_y = grid::unit(0.1, "cm"),
    style = north_arrow_orienteering(
      fill = c("white", "#253e49"),
      line_col = "#253e49",
      line_width = 0.6,
      text_col = "#253e49",
      text_size = 8
    )
  ) +

  # 7.6 图例与标度设置
  scale_linetype_manual(
    name = "Study Area",
    values = c("Minhang Boundary" = "solid")
  ) +
  scale_fill_manual(
    name = "Coffee AOI Clusters",
    values = cluster_cols,
    limits = names(cluster_cols)
  ) +
  scale_color_manual(
    name = "Coffee AOI Clusters",
    values = cluster_cols,
    limits = names(cluster_cols)
  ) +
  scale_shape_manual(
    name = "Unclustered",
  values = c("Noise" = 3),
  limits = "Noise",
  breaks = if (any(minhang_cafe$cluster_id == 0)) "Noise" else character()
  ) +

  # 调整图例顺序，并统一图例显示样式
  guides(
    linetype = guide_legend(
      override.aes = list(color = "#496672", linewidth = 0.8),
      order = 1
    ),
    fill = guide_legend(
      order = 2,
      override.aes = list(alpha = 0.6)
    ),
    color = guide_legend(order = 2),
    shape = guide_legend(order = 3)
  ) +

  # 7.7 版式设置
  theme_bw() +
  theme(
    legend.position = "right",
    legend.title = element_text(size = 10, face = "bold"),
    legend.text = element_text(size = 9),
    plot.title = element_text(face = "bold", hjust = 0.5),
    plot.subtitle = element_text(color = "grey40", hjust = 0.5, size = 10),
    panel.grid.major = element_line(colour = "grey90", linetype = "dashed")
  ) +

  # 7.8 图题与附注
  labs(
    title = "Coffee POI Cluster AOIs in Minhang District",
    subtitle = sprintf("DBSCAN: eps = %s m; minPts = %s; hull ratio = %s",
        eps_m, min_pts, hull_ratio),
    caption = paste(
      "Tiles: Esri, HERE, Garmin, (c) OpenStreetMap contributors,",
      "and the GIS user community"
    )
  )

凸包与凹包怎样选？

凸包是包住全部输入点的最小凸多边形；凹包允许边界向内凹陷，也应覆盖输入点。本例hull_ratio = 0.8控制轮廓形态，取1时对应凸包；它不表示保留80%的点或面积。allow_holes = FALSE不保留内部洞，轮廓内的空白仍应结合原始点图检查。

相同点集的凸包和凹包并排对照，点位保持不变，两种轮廓均覆盖全部点，凹包显示内凹边界。

同一组模拟点的凸包与凹包：两者均覆盖全部输入点，凹包允许边界内凹；坐标仅为示意单位。

凹包表示已归簇点的外轮廓，形状受参数控制。解释时应结合原始点位与研究区边界，轮廓本身不代表真实商业边界。

步骤 8：参数对照：改变半径

以下是一组独立的半径对照：固定同一批POI、研究区与米制投影，保持minPts = 10，只改变邻域半径。下表由当前点坐标重新计算；最大簇占比的分母为全部区内点。

parameter_comparison <- lapply(c(300, 500, 700), function(candidate_eps_m) {
  candidate_fit <- dbscan::dbscan(coords, eps = candidate_eps_m, minPts = 10)
  point_labels <- candidate_fit$cluster
  cluster_sizes <- table(point_labels[point_labels > 0])
  tibble::tibble(
    eps_m = candidate_eps_m,
    min_pts = 10L,
    n_clusters = length(cluster_sizes),
    noise_percent = 100 * mean(point_labels == 0),
    largest_percent = 100 * max(c(0, cluster_sizes)) / length(point_labels)
  )
}) %>% dplyr::bind_rows()

parameter_comparison %>%
  dplyr::mutate(
    noise_percent = sprintf("%.2f%%", noise_percent),
    largest_percent = sprintf("%.2f%%", largest_percent)
  ) %>%
  knitr::kable(
    format = "html", align = c("r", "r", "r", "r", "r"),
    col.names = c("eps（米）", "minPts", "簇数", "噪声比例", "最大簇占比")
  )

eps（米）	minPts	簇数	噪声比例	最大簇占比
300	10	6	81.26%	6.13%
500	10	9	66.73%	7.71%
700	10	7	53.24%	13.31%

半径从300米增至700米时，簇数依次为6 → 9 → 7。部分原先未归簇的点达到密度门槛，形成新簇；原有簇也可能连接起来。因此，簇数不必随半径增大而持续减少。可回到上方真实POI演示，查看新簇出现、原簇合并及噪声转入簇的位置。

怎样选参数，而不只看图是否好看？

先说明希望识别的集聚尺度，再确定候选值。定位精度、设施间距与研究目的都应参与判断；噪声最少或簇数最多，不能单独作为择优标准。

先选点数门槛：结合数据规模及研究希望识别的局部集聚，设置候选minPts。
查看邻居距离：当minPts = 10时，用dbscan::kNNdistplot(coords, k = 9)查看第9个其他邻居的距离排序曲线。转折可提示候选半径，实际数据也可能没有清楚的转折。
每次只改变一项：同时查看点的位置、簇成员与噪声比例，报告哪些结构保留、哪些结构改变。

要核查的设置	候选与固定条件	记录什么？
邻域半径`eps`	上表已比较300、500、700米；固定`minPts = 10`。	簇数、噪声比例、最大簇占比与合并／分裂位置。
点数门槛`minPts`	可在演示中比较5、10、15；固定`eps = 500`米。	核心点、边界点和噪声怎样变化。
凹包`ratio`	固定已识别的簇，再比较0.5、0.8、1。	AOI面积与轮廓；簇号和噪声比例应保持不变。

本例使用直线距离。隔河或被道路屏障分开的点仍可能很近；若问题涉及步行联系，应进一步选择适合的距离定义。边界附近的截点影响可联系专题2的通用反思。

方法阅读：dbscan官方入门可配合参数对照阅读；Hahsler等（2019）的R实现论文解释算法及其实现。sf文档用于查询凹包形态参数与版本要求。

选学：把事件发生的时间加入邻域

同一街口一小时内连续发生的事故，与相隔数月的事故，可以形成不同的研究问题。常见的ST-DBSCAN实现同时设置空间和时间阈值：两个事件只有在空间距离与时间差均满足条件时，才互为邻居，再依据点数判断核心点。例如，空间阈值为300米、时间阈值为30分钟。

空间阈值须与投影坐标同单位，时间阈值须与输入时间同单位；先统一时区，再将时间换算成明确的秒或分钟。阈值约束相邻事件，链式连接后的簇仍可持续更长时间。

Birant与Kut（2007）的原文还讨论了非空间属性等扩展。双阈值软件实现与原算法的全部细节应分别核对。阅读：原始论文；stdbscan官方文档，查看eps_spatial、eps_temporal与min_pts。

选学：地块邻域与不同密度下的聚类

DBSCAN的邻域思想也可用于面对象或不同密度的点集。更换对象或算法时，需要重新说明邻域关系与密度判据。

方法	怎样用于其他研究？	阅读
GDBSCAN	将邻域关系与密度判据推广到空间面及属性，可研究地块群；相交、边界间距、属性相似性需按问题选定。	Sander等（1998）原文包含二维多边形应用。
HDBSCAN	在多个密度层级中提取较稳定的簇，可探索商业中心与稀疏外围设施并存的点集，仍需选择参数并检查结果。	Campello等（2013）论文；R文档。
OPTICS	通过可达距离排序展示密度结构，可检查主商业带与支路小簇；输出簇号还需选择提取规则。	Ankerst等（1999）论文；R文档。

地块的几何邻近与社会经济属性相似回答不同问题。按属性得到同类地块后，若研究要求连续片区，还须检查连通性；将地块换成中心点会改变邻域含义。

任务 2：小组实践与展示

选择一份空间点数据，自定研究问题与分析范围，使用本章学过的方法完成一份本地 R Markdown 报告。练习重点是参数敏感性分析：关键参数改变后，主要空间模式是否仍然成立？

练习要求

说明数据来源、记录含义、时间与研究范围，完成坐标系检查及必要的数据清理。
自主选择分析方法，说明主要参数的含义与取值理由。在合理范围内改变其中一项关键参数，保持其余条件一致，观察结果的变化。
Knit 生成可在浏览器中打开的本地 HTML，用图表展示主要发现及参数调整后的结果，解释哪些模式比较稳定、哪些对参数较敏感。分析代码与文字由小组自行组织。

参数可以怎样调整？

根据所选方法，任选一项开展试验。每次记录参数值，并用相同的地图范围和可比的图例呈现结果。

所选方法	可调整的一项设定	观察重点
六边形聚合	网格面积或网格起点	哪些点重新归格，高低值分布是否改变
KDE／AOI	带宽或 AOI 提取阈值	高值区的位置、平滑程度或边界范围
DBSCAN	`eps` 或 `minPts`	簇的连接与分离、归簇点与噪声点的变化
Moran’s I／LISA	邻接规则或空间权重设定	全局统计量及局部聚类位置是否稳定

结果解释需联系参数所代表的尺度或邻域含义。

备选数据与选题提示

课程数据页中的商业 POI、公共交通设施及本章已有数据均可使用，也可从下列官方平台自行选取数据。选题提示只提供起点，研究范围、方法与参数由小组确定。

数据入口	可研究的问题	数据处理提示
上海 POI 与交通设施	某类设施集中在哪里？活动点与交通设施的距离有什么特点？	可沿用本章数据；直线距离用于描述设施邻近性
纽约交通碰撞	行人或自行车碰撞集中在哪些位置？集聚范围是否随参数改变？	有事件 ID、日期、坐标及伤亡字段；按时期与事件类型筛选，检查缺失坐标
纽约市现行树木点：ForMS	树木数量、胸径或状况在空间上有什么差异？	持续更新的管理记录；对象与筛选口径见下方“下载与预读”
英国警方犯罪数据	某类犯罪记录集中在哪里？带宽或邻域半径变化后，哪些区域保持稳定？	可按月份与警区下载 CSV；发布坐标为匿名近似位置，定位说明见官方文档
USGS 地震目录	震中记录集中在哪里？KDE带宽或DBSCAN半径改变后，集聚范围如何变化？	固定区域、时期与震级范围；二维分析使用震中位置，深度不计入平面距离
Citi Bike 骑行与站点	出发事件集中在哪里？站点使用强度有什么空间差异？	区分站点与骑行事件；月度文件较大，可限定日期。按数据使用协议从官方获取数据

地震数据也可按六边形汇总事件数，再比较网格尺度。以上分析描述目录中已记录事件的空间分布，不能直接用于地震风险判断或预测。

纽约交通碰撞 · 下载步骤

打开官方数据页，核对数据集名称与 Data Last Updated。每行代表一次交通碰撞。
全量数据较大，可先在 Data 页按研究范围和时段筛选。
点击 Export，在 Download file 中选择 CSV，核对导出范围与行数，再点击 Download 保存文件。

预读提示：网页预览与 API 返回量

网页一次显示多少行，与导出多少行是两件事。直接读取 API 地址也可能只取得一页，需要按接口文档明确范围并分页读取。使用网页导出时，核对筛选条件与导出行数即可。

清理时按 collision_id 核对重复导入；相同坐标的不同碰撞仍保留为不同事件。

纽约交通碰撞官方数据页截图，顶部有Data页签和Export按钮，下方列出Data Last Updated。

交通碰撞官方数据页：核对数据名称与数据更新时间。

纽约交通碰撞的Export dataset窗口截图，Download file下选择CSV，显示导出行数与Download按钮。

CSV 导出窗口：检查格式与导出范围后下载。

官网截图采集于 2026-10-21。

ForMS 树木与种植位 · 下载与预读

官方渠道：Forestry Tree Points（树木点表）；Forestry Planting Spaces（种植位表）。若限定研究行道树，需要两表共同判断。

分别打开所需数据页，核对名称和 Data Last Updated，按研究范围筛选。
点击 Export，将 Export format 改为 GeoJSON，核对导出范围后下载。可分别保存为 tree_points.geojson、planting_spaces.geojson。也可使用官方树木点 GeoJSON 导出和种植位 GeoJSON 导出链接获取全表。
沿用第四章的 sf::st_read() 读取本地文件，先查看字段、记录数和坐标系，再确定筛选条件。

预读提示：树木对象与种植位置

ForMS 是持续更新的管理记录，同一位置可能先后对应退役树木和补种树木。两表记录的对象不同：

数据表	先核对什么？
树木点表	`tpstructure` 表示结构，如 `Full`、`Stump`、`Retired`；`tpcondition` 表示状况。结构为 `Full` 的记录也可能标为 `Dead`
种植位表	`pssite` 区分街道、公园等位置；空种植位不能计作树木。限定行道树时核对 `pssite == "Street"`

树木表的 plantingspaceglobalid 对应种植位表的 globalid。关联方法沿用第二章；两表都读为 sf 时，先用 st_drop_geometry() 取出种植位属性，保留树木点的几何。记录下载日期及结构、状况的筛选口径。

Forestry Tree Points官方页面截图，显示树木点数据名称、更新时间和Export按钮。

树木点表：Forestry Tree Points。

ForMS树木点表的导出窗口截图，Export format已经选为GeoJSON，显示导出行数及Download按钮。

导出时选择 GeoJSON，保留点几何。

Forestry Planting Spaces官方页面截图，显示种植位表名称、更新时间和Export按钮。

种植位表：Forestry Planting Spaces，按相同步骤导出。

官网截图采集于 2026-10-21。

英国警方犯罪数据 · 下载步骤

图示以伦敦 2026 年 3 至 4 月的犯罪数据演示下载流程；月份与警区可根据选题调整。

预读提示：匿名位置与空案件编号

公开坐标经过匿名化处理，多个事件可映射到同一参考点。Crime ID 也可能为空，例如反社会行为记录；空编号不表示同一事件，不能直接按坐标或空编号去重。合并月份时保留来源文件名，并核对有效事件编号。

八格教学示意图：在英国警方网站选择2026年3至4月、Metropolitan Police Service与City of London Police，勾选犯罪数据，生成并下载ZIP，再解压取得CSV。

英国警方犯罪数据下载示意图。图片由课程教师提供。

Citi Bike 站点 · 下载步骤

打开 Citi Bike System Data，找到 Real-Time Data，点击 Get the GBFS feed here。
在 GBFS 目录的 en 下找到 station_information，打开对应的 url。也可直接访问官方站点 JSON。
使用浏览器的“另存为”保存 .json 文件，文件名可包含下载日期。站点记录位于 data → stations，坐标字段为 lon、lat；station_id 保留为字符型。

预读提示：GBFS 的嵌套 JSON

GBFS 是共享出行实时数据的发布规范。station_information 的普通 JSON 按层级组织，站点表位于 data → stations，一行对应一个站点，可用于研究站点分布。研究骑行使用强度时，需要使用官网的月度 Trip Histories，按站点汇总出发或到达次数，同站的多次骑行应保留。

预读时可使用 jsonlite::fromJSON()，需先安装 jsonlite 包。若将读取结果命名为 station_json，取出 station_json$data$stations 后，再沿用第四章的经纬度转 sf 流程。

station_status 提供某一时刻的可用车辆、车桩等状态，可按 station_id 连接；单次状态不能表示一段时期的骑行次数。GBFS 还含纽约市外站点，分析时按研究区边界裁取。

Citi Bike官方System Data页面截图，在Real-Time Data标题下提供Get the GBFS feed here链接。

站点入口位于 Real-Time Data。

Citi Bike GBFS目录的真实截图，英文en部分列出station_information及其站点JSON地址。

GBFS 目录：查找 station_information 对应的地址。

官方station_information JSON截图，data的stations数组内包含station_id、name、lon、lat和capacity。

站点 JSON：data → stations 保存站点记录。

官网截图采集于 2026-10-21。

概念	含义	在本例中如何使用
空间权重 \(W\)	明确单元间的邻接关系及权重。	将“周围”落实为可计算的邻居集合。
全局 Moran’s I	概括整个研究区中属性值与其空间滞后的关联。	判断中位报价相近的六边形是否倾向于相邻。
LISA 与局部 Moran	局部空间关联指标（local indicators of spatial association, LISA）是一类指标；本例采用其中的局部 Moran。	区分显著的同向集聚与空间离群格。

顺序	怎样读	本图中的含义
1 · 横轴	读取该格的标准化报价。	\(z=0\) 对应362个参检格的中位报价均值；\(z>0\) 为高于这一均值。
2 · 纵轴	读取邻居标准化报价的行加权平均。	\(Wz>0\) 表示邻域平均报价高于同一参照，不要求每个邻居都偏高。
3 · 象限	对照本格与邻域的高低。	HH122格、LH56格、LL165格、HL19格；这些是正负组合，还没有筛局部显著性。
4 · 全局线与检验	查看整体方向，再看随机置换的证据。	过原点线向右上倾斜，斜率为0.620857；9,999次单侧置换的\(p=0.0001\)。

比较项	正六边形	正方形
完整密铺内部的邻居	6 个共边邻居，中心距离相等。	4 个共边邻居；加入角邻居后为 8 个，角邻居距离是边邻居的 \(\sqrt{2}\) 倍。
形状与数据衔接	同面积时更紧凑，适合比较六个方向的邻域，仍有固定方向。	便于与栅格像元、行列索引衔接，也可旋转。
本例需要核对	少样本格退出后，邻居可能不足 6 个；实测 Queen 与 Rook 列表相同。	Queen 与 Rook 通常不同，需随几何重新定义。

第五章： 专题 1： 点数据分析

章节作者：刘贇喆 | 华东师范大学 社会发展学院

更新时间：2026-10-26

本章学习内容

学习导航

1.0 空间点数据

1.0 空间点数据概论

1.1 常见的空间点记录

1. 地点与设施点（POI）

2. 行为与事件点

3. 轨迹采样点与起讫点

4. 从其他空间对象得到的代表点

2.0 点数据实战 1

2.0 专题研究 1：曼哈顿 Airbnb 每夜报价

专题研究 1：从房源点到空间关联

分析导航

研究课题

2.1 研究问题与 ESDA 框架

本例的研究问题

分析单位

研究区域与数据

2.2 研究区与三期数据

认识研究地点

发布日期、入住日期与每夜报价

数据预处理

2.3 价格分布、房源均价与筛选

网格化与聚合

2.4 六边形网格与中位报价

全局空间自相关

2.5 全局 Moran’s I

局部空间自相关

2.6 LISA：局部关联在哪里

为什么整体显著，多数格却未通过局部检验？

3.0 专题研究 1：回顾与反思

3.0 专题研究 1：回顾与反思

3.1 研究回顾：问题与回答

3.2 知识回顾：从点到区域

为什么把点汇总到网格？

为什么选六边形？方格可以吗？

3.3 研究设置的不足与改进

样本与报价口径

价格筛选与空间覆盖

空间聚合与尺度选择：当前模式是否稳定？

怎样组织一次敏感性分析？

对北部 LL 的解释还缺哪些资料？

3.4 地理学方法拓展与阅读

聚合之后，怎样定义“周围”？

数值变换与显示分级会怎样影响结果？

局部检验与区域解释各有什么边界？

拓展阅读：从概念到研究实施

4.0 专题研究 2

4.0 专题研究 2: 社交媒体点

专题研究 2: 上海六区微博记录空间分析

分析导航

课题与目标

4.1 研究课题与分析目标

研究区域与数据

4.2 研究区域与数据获取

数据预处理

4.3 数据预处理与初步展示

核密度估计

4.4 核密度估计：从点到连续表面

边界提取与可视化

4.5 高密度区提取与 AOI 边界

设施邻近性对比

4.6 到公交站、地铁站的距离比较

5.0 专题研究 2：回顾与反思

5.0 专题研究 2：回顾与反思

5.1 研究回顾：问题与回答

5.2 知识回顾：从记录到高密度区

原始点、KDE、AOI与设施距离分别回答什么？

为什么相同分位规则，记录覆盖比例却不同？

哪些设置只改图，哪些设置会改变结果？

观测单位、权重与显著性怎样一起读？

5.3 研究设置的不足与改进

先把“可能不稳定”变成可以复现的对照

样本、时间与行政边界分别限制什么？

AOI由同一批记录估计，不确定性还缺哪一层？

5.4 概念拓展：从平台记录到城市问题

5.4.1 社交媒体数据：一条记录怎样进入研究？

第五章：专题 1：点数据分析

章节作者：刘贇喆 | 华东师范大学社会发展学院