-- 安装栅格扩展
CREATE EXTENSION IF NOT EXISTS postgis_raster;

-- 验证
SELECT PostGIS_Raster_Lib_Version();

# 基本导入
raster2pgsql -s 4326 -I -C -M /data/dem.tif public.dem | psql -d gisdb

# 参数说明:
# -s SRID        坐标系
# -I             创建空间索引
# -C             添加约束（范围、SRID等）
# -M             运行 VACUUM ANALYZE
# -t TILE_SIZE   分块大小（如 256x256）
# -l OVERVIEW_LEVEL 创建概视图层级

# 分块导入（推荐用于大文件）
raster2pgsql -s 4326 -I -C -M -t 256x256 /data/dem.tif public.dem_tiles | psql -d gisdb

# 多波段影像（如 RGB 卫星影像）
raster2pgsql -s 4326 -I -C -t 256x256 -b 1 -b 2 -b 3 /data/sentinel2.tif public.sentinel | psql -d gisdb

# 多文件批量导入
for tif in /data/tiles/*.tif; do
    raster2pgsql -s 4326 -I -C -a -t 256x256 "$tif" public.dem_tiles | psql -d gisdb
done

-- 从已知参数创建空栅格
SELECT ST_MakeEmptyRaster(
    100,      -- 宽度（像素数）
    100,      -- 高度
    0, 0,     -- 左上角坐标 (X, Y)
    10,       -- 像素宽度（分辨率）
    -10,      -- 像素高度（负值表示从上到下）
    0,        -- 旋转
    0,        -- 旋转
    4326      -- SRID
);

-- 从现有栅格创建同规格空栅格
SELECT ST_MakeEmptyRaster(rast) FROM dem LIMIT 1;

-- 从外部文件创建
SELECT ST_FromGDALRaster(pg_read_binary_file('/data/small_dem.tif'));

-- 获取栅格元数据
SELECT
    (md).width AS pixel_width,
    (md).height AS pixel_height,
    (md).upperleftx AS ul_x,
    (md).upperlefty AS ul_y,
    (md).scalex AS pixel_size_x,
    (md).scaley AS pixel_size_y,
    (md).numbands AS num_bands,
    (md).srid
FROM (SELECT ST_Metadata(rast) AS md FROM dem LIMIT 1) t;

-- 查看波段信息
SELECT
    (bd).nodata,
    (bd).pixeltype,
    (bd).hasnodata,
    (bd).isoutdb,
    (bd).path
FROM (SELECT ST_BandMetadata(rast, 1) AS bd FROM dem LIMIT 1) t;

-- 提取指定点的像素值
SELECT ST_Value(rast, 1, ST_SetSRID(ST_MakePoint(116.4, 39.9), 4326)) AS elevation
FROM dem
WHERE ST_Intersects(rast, ST_SetSRID(ST_MakePoint(116.4, 39.9), 4326));

-- 批量提取多个点的值
WITH points AS (
    SELECT id, ST_SetSRID(ST_MakePoint(lng, lat), 4326) AS geom
    FROM sample_points
)
SELECT
    p.id,
    ST_Value(d.rast, 1, p.geom) AS elevation
FROM points p
JOIN dem d ON ST_Intersects(d.rast, p.geom);

-- 按区域裁剪栅格
SELECT ST_Clip(rast, 1,
    ST_GeomFromText('POLYGON((116.3 39.8, 116.5 39.8, 116.5 40.0, 116.3 40.0, 116.3 39.8))', 4326)
) AS clipped_rast
FROM dem
WHERE ST_Intersects(rast,
    ST_GeomFromText('POLYGON((116.3 39.8, 116.5 39.8, 116.5 40.0, 116.3 40.0, 116.3 39.8))', 4326)
);

-- 将高程从米转换为英尺
SELECT ST_MapAlgebra(rast, 1,
    '[rast] * 3.28084'::text,
    '32BF'  -- 输出像素类型
) AS rast_feet
FROM dem;

-- 高程分级
SELECT ST_MapAlgebra(rast, 1,
    $$
    CASE
        WHEN [rast] < 100 THEN 1      -- 平原
        WHEN [rast] < 500 THEN 2      -- 丘陵
        WHEN [rast] < 1000 THEN 3     -- 低山
        WHEN [rast] < 3000 THEN 4     -- 中山
        ELSE 5                         -- 高山
    END
    $$::text,
    '8BUI'  -- 8位无符号整数
) AS rast_classified
FROM dem;

-- 计算两个时期的 NDVI 差异
SELECT ST_MapAlgebra(
    rast_2024, 1,
    rast_2020, 1,
    '[rast1] - [rast2]'::text,
    '32BF'
) AS ndvi_change
FROM ndvi_pairs;

-- 洪水淹没分析（高程 < 水位的区域）
SELECT ST_MapAlgebra(
    dem.rast, 1,
    ST_AddBand(ST_MakeEmptyRaster(dem.rast), '8BUI'::text, 0, -1),
    $$
    CASE
        WHEN [rast1] < 50.0 THEN 1   -- 被淹没
        ELSE 0                         -- 未淹没
    END
    $$::text
) AS flood_rast
FROM dem
WHERE ST_Intersects(dem.rast, study_area);

-- 使用自定义 PL/pgSQL 函数作为回调
CREATE OR REPLACE FUNCTION calc_slope(
    pixel_value DOUBLE PRECISION,
    pos INTEGER[],
    VARIADIC userargs TEXT[]
) RETURNS DOUBLE PRECISION AS $$
BEGIN
    -- 简化的坡度计算（实际需要邻域像素）
    IF pixel_value IS NULL THEN RETURN NULL; END IF;
    RETURN ABS(pixel_value) * 0.01;  -- 示例：简单转换
END;
$$ LANGUAGE plpgsql IMMUTABLE;

SELECT ST_MapAlgebra(rast, 1, 'calc_slope(double precision, integer[], text[])'::regprocedure, '32BF')
FROM dem;

-- 计算坡度（度）
SELECT ST_Slope(
    rast, 1,
    '32BF',   -- 输出类型
    'DEGREES'  -- 单位：度
) AS slope_rast
FROM dem;

-- 提取坡度值到点
SELECT ST_Value(
    ST_Slope(rast, 1, '32BF', 'DEGREES'),
    1,
    ST_SetSRID(ST_MakePoint(116.4, 39.9), 4326)
) AS slope_degrees
FROM dem
WHERE ST_Intersects(rast, ST_SetSRID(ST_MakePoint(116.4, 39.9), 4326));

-- 计算坡向（度，0=北, 90=东, 180=南, 270=西）
SELECT ST_Aspect(
    rast, 1,
    '32BF',
    'DEGREES'
) AS aspect_rast
FROM dem;

-- 生成山体阴影（Hillshade）
SELECT ST_Hillshade(
    rast, 1,
    '32BF',
    315,  -- 方位角（光源方向）
    45    -- 高度角（光源高度）
) AS hillshade_rast
FROM dem;

-- 计算每个行政区内的平均高程
SELECT
    d.name AS district_name,
    (stats).mean AS avg_elevation,
    (stats).min AS min_elevation,
    (stats).max AS max_elevation,
    (stats).stddev AS elevation_stddev
FROM districts d
CROSS JOIN LATERAL (
    SELECT ST_SummaryStats(ST_Clip(e.rast, d.geom)) AS stats
    FROM dem e
    WHERE ST_Intersects(e.rast, d.geom)
) s;

-- 全局统计
SELECT
    (stats).count AS pixel_count,
    (stats).sum AS total,
    (stats).mean AS average,
    (stats).stddev AS standard_deviation,
    (stats).min AS minimum,
    (stats).max AS maximum
FROM (SELECT ST_SummaryStats(rast) AS stats FROM dem) t;

-- 高程直方图
SELECT
    (hist).value AS elevation,
    (hist).count AS pixel_count
FROM (
    SELECT ST_Histogram(rast, 1, 10) AS hist
    FROM dem
) t;

-- 将栅格多边形化（转为矢量面）
SELECT
    (gv).val AS elevation,
    (gv).geom AS geom
FROM (
    SELECT ST_DumpAsPolygons(ST_Clip(rast, study_area)) AS gv
    FROM dem
    WHERE ST_Intersects(rast, study_area)
) t;

-- 将矢量面栅格化
SELECT ST_AsRaster(
    geom,           -- 矢量几何
    1000, 1000,     -- 宽度、高度（像素）
    '8BUI',         -- 像素类型
    255,            -- 填充值
    0               -- 无数据值
) AS rast
FROM districts
WHERE name = '北京市';

-- 创建 XYZ 瓦片函数
CREATE OR REPLACE FUNCTION generate_tile(
    z INTEGER, x INTEGER, y INTEGER
) RETURNS BYTEA AS $$
DECLARE
    tile BYTEA;
BEGIN
    SELECT ST_AsPNG(ST_ColorMap(
        ST_Clip(rast, tile_geom),
        1,
        '100% 0 128 0, 50% 255 255 0, 0% 255 0 0'
    ))
    INTO tile
    FROM dem,
         ST_TileEnvelope(z, x, y) AS tile_geom
    WHERE ST_Intersects(rast, tile_geom);

    RETURN tile;
END;
$$ LANGUAGE plpgsql STABLE;

-- 调用示例
SELECT generate_tile(10, 512, 256);

-- 为大栅格创建多级概视图，加速缩略图显示
SELECT ST_CreateOverview('public.dem'::regclass, 'rast', 2);
SELECT ST_CreateOverview('public.dem'::regclass, 'rast', 4);
SELECT ST_CreateOverview('public.dem'::regclass, 'rast', 8);

-- 查看栅格约束
SELECT
    r_table_name, r_raster_column,
    srid, scale_x, scale_y,
    blocksize_x, blocksize_y,
    same_alignment, regular_blocking,
    num_bands, pixel_types, nodata_values
FROM raster_columns
WHERE r_table_name = 'dem';

-- 添加对齐约束（确保所有瓦片对齐）
SELECT AddRasterConstraints('public'::name, 'dem'::name, 'rast'::name);