-- ⚠️ GCJ-02 和 BD-09 不是标准坐标系
-- PostGIS 没有内置的 SRID，不能直接用 ST_Transform 转换
-- 需要使用第三方算法或偏移修正

-- 偏移量示意（北京地区，约 300-600 米）
-- WGS84: (116.4074, 39.9042)
-- GCJ02: (116.4139, 39.9087)  -- 约偏移 500 米
-- BD09:  (116.4203, 39.9142)  -- 在 GCJ02 基础上再偏移

-- 查看 SRID 4326 的详细定义
SELECT srid, auth_name, auth_srid, srtext
FROM spatial_ref_sys WHERE srid = 4326;

-- 搜索 UTM Zone 50N（覆盖中国东部沿海 114°-120°E）
SELECT srid, auth_name, auth_srid, srtext
FROM spatial_ref_sys
WHERE srtext ILIKE '%utm%zone%50%'
  AND srtext ILIKE '%north%';

-- 搜索 CGCS2000 三度带
SELECT srid, srtext
FROM spatial_ref_sys
WHERE srtext ILIKE '%CGCS%2000%3%degree%'
LIMIT 10;

-- 速查表
WITH china_sr (srid, name, type, unit, coverage) AS (
    VALUES
        (4326,  'WGS 84',                              '地理', '度',   '全球'),
        (4490,  'CGCS 2000',                           '地理', '度',   '中国'),
        (4547,  'CGCS 2000 / 3-degree GK zone 39',    '投影', '米',   '117°带 (北京)'),
        (4548,  'CGCS 2000 / 3-degree GK zone 40',    '投影', '米',   '120°带 (上海)'),
        (4549,  'CGCS 2000 / 3-degree GK zone 41',    '投影', '米',   '123°带'),
        (4538,  'CGCS 2000 / 3-degree GK CM 117E',    '投影', '米',   '中央经线117°'),
        (4539,  'CGCS 2000 / 3-degree GK CM 120E',    '投影', '米',   '中央经线120°'),
        (32650, 'WGS 84 / UTM zone 50N',               '投影', '米',   '114°-120°E'),
        (32651, 'WGS 84 / UTM zone 51N',               '投影', '米',   '120°-126°E'),
        (3857,  'Web Mercator',                         '投影', '米',   '全球 Web 地图')
)
SELECT * FROM china_sr;

-- 将北京坐标从 WGS84 (4326) 转换为 CGCS2000 三度带 39 带 (4547)
SELECT
    ST_X(geom) AS wgs84_lng,
    ST_Y(geom) AS wgs84_lat,
    ST_X(ST_Transform(geom, 4547)) AS cgcs_x,
    ST_Y(ST_Transform(geom, 4547)) AS cgcs_y
FROM (SELECT ST_SetSRID(ST_MakePoint(116.4074, 39.9042), 4326) AS geom) t;

-- 结果:
-- wgs84_lng: 116.4074, wgs84_lat: 39.9042
-- cgcs_x: 447963.89,  cgcs_y: 4418554.07 (米)

-- 将整个表的数据从 4326 转换到 4547
UPDATE stores
SET geom = ST_Transform(geom, 4547);

-- ⚠️ 转换后记得更新 SRID
-- 方式 1: 用 ST_SetSRID 包装
UPDATE stores
SET geom = ST_SetSRID(ST_Transform(geom, 4547), 4547);

-- 方式 2: 使用 UpdateGeometrySRID（更新元数据）
SELECT UpdateGeometrySRID('stores', 'geom', 4547);
-- 注意: 这只更新元数据，不转换坐标！

-- 创建投影后的视图，原表保持 WGS84
CREATE VIEW stores_projected AS
SELECT
    id, name, address, store_type,
    ST_Transform(geom, 4547) AS geom
FROM stores;

-- 查询投影后的视图
SELECT name, ST_X(geom) AS x_m, ST_Y(geom) AS y_m
FROM stores_projected;

-- 反复转换可能导致精度损失
WITH point AS (
    SELECT ST_SetSRID(ST_MakePoint(116.4074, 39.9042), 4326) AS geom
)
SELECT
    ST_X(geom) AS original_lng,
    ST_X(ST_Transform(ST_Transform(geom, 4547), 4326)) AS roundtrip_lng,
    ST_X(geom) - ST_X(ST_Transform(ST_Transform(geom, 4547), 4326)) AS error_degrees
FROM point;
-- 误差通常在 1e-10 度以内，可忽略

-- ❌ 错误: 对已经是投影坐标的几何再次设置 4326
UPDATE table SET geom = ST_SetSRID(ST_MakePoint(447963, 4418554), 4326);
-- 这会创建一个经度 447963° 的"坐标"，完全错误！

-- ✅ 正确: 明确数据的实际 SRID
UPDATE table SET geom = ST_SetSRID(ST_MakePoint(447963, 4418554), 4547);

-- 计算指定经度对应的 CGCS2000 三度带号
-- 带号 = FLOOR((经度 - 1.5) / 3) + 25
-- SRID = 4500 + 带号

-- 示例：北京 (116.4074°E)
-- 带号 = FLOOR((116.4074 - 1.5) / 3) + 25 = FLOOR(38.3) + 25 = 63
-- SRID = 4500 + 63 - 16 = 4547

-- 创建函数自动计算
CREATE OR REPLACE FUNCTION get_cgcs2000_srid(lng DOUBLE PRECISION)
RETURNS INTEGER AS $$
BEGIN
    RETURN 4500 + FLOOR((lng - 1.5) / 3)::INTEGER + 25 - 16;
END;
$$ LANGUAGE plpgsql IMMUTABLE;

-- 测试
SELECT get_cgcs2000_srid(116.4074);  -- 北京: 4547
SELECT get_cgcs2000_srid(121.4737);  -- 上海: 4548
SELECT get_cgcs2000_srid(104.0668);  -- 成都: 4545

-- 创建测试: 北京到上海的距离
WITH points AS (
    SELECT
        ST_SetSRID(ST_MakePoint(116.4074, 39.9042), 4326) AS beijing,
        ST_SetSRID(ST_MakePoint(121.4737, 31.2304), 4326) AS shanghai
)
SELECT
    -- 方法 1: Geometry 距离 (度) — 无意义
    ROUND(ST_Distance(beijing, shanghai)::numeric, 4) AS dist_degrees,

    -- 方法 2: Geography 距离 (米) — 精确球面距离
    ROUND((ST_Distance(beijing::geography, shanghai::geography) / 1000)::numeric, 1) AS dist_geo_km,

    -- 方法 3: 投影后距离 (米) — 精确但取决于投影
    ROUND((ST_Distance(
        ST_Transform(beijing, 4547),
        ST_Transform(shanghai, 4548)
    ) / 1000)::numeric, 1) AS dist_proj_km,

    -- 方法 4: Haversine 公式 — 参考值
    ROUND((6371 * acos(
        cos(radians(39.9042)) * cos(radians(31.2304)) *
        cos(radians(121.4737) - radians(116.4074)) +
        sin(radians(39.9042)) * sin(radians(31.2304))
    ))::numeric, 1) AS dist_haversine_km
FROM points;

-- 结果对比:
-- dist_degrees:  ~10.77 (无意义)
-- dist_geo_km:   ~1068.1 (精确)
-- dist_proj_km:  ~1070.5 (投影误差)
-- dist_haversine: ~1068.0 (参考)

-- 添加自定义投影坐标系
-- 示例: 某城市的地方坐标系
INSERT INTO spatial_ref_sys (srid, auth_name, auth_srid, proj4text, srtext)
VALUES (
    920001,
    'CUSTOM',
    920001,
    '+proj=tmerc +lat_0=39.9042 +lon_0=116.4074 +k=1 +x_0=500000 +y_0=0 +datum=WGS84 +units=m +no_defs',
    'PROJCS["Beijing Local",GEOGCS["WGS 84",DATUM["WGS_1984",SPHEROID["WGS 84",6378137,298.257223563]],PRIMEM["Greenwich",0],UNIT["degree",0.0174532925199433]],PROJECTION["Transverse_Mercator"],PARAMETER["latitude_of_origin",39.9042],PARAMETER["central_meridian",116.4074],PARAMETER["scale_factor",1],PARAMETER["false_easting",500000],PARAMETER["false_northing",0],UNIT["metre",1]]'
);

-- 使用自定义坐标系
SELECT ST_Transform(
    ST_SetSRID(ST_MakePoint(116.4074, 39.9042), 4326),
    920001
);

-- 诊断表中几何的 SRID 分布
SELECT
    ST_SRID(geom) AS srid,
    count(*) AS count,
    ST_XMin(ST_Extent(geom)) AS xmin,
    ST_YMin(ST_Extent(geom)) AS ymin,
    ST_XMax(ST_Extent(geom)) AS xmax,
    ST_YMax(ST_Extent(geom)) AS ymax
FROM my_table
GROUP BY ST_SRID(geom);

-- 检查坐标范围是否合理
-- WGS84 经度: -180 ~ 180, 纬度: -90 ~ 90
SELECT
    CASE
        WHEN ST_XMin(ext) < -180 OR ST_XMax(ext) > 180 THEN '经度超出范围！'
        WHEN ST_YMin(ext) < -90 OR ST_YMax(ext) > 90 THEN '纬度超出范围！'
        ELSE '坐标范围正常'
    END AS diagnosis,
    ext
FROM (SELECT ST_Extent(geom) AS ext FROM my_table) t;