企业官网建设流程全解析

简介：《空间分析数据》面向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", driver="GPKG")

逻辑说明：site_id是业务唯一键，重复点先按它去重，避免同一设施被重复统计。make_valid针对无效面最有用，点几何一般不受影响；preserve_topology不是这里必须的参数，因为make_valid的目标是得到合法几何。参数driver="GPKG"表示输出为 GeoPackage，适合保留 CRS 和属性。修复后不要直接覆盖原表，保留一份 clean 版本，后面排查有后悔药。

2.3 空间连接：sjoin 的 predicate、半径和候选集

空间连接是把“点属于哪个面”或“点附近有哪些点”落成表的关键步骤。predicate决定关系：within是点落在面内，intersects是相交，contains是包含。how="left"会保留未匹配的点，适合检查覆盖缺口；how="inner"只留匹配上的记录，适合做统计但容易掩盖边界问题。半径连接要盯住 CRS 单位，米制投影下max_distance=1000才是 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"]], how="left", predicate="within" ) near = gpd.sjoin_nearest( points, zones[["zone_id", "geometry"]], how="left", max_distance=1000, distance_col="dist_m" )

逻辑说明：inside里未匹配到zone_id的点，要么在边界外，要么面几何有问题，先查再删。sjoin_nearest的max_distance使用当前 CRS 的单位，6933 下是米；如果还在 4326 下，1000 会被当成度。distance_col保留实际距离，方便后续按 300 米、500 米、1000 米做敏感性对比。空间连接前先对两个图层做to_crs对齐，别指望库自动猜。

2.4 落库：geometry 与 geography 的边界

单机脚本跑通后，空间分析数据通常要落到空间数据库，供多人查询和增量更新。这里最常纠结的是geometry和geography：geometry快，适合投影坐标；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 或距离阈值。权重选错，Moran's I 和热点分析就会给出看似显著、实际不可解释的结果。这里要先问关系假设，再谈算法。

3.1 Queen、Rook、KNN 与距离阈值：先问关系假设

面数据的 Queen 邻接把共点也算邻居，Rook 只算共边；点数据更常用 KNN 或距离带。Queen 适合边界接触也算邻接的业务单元，Rook 适合只认共享边界的场景。KNN 的k不是越大越好，k 太小容易孤岛，k 太大把远处无关单元拉进来。距离阈值要有业务含义，比如 1 公里覆盖圈，而不是为了让结果显著随便调。权重还要决定是否行标准化：行标准化后每个单元的邻居权重和为 1，适合比较不同邻居数量的单元。

3.2 用 PySAL 建权重并检查孤岛

构建权重后第一件事不是跑统计，而是看孤岛。孤岛是没有邻居的单元，它们会让局部统计出现空值，也会让全局指标偏移。下面的代码对面数据建 Queen 权重，对点数据建 KNN 权重，并做一次 Moran's 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_index=True) queen.transform = "r" print("孤岛数量:", len(queen.islands)) points = gpd.read_file("data/site_points_clean.gpkg").to_crs(6933) knn = KNN.from_dataframe(points, k=6) knn.transform = "r" mi = Moran(zones["density"].values, queen, permutations=999) print(mi.I, mi.p_sim)

逻辑说明：use_index=True让权重矩阵使用 GeoDataFrame 的索引，后续对齐属性更稳。transform="r"是行标准化，适合邻居数量不一致的单元。permutations=999做随机置换检验，p_sim是模拟显著性；样本很小时不要只盯 p 值。k=6只是起点，实际要按业务半径和点密度调，调完记录参数，别让下一次跑批变成黑匣子。

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, geometry="geometry", crs=points.crs)

逻辑说明：grid_size=500表示 500 米方格，先除以格宽再取整，得到行列号。聚合时size统计点数量，mean统计指标均值；如果业务上更关心总和，把mean换成sum。生成格网中心点时乘回格宽并加半格，保证点落在格子中心。H3 和 GeoHash 的实现方式不同，但原则一样：先定层级或精度，再做多尺度对比，不要只交一个尺度。

3.4 热点分析：Getis-Ord Gi* 前先确认权重和样本量

热点分析不是给每个点打分，而是看局部高值或低值是否在空间上聚集。Getis-Ord Gi* 的star=True表示包含自身，适合找高值簇；star=False排除自身，解释会变。权重仍用 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, k=8) w.transform = "r" gi = G_Local( cells["point_count"].astype(float).values, w, star=True, permutations=999 ) cells["gi_z"] = gi.Zs cells["gi_p"] = gi.p_sim

逻辑说明：先过滤point_count >= 3，避免稀疏格网制造伪热点。k=8是邻域数量，格网越密可以适当增大，但要在报告里固定。star=True让自身参与计算，permutations=999做随机检验。gi_z是局部 z 分数，gi_p是模拟 p 值；解释时不要只看红色高值，还要看它是否通过显著性阈值，以及换一个格网尺寸后是否仍成立。

4. 空间分析数据翻车排查：五个高频现象和对应处置

空间分析数据的坑很少以报错形式出现，更多是结果看起来合理、实际已经偏了。下面五条是血泪经验里最常见的一类：空间连接丢记录、缓冲区单位错、查询慢、尺度反转、重复点污染。每条都按“现象、原因、解决”拆开，排查时先复现现象，再定位原因，最后用小样本验证修复。

4.1 现象：空间连接后大量空值；原因：几何无效或 CRS 不一致；解决：先对齐再修复

现象是sjoin之后zone_id空值比例异常高，但点位明明落在面内。原因通常有两个：两个图层 CRS 不一致，或者面几何自交、环方向异常导致within判断失败。解决顺序是先to_crs对齐，再用make_valid修复面几何，最后用how="left"保留未匹配记录逐条查。不要一上来就删空值，空值本身就是排查线索。修复后重新跑连接，并对比修复前后的匹配率，匹配率突增说明原几何确实有问题。

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 米，比较热点位置和 Moran's 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(f"data/tiles/points_{ix}_{iy}.gpkg", driver="GPKG")

逻辑说明：total_bounds给出整体范围，tile=50000表示 50 公里。intersects保留与瓦片相交的点，边界点可能同时落入两个瓦片，后续统计要按主键去重。输出文件名带ix、iy，方便重跑单个瓦片。参数不要照搬，点密集就缩小瓦片，面很大就放大瓦片，原则是单任务内存可控且不频繁切碎几何。

5.2 几何简化：拓扑简化与容差选择的边界

几何简化能减少传输和渲染成本，但简化后的几何不能直接用于面积统计。simplify(tolerance=10, preserve_topology=True)表示容差 10 米并保持拓扑，适合展示和粗筛；面积、长度仍用原始几何算。容差太大会让相邻面出现缝隙或重叠，空间连接结果随之改变。我的习惯是保留两列：原始几何用于计算，简化几何用于出图和前端加载，并在元数据里写清容差。

zones = gpd.read_file("data/zone_boundary.gpkg").to_crs(6933) zones["geom_simplified"] = zones.geometry.simplify( tolerance=10, preserve_topology=True ) zones["area_m2"] = zones.geometry.area # 面积仍用原始几何

逻辑说明：tolerance=10的单位是米，来自 6933；如果 CRS 是度，10 就完全不是一个量级。preserve_topology=True能减少自交和缝隙，但不能保证简化后完全满足业务拓扑，关键图层仍要做有效性检查。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 是否统一；第二项是属性校验：主键重复率、空间连接空值率、指标极值是否越界；第三项是统计校验：固定样本跑 Moran's I 或热点结果，和上一次跑批对比。三项都过，再谈出图。

固定样本不需要很大，抽几百条覆盖边界、中心、稀疏区即可。把样本主键、CRS、格网尺寸、权重参数写进跑批配置，下一次换人跑也能复现。我现在的习惯是：任何空间分析数据交付前，先问一句“换一台机器、换一个格网，这个结论还成立吗”，答不上来就先补校验。这个方向值得投入，但别把图画得漂亮当成结果可靠。希望帮到你。

本文还有配套的精品资源，点击获取