简介「空间分析数据」压缩包面向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规范要求 WGS84EPSG: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:3857Web 墨卡托把地球投影成正方形单位是「米」网页地图瓦片全靠它。但它为了让全球看起来是个正方形在高纬度做了拉伸尺度因子是 1/cos(纬度)纬度cos(纬度)3857 上量到 1 米实际地面距离0°1.0001.00 m30°0.8661.15 m40°0.7661.31 m60°0.5002.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.0010.1 量级业务代码里乘个 111000 当米用中纬度会直接错 20% 以上。写法二把点转成geography计算时走椭球面公式几十米到几千公里的距离精度足够代价是比平面计算慢且只能用于点线面之间的距离和包含不支持部分投影类函数。写法三把数据搬到平面坐标系之后所有距离、面积、缓冲都是纯平面运算速度最快适合频繁计算但只在单个 UTM 带内准跨带要分区。参数上::geography这一步转换本身有开销如果一张表有几百万行、又要反复按距离筛选我会额外加一列geog geography(Point,4326)并建 GiST 索引而不是每次现转。UTM 带号那里注意经度 180° 边界floor((180180)/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表示写入时算好落盘读的时候零成本代价是占空间、写入略慢。生成列要求表达式IMMUTABLEST_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_cellv3 里叫 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是分辨率取值 015越大格子越小选哪个取决于业务粒度不是越细越好——格子太细每格样本量不足后续做聚集性检验时统计功效不够这个在第 6 章会具体说。六边形相比方格的核心优势是邻居一致性六个邻居到中心的距离完全相同没有对角邻居的歧义做邻接权重矩阵时天然对应 k6。这套 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。三套约定混在一起写错一次很难从数值上看出来。解决入库时加范围断言——国境内的经度大概在 73135、纬度在 354按实际覆盖范围收紧越界直接落拒绝表别进主表。开发阶段随机抽 10 条反查行政区肉眼看一遍比读代码快。5.5 网格聚合总计数和明细对不上现象按网格聚合后的订单总数比原始表少了一部分或者多了一部分差值不大两三个百分点。原因分桶时用了闭区间取整round或者ceil边界上的点被两个桶各算一次或者一个都没算。跨 180° 经线的数据如果直接按floor(lon/x)分桶还会出现负索引错位。解决统一用floor保证半开区间[a, b)点只属于一个桶。跨经线的数据先把经度归一化到[0, 360)再分桶或者干脆用 H3 这类基于球面的编码避开平面坐标的边界问题。聚合完一定做一次SUM(桶内计数) 原始行数的校验这个断言几乎零成本却能挡住一整类错误。6. 结论是聚集还是巧合用全局莫兰指数给空间分析数据做一次体检地图上看到一片深色我们会本能地说「这里聚集」。但深色可能只是样本少也可能只是整体分布的自然波动。要让「聚集」这个结论站得住得给它配一个统计量——全局莫兰指数Global Morans 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, k6) w.transform r # 行标准化避免边缘格子邻居少而权重偏低 mi Moran(y, w, permutations999) print(fI{mi.I:.4f} E[I]{mi.EI:.4f} z{mi.z_sim:.2f} p_sim{mi.p_sim:.4f})三个参数值得改着试。k的取值依赖几何结构六边形取 6方格用 Queen 邻接一般取 8纯距离带宽则在密度不均时会出问题中心密、边缘稀同一带宽在两地含义不同。transformr做行标准化让每个格子的邻居权重和为 1否则边缘格子邻居少、贡献天然偏低。permutations至少 999要报更小的 p 值就加到 9999。这套东西最大的价值是当分辨率的敏感性检查。同一个指标在 H3 res 7、8、9 上各跑一遍如果三个尺度上 I 都显著为正结论是稳的如果只在某一个分辨率上显著那多半是尺度选出来的假象——格子大到把两个功能区合成一个或者小到每格样本只有个位数都会凭空造出聚集。我吃过这个亏一份区域订单密度分析res 9 上 I 显著换到 res 8 就掉到不显著回头一查是细格子里大量零值把方差抬高了。我现在的习惯是任何一个带距离、面积、「附近」字样的结论交付前先问三句单位是什么、坐标系是哪个、索引走没走。三问答不上来图再好看也不往下签字。希望帮到你。本文还有配套的精品资源点击获取