简介：「空间分析数据」压缩包面向GIS初学者与需要巩固ArcGIS实操的从业者，围绕空间位置、空间关系、空间统计及栅格与矢量分析等主题，提供可直接加载练习的数据素材，帮助解决缺少案例数据、难以走通完整分析流程的问题。包内共1555个文件，约98.08MB，包含adf、dat、nit等coverage数据，shp、shx、dbf等矢量图层，以及mxd地图文档、mdb地理数据库、img栅格影像和lyr图层等，可支撑数据浏览、制图、空间查询与工具箱分析。已有206人学习，适合按主题拆解练习并对照结果理解空间关系与统计含义。通过加载和操作这些数据，读者能熟悉ArcGIS多源数据的组织方式，完成缓冲区、叠加、网络与热点分析等典型任务，并进一步体会GIS在环境保护、城市规划、交通管理和灾害应对中的应用。
1. 一份门店表差点算错 380 公里：空间分析数据先要过的三道关
前些年接手一个门店选址的复盘，原始表里只有经度、纬度两列。有人写了句「距离小于 5 公里」，跑完发现某区域只剩 3 个候选，而业务方手上明明有 40 多家。把中间结果拉出来一看，两家实际相距不到 800 米的门店，算出来的「距离」是 0.0072——单位不是米，是度。
这就是空间分析数据最典型的样子：数据没错、代码不报错、结论全错。它和普通业务表的区别只有三点——每条记录带位置，位置必须挂在某个坐标系上，值不值钱取决于它和别的位置之间的拓扑关系（相交、包含、邻近、聚集）。所以落地时要过的也就三道关：位置怎么存、坐标系怎么统一、关系怎么算得快。做选址、物流调度、地理围栏、区域经营分析的人，基本都绕不开这三关。
2. 空间分析数据怎么存：GeoJSON、WKB 与 PostGIS 的取舍和建表命令
2.1 五种载体各自负责什么
同样的一个多边形，可以躺在五种不同的容器里，选错了后面每一步都别扭。判断标准其实就三个：要不要带坐标系、要不要建索引、要不要人能直接打开看。
| 载体 | 存储形态 | 自带坐标系 | 可建空间索引 | 典型用途 |
|---|---|---|---|---|
| GeoJSON | 文本 JSON | 规范要求 WGS84（EPSG:4326） | 否 | 接口传输、问题排查 |
| WKB / EWKB | 二进制 | EWKB 带 SRID | 否 | 数据库内部、序列化 |
| Shapefile | 多文件（.shp/.shx/.dbf/.prj） | 靠 .prj 旁挂 | 否 | 与老系统交换 |
| GeoPackage | 单文件 SQLite | 是 | 是 | 桌面端交付、离线包 |
| PostGIS geometry | 表字段 | 是 | GiST / SP-GiST | 查询、连接、聚合 |
GeoJSON 最容易被滥用。它可读性好，出问题时打开就能看到坐标，我排错时第一件事往往就是把可疑记录ST_AsGeoJSON出来瞄一眼。但它有两个硬伤：坐标是十进制文本，几百万个点体积能翻三四倍；没有索引，任何「附近」查询都是全量扫描。它的正确位置是「交换格式」和「调试格式」，不是「存储格式」。
Shapefile 的坑更隐蔽。字段名限 10 个字符、单文件 2GB 上限、属性表编码经常是 GBK 或 Latin-1 而非 UTF-8，中文属性读进来变乱码是常事。它现在只适合当历史数据的一次性入口，导完就该转走。
真正扛计算的是 PostGIS。geometry和geography两种字段类型，前者按平面算、后者按椭球面算，这个区别在第 3 章会展开——它是后面所有距离、面积、缓冲问题的总开关。
2.2 建表、灌数、加 GiST 索引的最小命令
先把表建对。字段类型里那个geometry(Point, 4326)的写法，等于给整张表上了一道保险：以后想往里塞 3857 的坐标，数据库直接拒绝，而不是默默算出一堆错结果。
-- 点表：门店 / 设备 / 订单地址 CREATE TABLE poi ( id bigserial PRIMARY KEY, name text, category text, lon double precision, lat double precision, geom geometry(Point, 4326) -- 类型 + SRID 双约束 ); -- 面表：商圈 / 行政区 / 配送范围 CREATE TABLE boundary ( id bigserial PRIMARY KEY, name text, geom geometry(MultiPolygon, 4326) ); -- 写入：先 ST_MakePoint，再 ST_SetSRID，顺序不能反 INSERT INTO poi (name, category, lon, lat, geom) VALUES ('样例点A', 'store', 116.40, 39.90, ST_SetSRID(ST_MakePoint(116.40, 39.90), 4326)); -- 空间索引：GiST 建在 geometry 上，实际索引的是外接矩形 CREATE INDEX idx_poi_geom ON poi USING GIST (geom); CREATE INDEX idx_boundary_geom ON boundary USING GIST (geom);批量导入走COPY比逐条 INSERT 快一个数量级，但要注意先落成文本列、再用表达式转 geom，避免中途某一行脏数据把整批打断：
# poi.csv 表头：name,category,lon,lat psql -d spatial -c "\copy poi_raw(name,category,lon,lat) FROM 'poi.csv' WITH (FORMAT csv, HEADER true)" psql -d spatial -c " INSERT INTO poi (name, category, lon, lat, geom) SELECT name, category, lon, lat, ST_SetSRID(ST_MakePoint(lon, lat), 4326) FROM poi_raw WHERE lon BETWEEN -180 AND 180 AND lat BETWEEN -90 AND 90; "参数上有三处值得留意。ST_MakePoint接收的是 (x, y)，在 4326 下 x 是经度、y 是纬度，写成 (lat, lon) 是很常见的翻车点。ST_SetSRID只是打标签，不做任何坐标换算，标错了数据库不会拦你，但后面一算就全错；真要换坐标系用ST_Transform。GiST 索引默认索引的是外接矩形（bbox），所以它只能做粗筛，精算还得靠后面的ST_Intersects这类函数，两段式是设计如此，不是浪费。
2.3 入库前后各做一次校验
空间分析数据最贵的地方在于错误不会报异常。所以我把校验做成固定动作，灌完立刻跑，别等出报表才回头查。
-- 1) 空值、越界、重复一次看全 SELECT count(*) AS total, count(*) FILTER (WHERE geom IS NULL) AS null_geom, count(*) FILTER (WHERE ST_X(geom) NOT BETWEEN -180 AND 180 OR ST_Y(geom) NOT BETWEEN -90 AND 90) AS out_of_range, count(DISTINCT geom::text) AS distinct_geom FROM poi; -- 2) 面数据有效性巡检，TopologyException 的根都在这里 SELECT id, name, ST_IsValidReason(geom) AS reason FROM boundary WHERE NOT ST_IsValid(geom); -- 3) 抽样反查：随机 10 个点落进哪个面 SELECT p.id, b.name FROM poi p LEFT JOIN boundary b ON ST_Intersects(p.geom, b.geom) ORDER BY random() LIMIT 10;第一步的out_of_range必须为 0，不为 0 说明经纬度顺序或者单位有问题；distinct_geom远小于total说明大量点堆在同一个默认坐标上，多半是解析失败被兜底成了 (0,0)。第二步只要出现自相交、环方向错误，ST_Contains在运行时就可能抛异常，或者更糟——静默返回空集。第三步是给人看的，落在海里或者落在几个不相邻的面里，肉眼一秒就能发现，比看聚合数字有效得多。
3. 坐标系选错，空间分析数据全白算：EPSG:4326 与 3857 的实测差距
3.1 度、米与形变：两个常用坐标系差在哪
EPSG:4326 是经纬度，单位是「度」，一度纬度大约 111.32 公里，一度经度是 111.32 × cos(纬度) 公里。这意味着同样 1 度的差，在高纬度和低纬度对应的地面距离完全不同——它就是不能直接拿来算距离和面积。
EPSG:3857（Web 墨卡托）把地球投影成正方形，单位是「米」，网页地图瓦片全靠它。但它为了让全球看起来是个正方形，在高纬度做了拉伸，尺度因子是 1/cos(纬度)：
| 纬度 | cos(纬度) | 3857 上量到 1 米，实际地面距离 |
|---|---|---|
| 0° | 1.000 | 1.00 m |
| 30° | 0.866 | 1.15 m |
| 40° | 0.766 | 1.31 m |
| 60° | 0.500 | 2.00 m |
中纬度城市做 5 公里半径筛选，用 3857 会系统性偏大 30% 左右。这个偏差不会报警，报表上就是「覆盖门店数比预期少了一截」，然后开始怀疑数据源。
我的分工是这样的：存储统一 4326，因为几乎所有数据源、GeoJSON 规范、GPS 原始值都是它；展示用 3857，对齐底图；计算走 geography 或本地投影坐标。三者不混用，混用必翻车。
3.2 三种算法实测：度、球面米、投影米
同一对点，三种写法给出三个数。下面这段可以直接拿去跑：
-- 写法一（错）：4326 上直接 ST_Distance，返回单位是度 SELECT ST_Distance(a.geom, b.geom) AS deg FROM poi a, poi b WHERE a.id = 1 AND b.id = 2; -- 写法二：转 geography，按 WGS84 椭球算，返回米 SELECT ST_Distance(a.geom::geography, b.geom::geography) AS m_geog FROM poi a, poi b WHERE a.id = 1 AND b.id = 2; -- 写法三：投影到本地 UTM 带再算，返回米 -- UTM 带号按经度推：zone = floor((lon + 180) / 6) + 1，北半球加 32600，南半球加 32700 WITH z AS ( SELECT floor((ST_X(a.geom) + 180) / 6)::int + 1 + 32600 AS srid FROM poi a WHERE a.id = 1 ) SELECT ST_Distance( ST_Transform(a.geom, (SELECT srid FROM z)), ST_Transform(b.geom, (SELECT srid FROM z)) ) AS m_proj FROM poi a, poi b WHERE a.id = 1 AND b.id = 2;逻辑说明：写法一返回的是角度差，数值通常在 0.001～0.1 量级，业务代码里乘个 111000 当米用，中纬度会直接错 20% 以上。写法二把点转成geography，计算时走椭球面公式，几十米到几千公里的距离精度足够，代价是比平面计算慢，且只能用于点线面之间的距离和包含，不支持部分投影类函数。写法三把数据搬到平面坐标系，之后所有距离、面积、缓冲都是纯平面运算，速度最快，适合频繁计算，但只在单个 UTM 带内准，跨带要分区。
参数上，::geography这一步转换本身有开销，如果一张表有几百万行、又要反复按距离筛选，我会额外加一列geog geography(Point,4326)并建 GiST 索引，而不是每次现转。UTM 带号那里注意经度 180° 边界：floor((180+180)/6)+1 = 61，得改成 60。
3.3 三条硬规则，写进团队规范
第一条：任何带「米」「公里」字样的条件，先确认列类型。是geometry就必须先ST_Transform或::geography，不允许直接比较。我把这条写进了代码评审清单，比事后查快得多。
第二条：面积统计不能用 4326。ST_Area(geom)在 4326 上返回的是「平方度」，这个数字没有任何物理意义，但它长得像个正常数值，所以特别容易被当成平方米贴进报表。要么ST_Area(geom::geography)，要么投影后再算。
第三条：跨带别硬投影。全国范围的数据投影到某一个 UTM 带，边缘区域形变能到百分之几。这种场景用geography，或者按纬度分带、分区计算再汇总。曾经有个区域覆盖率的指标，就是因为全国数据统一投到一个带里，边疆区域的面积整体偏大，指标看着漂亮了 6 个百分点，最后是靠换geography重算才发现的。
4. 空间连接从 40 分钟到 2 秒：索引、最近邻与网格聚合的写法
4.1 加了索引还是慢，先看执行计划里的 Seq Scan
最常见的疑惑是「GiST 索引建了啊，怎么还是慢」。九成原因是索引根本没被用上，看执行计划就知道：
EXPLAIN (ANALYZE, BUFFERS) SELECT count(*) FROM poi p JOIN boundary b ON ST_Intersects(p.geom, b.geom) WHERE b.name = '样例商圈';如果输出里出现Seq Scan on poi，且 rows 是整表行数，说明优化器放弃索引了。三个高频原因：字段被函数包住（ST_Transform(p.geom, 3857) && ...）、两边 SRID 不一致导致隐式转换、统计信息过期（大批量导入后没跑ANALYZE）。
前两个的解法是让索引和条件长得一样。需要频繁在 3857 上做查询，就物化一列：
-- 生成列 + 表达式索引：条件和索引完全对齐，优化器才认 ALTER TABLE poi ADD COLUMN geom_3857 geometry(Point, 3857) GENERATED ALWAYS AS (ST_Transform(geom, 3857)) STORED; CREATE INDEX idx_poi_geom_3857 ON poi USING GIST (geom_3857); ANALYZE poi;STORED表示写入时算好落盘，读的时候零成本；代价是占空间、写入略慢。生成列要求表达式IMMUTABLE，ST_Transform(geom, integer)满足这个条件。版本较老、不支持生成列的环境，退回触发器或者写入时直接算好，效果一样。
4.2 三种典型查询：点面判定、最近邻、半径筛选
点落在哪个面里，用ST_Intersects，它内部先比 bbox 再用索引，最后才做精确判定：
SELECT p.id, p.name, b.name AS area_name FROM poi p JOIN boundary b ON ST_Intersects(p.geom, b.geom) WHERE p.category = 'store';最近的 5 个点，用 KNN 的<->算子，这是少数几个能直接吃 GiST 索引的排序：
SELECT id, name FROM poi WHERE category = 'store' ORDER BY geom <-> ST_SetSRID(ST_MakePoint(116.40, 39.90), 4326) LIMIT 5;注意<->比较的是外接矩形中心到查询点的距离，不是真实球面距离。点数密集时排序会有细微偏差，稳妥做法是先取LIMIT 20做粗排，外面套一层用ST_Distance精排取前 5。
3 公里内有多少点，用ST_DWithin，它的关键优势是能走索引，而ST_Distance(...) < 3000不能：
SELECT count(*) FROM poi WHERE ST_DWithin( geom::geography, ST_SetSRID(ST_MakePoint(116.40, 39.90), 4326)::geography, 3000 -- 单位：米，因为左边是 geography );ST_DWithin的第三个参数单位跟随字段类型：geography是米，geometry是坐标系单位。这是又一个「不报错但算错」的地方，写的时候最好在注释里把单位标死。
4.3 百万级点集别硬做连接：网格与 H3 分桶
两张百万行的表做空间连接，即使有索引也快不到哪去，因为每个点都要真去探一次索引。这种规模的聚合，正确姿势是先降维——把点按空间分桶，桶内计数，再拿桶去做分析。
规则网格最简单：
-- 按 0.01 度（约 1.1 公里）分桶 SELECT floor(ST_X(geom) / 0.01)::int AS gx, floor(ST_Y(geom) / 0.01)::int AS gy, count(*) AS cnt FROM poi GROUP BY 1, 2;floor而不是round，是为了保证每个点只落进一个桶——闭区间取整会让边界点被相邻两个桶各算一次，总计数对不上。这个问题排查起来很费劲，因为每个桶单独看都「差不多对」。
规则网格的毛病是边界效应：一个点离邻格中心比离本格中心还近，八个邻居到中心的距离也不一样（正交 1 格、对角 1.414 格）。做密度聚类时这会制造假的「热点边缘」。所以点位聚合我一般用 H3 六边形：
| 分辨率 | 平均单元面积 | 适合的聚合粒度 |
|---|---|---|
| 7 | 约 5.16 km² | 城市级热力、大区对比 |
| 8 | 约 0.74 km² | 商圈级、配送站覆盖 |
| 9 | 约 0.11 km² | 街区级、门店辐射 |
import h3 import pandas as pd df = pd.read_parquet("orders.parquet") # 列：order_id, lon, lat # h3 v4 的函数名是 latlng_to_cell，v3 里叫 geo_to_h3，换版本会直接报错 df["cell"] = [h3.latlng_to_cell(lat, lon, 8) for lat, lon in zip(df.lat, df.lon)] agg = (df.groupby("cell") .agg(orders=("order_id", "count")) .reset_index()) # 取六边形顶点，方便落库或画图 agg["boundary"] = agg["cell"].map(lambda c: h3.cell_to_boundary(c))逻辑说明：latlng_to_cell把经纬度映射到一个 64 位整数 ID，同格内的点 ID 相同，聚合退化成普通groupby，没有几何运算，千万行也就几十秒。参数res是分辨率，取值 0～15，越大格子越小；选哪个取决于业务粒度，不是越细越好——格子太细每格样本量不足，后续做聚集性检验时统计功效不够，这个在第 6 章会具体说。
六边形相比方格的核心优势是邻居一致性：六个邻居到中心的距离完全相同，没有对角邻居的歧义，做邻接权重矩阵时天然对应 k=6。这套 ID 还能直接当分区键用，把「同一个格子附近的数据」物理上放到一起，后续连接和缓存都更友好。
5. 空间分析数据排查清单：5 个反复翻车的现场
5.1 距离算出来是 0.007，不是 7 公里
现象：筛选「5 公里内」的门店，结果要么只有几家，要么几十万行全中。
原因：字段是geometry(...,4326)，ST_Distance返回的是度。0.007 度约等于 780 米，代码里乘 111000 在低纬度勉强能用，中高纬度直接错两三成。
解决：先确认字段类型，再决定写法。要么::geography走米，要么ST_Transform到本地投影。改完在两个已知距离的点上回归一次，把结果写进单测——这类错误如果不在代码里固定下来，下次换个人写还会再来一遍。
5.2 索引建了，执行计划还是全表扫
现象：表上明明有 GiST 索引，带空间条件的查询要跑几十秒，EXPLAIN里是Seq Scan。
原因：条件里的字段被函数包住了，或者两边 SRID 不一致触发隐式转换，优化器认为索引用不上。另一种情况是刚批量导入几百万行、统计信息还是旧的，优化器估错了行数。
解决：先ANALYZE一次；还在扫就检查条件表达式，把ST_Transform(geom, 3857) && ...改成物化列geom_3857 && ...并给它建索引。判断标准很简单：EXPLAIN里出现Index Scan using idx_xxx才算过。
5.3 多边形一相交就抛 TopologyException
现象：两个面做叠置或者包含判断，跑到某几行报GEOSTopologyException: Self-intersection，整套任务中断。
原因：面数据是多源拼接来的，存在自相交、重复点、环方向错误。这些几何从数据源看「画得挺好」，但拓扑上不合法。
解决：入库时就用ST_IsValidReason巡检，把不合法的挑出来单独处理；确实修不了的用ST_MakeValid(geom)兜底，修完再ST_Multi统一成 MultiPolygon，避免类型不匹配。注意ST_MakeValid有可能把多边形拆成 GeometryCollection，落库前要过滤掉非面要素。
5.4 点全跑到海里去了
现象：地图上撒点，一大片落在海上或者地图边缘，但原始 CSV 里的经纬度数值看着完全正常。
原因：经纬序反了。GeoJSON 规范是[lon, lat]，很多库的构造参数是(lat, lon)，ST_MakePoint(x, y)又是先 x 后 y。三套约定混在一起，写错一次很难从数值上看出来。
解决：入库时加范围断言——国境内的经度大概在 73～135、纬度在 3～54（按实际覆盖范围收紧），越界直接落拒绝表，别进主表。开发阶段随机抽 10 条反查行政区，肉眼看一遍，比读代码快。
5.5 网格聚合总计数和明细对不上
现象：按网格聚合后的订单总数，比原始表少了一部分或者多了一部分，差值不大，两三个百分点。
原因：分桶时用了闭区间取整（round或者ceil），边界上的点被两个桶各算一次，或者一个都没算。跨 180° 经线的数据如果直接按floor(lon/x)分桶，还会出现负索引错位。
解决：统一用floor，保证半开区间[a, b)，点只属于一个桶。跨经线的数据先把经度归一化到[0, 360)再分桶，或者干脆用 H3 这类基于球面的编码，避开平面坐标的边界问题。聚合完一定做一次SUM(桶内计数) = 原始行数的校验，这个断言几乎零成本，却能挡住一整类错误。
6. 结论是聚集还是巧合：用全局莫兰指数给空间分析数据做一次体检
地图上看到一片深色，我们会本能地说「这里聚集」。但深色可能只是样本少，也可能只是整体分布的自然波动。要让「聚集」这个结论站得住，得给它配一个统计量——全局莫兰指数（Global Moran's I）就是最常用的那个。
它的直觉很简单：把每个格子的值和邻居的值相乘求和，再按方差归一化。值相近且相邻，乘积为正，指数就高；高低交错，乘积为负，指数就低。取值大致落在 -1 到 1，期望值约等于 -1/(n-1)——注意不是 0，格子数少的时候这个差别不可忽略。真正要看的不是 I 的绝对值，而是置换检验出来的伪 p 值：把观测值随机打乱 999 次，看真实 I 在随机分布里排第几。
import geopandas as gpd from libpysal.weights import KNN from esda.moran import Moran gdf = gpd.read_file("hex_agg.gpkg") # 已按 H3 res 8 聚合 y = gdf["order_rate"].values # 六边形一阶邻居正好 6 个，k 取 6 与几何结构对应 w = KNN.from_dataframe(gdf, k=6) w.transform = "r" # 行标准化，避免边缘格子邻居少而权重偏低 mi = Moran(y, w, permutations=999) print(f"I={mi.I:.4f} E[I]={mi.EI:.4f} z={mi.z_sim:.2f} p_sim={mi.p_sim:.4f}")三个参数值得改着试。k的取值依赖几何结构：六边形取 6，方格用 Queen 邻接一般取 8，纯距离带宽则在密度不均时会出问题（中心密、边缘稀，同一带宽在两地含义不同）。transform="r"做行标准化，让每个格子的邻居权重和为 1，否则边缘格子邻居少、贡献天然偏低。permutations至少 999，要报更小的 p 值就加到 9999。
这套东西最大的价值是当分辨率的敏感性检查。同一个指标在 H3 res 7、8、9 上各跑一遍，如果三个尺度上 I 都显著为正，结论是稳的；如果只在某一个分辨率上显著，那多半是尺度选出来的假象——格子大到把两个功能区合成一个，或者小到每格样本只有个位数，都会凭空造出聚集。我吃过这个亏：一份区域订单密度分析，res 9 上 I 显著，换到 res 8 就掉到不显著，回头一查是细格子里大量零值把方差抬高了。
我现在的习惯是，任何一个带距离、面积、「附近」字样的结论，交付前先问三句：单位是什么、坐标系是哪个、索引走没走。三问答不上来，图再好看也不往下签字。希望帮到你。
本文还有配套的精品资源，点击获取
