简介《空间分析数据》面向GIS初学者与ArcGIS实践者围绕空间位置、空间关系、空间统计及栅格与矢量分析等主题帮助读者在真实数据中理解缓冲区、网络分析、聚类和热点分析等方法的实现过程。压缩包为rar格式共1555个文件约98.08MB包含shp、dbf、prj等矢量数据adf、dat等栅格文件以及mxd地图文档、mdb地理数据库和aux、xml等辅助文件覆盖数据加载、地图制图到空间分析工具运行的完整链路。已有206人学习适合课程实验、项目演练与技能查漏补缺读者可借此熟悉矢量与栅格数据的组织方式练习空间查询、制图表达和结果解读并探索城市规划、交通管理、灾害应对等应用场景积累可复用的空间思维与操作经验。资源目录按分析环节组织便于循序渐进地完成实验与复盘。1. 空间分析数据不是“带经纬度的表”它多算了哪几层关系一份带经纬度的点位表画成散点图往往挺像样一旦要算密度、最近邻、可达性或热点结果就开始飘。空间分析数据真正要落地的不只是两列坐标而是几何对象、坐标参考系、拓扑关系、尺度假设和空间权重这几层叠在一起的数据形态。点位、线、面在数据库里怎么存距离按米还是按度算邻居按接边还是按半径算格网切多大都会改变结论。它适合做选址评估、资源覆盖、事件聚集和区域对比的团队尤其适合已经有一堆业务表、却总在“图好看但数不对”上反复踩坑的人。先把这层关系讲清楚后面的代码和参数才有意义。2. 从点位表到可分析的空间分析数据坐标系、几何与索引的最小闭环拿到一份点位表最容易犯的错是直接对经度、纬度做距离计算。经度一度在不同纬度对应的地面距离并不相同纬度一度也随椭球模型变化度制坐标上的 0.01 并不等于固定米数。空间分析数据的第一道门槛不是算法而是坐标参考系存储层用经纬度没问题计算层必须明确“我现在算的是度还是米”。这一步没锁死后面缓冲区、最近邻、密度格网都会变成橡皮筋。2.1 先锁坐标系EPSG:4326 只适合存算距离前必须换投影我一般把原始点位按 EPSG:4326 存因为交换和入库通用做面积、距离、缓冲区时再转到等积或局部投影。EPSG:3857 常被拿来做底图但它会放大高纬度面积不适合算密度和面积。全球尺度可以用 EPSG:6933 这类等积投影局部尺度则按研究区选择局部等积投影。下面这段代码先把坐标系显式化再转到米制空间。import geopandas as gpd from pyproj import CRS gdf gpd.read_file(data/site_points.gpkg) if gdf.crs is None: gdf gdf.set_crs(4326) gdf_ll gdf.to_crs(4326) # 交换、入库用经纬度 gdf_m gdf_ll.to_crs(6933) # 等积投影用于面积和密度 print(CRS.from_epsg(6933).axis_info[0].unit_name) # 输出应为 metre逻辑说明set_crs只声明坐标系不做坐标变换to_crs才会重算坐标。参数上4326 的单位是度6933 的单位是米如果数据本身是局部投影不要硬转 6933应按研究区选择局部等积投影。判断标准很简单代码里出现buffer(1000)时当前 CRS 的单位必须是米否则 1000 就是 1000 度结果必然翻车。2.2 几何有效性点去重、线自交、面环方向与修复空间分析数据里的几何错误不一定报错但会让空间连接悄悄丢记录。面自交、环方向异常、空几何、重复点位都是常见来源。点数据虽然很少出现无效几何但重复点位会把计数和密度抬高线和面则要先做有效性检查。我通常先按业务主键去重再用make_valid修复几何最后把空几何剔除。import geopandas as gpd from shapely.validation import make_valid gdf gpd.read_file(data/site_points.gpkg).to_crs(6933) gdf gdf.drop_duplicates(subset[site_id]).copy() gdf[geometry] gdf.geometry.apply( lambda g: make_valid(g) if g is not None and not g.is_valid else g ) gdf gdf[gdf.geometry.notna() ~gdf.geometry.is_empty] gdf.to_file(data/site_points_clean.gpkg, driverGPKG)逻辑说明site_id是业务唯一键重复点先按它去重避免同一设施被重复统计。make_valid针对无效面最有用点几何一般不受影响preserve_topology不是这里必须的参数因为make_valid的目标是得到合法几何。参数driverGPKG表示输出为 GeoPackage适合保留 CRS 和属性。修复后不要直接覆盖原表保留一份 clean 版本后面排查有后悔药。2.3 空间连接sjoin 的 predicate、半径和候选集空间连接是把“点属于哪个面”或“点附近有哪些点”落成表的关键步骤。predicate决定关系within是点落在面内intersects是相交contains是包含。howleft会保留未匹配的点适合检查覆盖缺口howinner只留匹配上的记录适合做统计但容易掩盖边界问题。半径连接要盯住 CRS 单位米制投影下max_distance1000才是 1 公里。import geopandas as gpd points gpd.read_file(data/site_points_clean.gpkg).to_crs(6933) zones gpd.read_file(data/zone_boundary.gpkg).to_crs(6933) inside gpd.sjoin( points, zones[[zone_id, geometry]], howleft, predicatewithin ) near gpd.sjoin_nearest( points, zones[[zone_id, geometry]], howleft, max_distance1000, distance_coldist_m )逻辑说明inside里未匹配到zone_id的点要么在边界外要么面几何有问题先查再删。sjoin_nearest的max_distance使用当前 CRS 的单位6933 下是米如果还在 4326 下1000 会被当成度。distance_col保留实际距离方便后续按 300 米、500 米、1000 米做敏感性对比。空间连接前先对两个图层做to_crs对齐别指望库自动猜。2.4 落库geometry 与 geography 的边界单机脚本跑通后空间分析数据通常要落到空间数据库供多人查询和增量更新。这里最常纠结的是geometry和geographygeometry快适合投影坐标geography按椭球面算米制距离适合直接存经纬度又想做半径查询。我的做法是主表存 4326 的geometry需要米制半径时再转geography或者额外存一列投影几何做重计算。CREATE TABLE fact_site_point ( site_id text PRIMARY KEY, zone_id text, metric_value numeric, updated_at timestamptz, geom geometry(Point, 4326) ); CREATE INDEX idx_fact_site_point_geom ON fact_site_point USING gist (geom); -- 1 km 半径邻居geography 按米计算 SELECT a.site_id, b.site_id AS neighbor_id FROM fact_site_point a JOIN fact_site_point b ON ST_DWithin(a.geom::geography, b.geom::geography, 1000) WHERE a.site_id b.site_id;逻辑说明geometry(Point, 4326)把类型和 SRID 一起约束住避免不同坐标系混入同一列。GIST 索引是空间查询的入口没有它ST_DWithin容易退化成全表扫描。a.site_id b.site_id用来去掉重复点对和自配对。参数 1000 的单位由geography决定是米如果数据量大::geography的转换成本要评估可以增加投影几何列并在该列上建索引。3. 空间分析数据怎么表达“邻居”空间权重、格网聚合与热点检验空间分析数据和普通表格最大的差别是观测值之间不再独立。一个点的高值可能受邻近点影响一个面的属性也可能受接边单元影响。空间权重就是把“谁和谁相邻”写成矩阵面可以用 Queen 或 Rook 邻接点可以用 KNN 或距离阈值。权重选错Morans I 和热点分析就会给出看似显著、实际不可解释的结果。这里要先问关系假设再谈算法。3.1 Queen、Rook、KNN 与距离阈值先问关系假设面数据的 Queen 邻接把共点也算邻居Rook 只算共边点数据更常用 KNN 或距离带。Queen 适合边界接触也算邻接的业务单元Rook 适合只认共享边界的场景。KNN 的k不是越大越好k 太小容易孤岛k 太大把远处无关单元拉进来。距离阈值要有业务含义比如 1 公里覆盖圈而不是为了让结果显著随便调。权重还要决定是否行标准化行标准化后每个单元的邻居权重和为 1适合比较不同邻居数量的单元。3.2 用 PySAL 建权重并检查孤岛构建权重后第一件事不是跑统计而是看孤岛。孤岛是没有邻居的单元它们会让局部统计出现空值也会让全局指标偏移。下面的代码对面数据建 Queen 权重对点数据建 KNN 权重并做一次 Morans I 检查。import geopandas as gpd from libpysal.weights import Queen, KNN from esda.moran import Moran zones gpd.read_file(data/zone_boundary.gpkg).to_crs(6933) queen Queen.from_dataframe(zones, use_indexTrue) queen.transform r print(孤岛数量:, len(queen.islands)) points gpd.read_file(data/site_points_clean.gpkg).to_crs(6933) knn KNN.from_dataframe(points, k6) knn.transform r mi Moran(zones[density].values, queen, permutations999) print(mi.I, mi.p_sim)逻辑说明use_indexTrue让权重矩阵使用 GeoDataFrame 的索引后续对齐属性更稳。transformr是行标准化适合邻居数量不一致的单元。permutations999做随机置换检验p_sim是模拟显著性样本很小时不要只盯 p 值。k6只是起点实际要按业务半径和点密度调调完记录参数别让下一次跑批变成黑匣子。3.3 格网聚合H3、GeoHash 与业务单元不是一回事业务单元边界往往不规则面积差异大直接比较计数会受面积影响。格网聚合把空间切成等面积或近似等面积的单元适合做密度、热力和跨区域对比。H3 适合六边形层级聚合GeoHash 适合字符串索引和快速分桶规则方格最容易用代码实现。代价是格网边界会切断真实业务边界尺度一变结论可能变这就是空间分析里的尺度效应。import numpy as np import geopandas as gpd from shapely.geometry import Point points gpd.read_file(data/site_points_clean.gpkg).to_crs(6933) grid_size 500 # 米 points[gx] np.floor(points.geometry.x / grid_size).astype(int64) points[gy] np.floor(points.geometry.y / grid_size).astype(int64) cells ( points.groupby([gx, gy]) .agg(point_count(site_id, size), metric_mean(metric_value, mean)) .reset_index() ) cells[geometry] [ Point(x * grid_size grid_size / 2, y * grid_size grid_size / 2) for x, y in zip(cells[gx], cells[gy]) ] cells gpd.GeoDataFrame(cells, geometrygeometry, crspoints.crs)逻辑说明grid_size500表示 500 米方格先除以格宽再取整得到行列号。聚合时size统计点数量mean统计指标均值如果业务上更关心总和把mean换成sum。生成格网中心点时乘回格宽并加半格保证点落在格子中心。H3 和 GeoHash 的实现方式不同但原则一样先定层级或精度再做多尺度对比不要只交一个尺度。3.4 热点分析Getis-Ord Gi* 前先确认权重和样本量热点分析不是给每个点打分而是看局部高值或低值是否在空间上聚集。Getis-Ord Gi* 的starTrue表示包含自身适合找高值簇starFalse排除自身解释会变。权重仍用 KNN 或距离带样本量太小的格网不要强行算否则 p 值会被少量样本主导。多重比较也要注意几百个格网同时检验假阳性会累积报告里应写清显著性阈值和校正方式。from libpysal.weights import KNN from esda.getisord import G_Local cells cells[cells[point_count] 3].copy() w KNN.from_dataframe(cells, k8) w.transform r gi G_Local( cells[point_count].astype(float).values, w, starTrue, permutations999 ) cells[gi_z] gi.Zs cells[gi_p] gi.p_sim逻辑说明先过滤point_count 3避免稀疏格网制造伪热点。k8是邻域数量格网越密可以适当增大但要在报告里固定。starTrue让自身参与计算permutations999做随机检验。gi_z是局部 z 分数gi_p是模拟 p 值解释时不要只看红色高值还要看它是否通过显著性阈值以及换一个格网尺寸后是否仍成立。4. 空间分析数据翻车排查五个高频现象和对应处置空间分析数据的坑很少以报错形式出现更多是结果看起来合理、实际已经偏了。下面五条是血泪经验里最常见的一类空间连接丢记录、缓冲区单位错、查询慢、尺度反转、重复点污染。每条都按“现象、原因、解决”拆开排查时先复现现象再定位原因最后用小样本验证修复。4.1 现象空间连接后大量空值原因几何无效或 CRS 不一致解决先对齐再修复现象是sjoin之后zone_id空值比例异常高但点位明明落在面内。原因通常有两个两个图层 CRS 不一致或者面几何自交、环方向异常导致within判断失败。解决顺序是先to_crs对齐再用make_valid修复面几何最后用howleft保留未匹配记录逐条查。不要一上来就删空值空值本身就是排查线索。修复后重新跑连接并对比修复前后的匹配率匹配率突增说明原几何确实有问题。4.2 现象1000 米缓冲区像几百公里原因4326 下把度当米解决投影后再 buffer现象是buffer(1000)出来的圆巨大或者半径连接把整片区域都连上。原因很直接当前 CRS 是 4326单位是度1000 度当然离谱。解决是先把数据转到米制投影再执行buffer或sjoin_nearest。如果必须保留经纬度可以用geography做米制半径查询但要接受转换成本。排查时打印gdf.crs.axis_info[0].unit_name看到 degree 就不要做米制操作。这个坑很常见留下一条检查习惯比记住某个参数更管用。4.3 现象空间查询几十万点就慢原因索引没用上或函数包住几何列解决先用 bbox 走 GIST再精确过滤现象是空间数据库查询随数据量线性变慢EXPLAIN ANALYZE里出现全表扫描。原因常见于WHERE ST_DWithin(...)直接包住几何列或者没有 GIST 索引或者先用函数计算再过滤。解决是先用geom ST_MakeEnvelope(...)做边界框粗筛让 GIST 索引生效再用ST_DWithin做精确过滤。下面这条查询把粗筛和精确过滤分开便于看执行计划。EXPLAIN ANALYZE SELECT site_id FROM fact_site_point WHERE geom ST_MakeEnvelope(:xmin, :ymin, :xmax, :ymax, 4326) AND ST_DWithin( geom::geography, ST_SetSRID(ST_Point(:lon, :lat), 4326)::geography, 1000 );逻辑说明是边界框相交能走 GIST 索引ST_DWithin再做米制精确判断。参数:xmin、:ymin、:xmax、:ymax来自目标点的外扩范围1000是米。如果执行计划仍不走索引检查列类型是否为geometry、索引是否建在同一个列上、统计信息是否过期。别用ST_Transform包住列再比较那会让索引失效。4.4 现象换一个格网尺寸结论反转原因尺度效应和边界切割解决多尺度敏感性对比现象是 500 米格网显示高值聚集1 公里格网却变成随机分布。原因不一定是代码错而是尺度效应格网变大后高值和低值被平均边界切割也改变了归属。解决是至少跑三档尺度比如 250 米、500 米、1000 米比较热点位置和 Morans I 的稳定性。报告里不要只放一张最显著的结果图要说明尺度选择依据。若业务单元本身有含义优先用业务单元再用格网做交叉验证不要为了图好看换尺度。4.5 现象密度图局部异常高原因重复点或多记录合并解决唯一键去重并明确权重现象是某个格网密度突然高出周围几倍回到原始表发现有同一设施的多条记录或者一条业务记录对应多个服务能力。原因是没有按分析粒度去重也没有定义权重。解决是先确定分析单元一个设施一条还是一个服务能力一条前者按site_id去重后者保留记录但把metric_value作为权重。密度计算前检查主键重复率重复率高就先和业务确认不要直接用drop_duplicates把有效多值记录删掉。5. 空间分析数据从单机脚本到可重复跑批分区、简化与增量更新单机脚本能跑通不等于能稳定交付。数据量上来后全量读入内存、每次重算所有空间连接、几何精度过高都会让跑批时间失控。更稳的做法是把空间分析数据拆成可重复的层分区裁剪控制内存几何简化控制传输量物化中间表控制重复计算增量更新控制每次跑批的范围。这里的目标不是追求最快而是让每次结果可解释、可回滚。5.1 分区裁剪先按网格切任务不让一次全量进内存全量read_file在数据量大时很容易把内存打满尤其是面和面相交。我一般按固定大小的网格切任务每个任务只处理一个瓦片输出带瓦片编号的中间文件。瓦片大小要大于最大几何的跨度否则一个面会被切碎但也不能太大否则失去分区意义。下面的代码用 50 公里瓦片切点数据面数据同理只是裁剪时要保留与瓦片相交的面。import geopandas as gpd from shapely.geometry import box all_points gpd.read_file(data/site_points_clean.gpkg).to_crs(6933) xmin, ymin, xmax, ymax all_points.total_bounds tile 50000 # 50 km for ix in range(int(xmin // tile), int(xmax // tile) 1): for iy in range(int(ymin // tile), int(ymax // tile) 1): clip box(ix * tile, iy * tile, (ix 1) * tile, (iy 1) * tile) part all_points[all_points.intersects(clip)] if not part.empty: part.to_file(fdata/tiles/points_{ix}_{iy}.gpkg, driverGPKG)逻辑说明total_bounds给出整体范围tile50000表示 50 公里。intersects保留与瓦片相交的点边界点可能同时落入两个瓦片后续统计要按主键去重。输出文件名带ix、iy方便重跑单个瓦片。参数不要照搬点密集就缩小瓦片面很大就放大瓦片原则是单任务内存可控且不频繁切碎几何。5.2 几何简化拓扑简化与容差选择的边界几何简化能减少传输和渲染成本但简化后的几何不能直接用于面积统计。simplify(tolerance10, preserve_topologyTrue)表示容差 10 米并保持拓扑适合展示和粗筛面积、长度仍用原始几何算。容差太大会让相邻面出现缝隙或重叠空间连接结果随之改变。我的习惯是保留两列原始几何用于计算简化几何用于出图和前端加载并在元数据里写清容差。zones gpd.read_file(data/zone_boundary.gpkg).to_crs(6933) zones[geom_simplified] zones.geometry.simplify( tolerance10, preserve_topologyTrue ) zones[area_m2] zones.geometry.area # 面积仍用原始几何逻辑说明tolerance10的单位是米来自 6933如果 CRS 是度10 就完全不是一个量级。preserve_topologyTrue能减少自交和缝隙但不能保证简化后完全满足业务拓扑关键图层仍要做有效性检查。area_m2用原始几何计算避免简化误差进入指标。简化几何单独存列或单独出表不要覆盖原始几何。5.3 物化中间表把空间连接结果固化成分析宽表每次查询都实时做空间连接数据量大时会把数据库拖慢。更稳的做法是把点面关系物化成分析宽表后续聚合直接查宽表。物化表要保留业务主键、空间单元、指标值、更新时间和几何方便回溯。下面 SQL 把点位落到业务单元并建普通索引供按单元聚合。CREATE TABLE analysis_point_zone AS SELECT p.site_id, z.zone_id, p.metric_value, p.updated_at, p.geom FROM fact_site_point p LEFT JOIN dim_zone z ON ST_Within(p.geom, z.geom); CREATE INDEX idx_analysis_point_zone_zone ON analysis_point_zone (zone_id);逻辑说明LEFT JOIN保留未落入任何业务单元的点后续可以单独查边界问题。ST_Within要求点完全落在面内边界点可能不匹配如果业务上边界算归属改用ST_Intersects并加去重规则。物化表不是一劳永逸源表更新后要同步刷新刷新策略见下一节。索引建在zone_id上适合按业务单元汇总。5.4 增量更新按更新时间戳和空间分区双条件全量重建物化表简单但成本高增量更新要解决两个问题哪些记录变了以及删除和插入如何保持幂等。我一般用updated_at :since拉取变更再按主键先删后插如果表按空间分区还可以加瓦片编号缩小范围。下面的 SQL 是事务内的先删后插参数:since由上一次成功跑批的时间决定。BEGIN; DELETE FROM analysis_point_zone WHERE updated_at :since; INSERT INTO analysis_point_zone (site_id, zone_id, metric_value, updated_at, geom) SELECT p.site_id, z.zone_id, p.metric_value, p.updated_at, p.geom FROM fact_site_point p LEFT JOIN dim_zone z ON ST_Within(p.geom, z.geom) WHERE p.updated_at :since; COMMIT;逻辑说明先删后插保证同一主键不会重复事务保证删除和插入要么都成功要么都回滚。:since要使用数据库时间而不是脚本本地时间避免时区偏差。如果源表存在延迟到达的记录时间戳窗口要留重叠比如回退 10 分钟再跑一次。增量跑完后对账源表变更条数、宽表删除条数、宽表插入条数应当能对上对不上就查时间窗口和主键重复。6. 空间分析数据交付前的三项校验结果能复现比图画得漂亮更重要交付前我一般做三项校验不靠肉眼盯图。第一项是几何校验空几何、无效几何、CRS 是否统一第二项是属性校验主键重复率、空间连接空值率、指标极值是否越界第三项是统计校验固定样本跑 Morans I 或热点结果和上一次跑批对比。三项都过再谈出图。校验项最小动作通过标准几何检查 CRS、空几何、is_valid无空几何CRS 一致无效几何已修复属性统计主键重复率和连接空值率重复率为 0空值率有解释统计固定样本重跑 Morans I 和热点指标在阈值内热点位置稳定固定样本不需要很大抽几百条覆盖边界、中心、稀疏区即可。把样本主键、CRS、格网尺寸、权重参数写进跑批配置下一次换人跑也能复现。我现在的习惯是任何空间分析数据交付前先问一句“换一台机器、换一个格网这个结论还成立吗”答不上来就先补校验。这个方向值得投入但别把图画得漂亮当成结果可靠。希望帮到你。本文还有配套的精品资源点击获取
