-- 启用拓扑扩展
CREATE EXTENSION IF NOT EXISTS postgis_topology;

-- 创建拓扑模式
-- 参数：拓扑名称, SRID, 容差（容差内的点视为同一个点）
SELECT topology.CreateTopology('admin_topo', 4326, 0.00001);

-- 查看拓扑信息
SELECT * FROM topology.topology WHERE name = 'admin_topo';

-- 假设已有行政区划表 districts (id, name, geom)
-- 创建拓扑层
SELECT topology.AddTopoGeometryColumn('admin_topo', 'public', 'districts', 'topo_geom', 'MULTIPOLYGON');

-- 将几何数据导入拓扑
UPDATE districts
SET topo_geom = topology.toTopoGeom(geom, 'admin_topo', 1);

-- 查看拓扑要素数量
SELECT
    'nodes' AS type, count(*) FROM admin_topo.node
UNION ALL
SELECT
    'edges' AS type, count(*) FROM admin_topo.edge_data
UNION ALL
SELECT
    'faces' AS type, count(*) FROM admin_topo.face;

-- 添加新的行政区划
INSERT INTO districts (name, topo_geom)
VALUES (
    '新城区',
    topology.toTopoGeom(
        ST_GeomFromText('POLYGON((116.5 39.9, 116.6 39.9, 116.6 40.0, 116.5 40.0, 116.5 39.9))', 4326),
        'admin_topo',
        1
    )
);

-- 移动共享边界（自动更新相邻面）
-- 步骤 1: 找到要编辑的边
SELECT edge_id, ST_AsText(geom)
FROM admin_topo.edge_data
WHERE ST_Intersects(geom, ST_GeomFromText('POINT(116.55 39.95)', 4326));

-- 步骤 2: 使用 ST_MoveIsoEdge 或直接编辑节点
-- 注意: 拓扑编辑通常通过拓扑编辑工具（如 QGIS Topology Editor）完成

-- 删除一个行政区划
DELETE FROM districts WHERE name = '新城区';

-- 清理孤立的拓扑对象
SELECT topology.ST_RemoveIsoNode('admin_topo', node_id)
FROM admin_topo.node n
WHERE NOT EXISTS (
    SELECT 1 FROM admin_topo.relation r
    WHERE r.topogeo_id = n.node_id AND r.layer_id = 1
);

-- 从拓扑要素获取几何
SELECT
    d.name,
    d.topo_geom,
    topology.ST_GetFaceGeometry('admin_topo', face_id) AS face_geom
FROM districts d;

-- 拓扑关系查询
SELECT a.name AS district_a, b.name AS district_b
FROM districts a, districts b
WHERE a.id < b.id
  AND ST_Intersects(
      topology.ST_GetFaceGeometry('admin_topo', (a.topo_geom).topogeo_id),
      topology.ST_GetFaceGeometry('admin_topo', (b.topo_geom).topogeo_id)
  );

-- 查看两个相邻面共享的边
SELECT
    e.edge_id,
    ST_AsText(e.geom) AS edge_geom,
    e.start_node,
    e.end_node
FROM admin_topo.edge_data e
WHERE EXISTS (
    SELECT 1 FROM admin_topo.relation r1
    WHERE r1.element_id = e.edge_id AND r1.element_type = 2
      AND r1.topogeo_id = (SELECT (topo_geom).topogeo_id FROM districts WHERE name = '东城区')
)
AND EXISTS (
    SELECT 1 FROM admin_topo.relation r2
    WHERE r2.element_id = e.edge_id AND r2.element_type = 2
      AND r2.topogeo_id = (SELECT (topo_geom).topogeo_id FROM districts WHERE name = '西城区')
);

-- 检查同层面要素是否重叠
WITH pairs AS (
    SELECT
        a.id AS id_a, a.name AS name_a, a.geom AS geom_a,
        b.id AS id_b, b.name AS name_b, b.geom AS geom_b
    FROM districts a, districts b
    WHERE a.id < b.id
)
SELECT
    name_a, name_b,
    ST_Area(ST_Intersection(geom_a, geom_b)::geography) / 1000000 AS overlap_km2,
    ST_Intersection(geom_a, geom_b) AS intersection_geom
FROM pairs
WHERE ST_Overlaps(geom_a, geom_b);

-- 检查面要素之间的间隙
WITH union_geom AS (
    SELECT ST_Union(geom) AS geom FROM districts
),
boundary AS (
    SELECT ST_Boundary(
        ST_Envelope(geom)
    ) AS geom FROM union_geom
)
SELECT ST_Difference(b.geom, u.geom) AS gap_geom
FROM boundary b, union_geom u
WHERE NOT ST_IsEmpty(ST_Difference(b.geom, u.geom));

-- 修复重叠：将重叠部分分配给面积较大的要素
CREATE OR REPLACE FUNCTION fix_overlaps() RETURNS void AS $$
DECLARE
    rec RECORD;
BEGIN
    FOR rec IN
        SELECT a.id AS id_a, b.id AS id_b,
               ST_Intersection(a.geom, b.geom) AS overlap_geom
        FROM districts a, districts b
        WHERE a.id < b.id AND ST_Overlaps(a.geom, b.geom)
    LOOP
        -- 从较小的要素中减去重叠部分
        UPDATE districts
        SET geom = ST_Difference(geom, rec.overlap_geom)
        WHERE id = rec.id_b;
    END LOOP;
END;
$$ LANGUAGE plpgsql;

-- 验证拓扑的有效性
SELECT * FROM topology.ValidateTopology('admin_topo');

-- 验证结果包含:
-- 1. 无效的边（自相交等）
-- 2. 无效的面（不闭合等）
-- 3. 孤立的节点/边
-- 4. 不一致的关系

-- 验证单个拓扑几何的有效性
SELECT
    name,
    topology.ST_GetFaceGeometry('admin_topo', (topo_geom).topogeo_id) AS geom,
    ST_IsValid(topology.ST_GetFaceGeometry('admin_topo', (topo_geom).topogeo_id)) AS is_valid
FROM districts;

-- 导出为普通几何表
CREATE TABLE districts_simple AS
SELECT
    id,
    name,
    (topo_geom)::geometry AS geom
FROM districts;

-- 创建空间索引
CREATE INDEX idx_districts_simple_geom ON districts_simple USING GIST(geom);

-- 清空拓扑后重新导入
-- 1. 清空
SELECT topology.ClearTopoGeomColumn('public', 'districts', 'topo_geom');

-- 2. 删除旧拓扑
SELECT topology.DropTopology('admin_topo');

-- 3. 重新创建
SELECT topology.CreateTopology('admin_topo', 4326, 0.00001);
SELECT topology.AddTopoGeometryColumn('admin_topo', 'public', 'districts', 'topo_geom', 'MULTIPOLYGON');

-- 4. 重新导入
UPDATE districts
SET topo_geom = topology.toTopoGeom(geom, 'admin_topo', 1);

-- 创建地籍拓扑
SELECT topology.CreateTopology('cadastral_topo', 4490, 0.001);

CREATE TABLE land_parcels (
    parcel_id SERIAL PRIMARY KEY,
    owner_name TEXT,
    land_use VARCHAR(50),
    area_m2 NUMERIC(12,2),
    topo_geom topology.topogeometry
);

SELECT topology.AddTopoGeometryColumn('cadastral_topo', 'public', 'land_parcels', 'topo_geom', 'MULTIPOLYGON');

-- 导入地籍数据
INSERT INTO land_parcels (owner_name, land_use, topo_geom)
SELECT
    owner_name,
    land_use,
    topology.toTopoGeom(geom, 'cadastral_topo', 1)
FROM raw_parcels;

-- 查询相邻地块
WITH parcel AS (
    SELECT topo_geom FROM land_parcels WHERE parcel_id = 100
)
SELECT lp.parcel_id, lp.owner_name
FROM land_parcels lp, parcel p
WHERE lp.parcel_id != 100
  AND ST_Touches(
      (lp.topo_geom)::geometry,
      (p.topo_geom)::geometry
  );