在PostGIS中进行点数据的栅格化（Rasterization of Point Data in PostGIS）

说明

介绍在PotGIS中将点数据转换为栅格数据。

关键字： raster、point、PostGIS

环境准备

• Postgresql版本：PostgreSQL 14.0, 64-bit
• PostGIS版本：POSTGIS="3.3.2"
• QGIS版本：3.28.3-Firenze

基本步骤

一、数据准备

测试数据中有一张点数据表，坐标系3857。

CREATE TABLE IF NOT EXISTS public.test_point
(
    geom geometry(Point,3857),
    id bigint,
    h numeric,--指定要转换为栅格值的列，如高程、温度、风速等
    CONSTRAINT test_point_pkey PRIMARY KEY (id)
)

TABLESPACE pg_default;

ALTER TABLE public.test_point
    OWNER to postgres;

二、二维点数据转换为三维点

DROP TABLE IF EXISTS test_point_3d;
CREATE TABLE test_point_3d
    as select h,ST_Force3D(ST_Transform(geom,3857),tp.h) as geom 
    from test_point tp;       

三、将三维点数据转换为栅格

WITH inputs AS 
    (SELECT 16::INTEGER AS pixelsize,'linear:radius:100' AS algorithm, 
    ST_Collect(geom) AS geom, ST_Expand(ST_Collect(geom), 10) AS ext
    FROM test_point_3d ),sizes AS 
    (SELECT ceil((ST_XMax(ext) - ST_XMin(ext))/pixelsize)::integer AS width,
         ceil((ST_YMax(ext) - ST_YMin(ext))/pixelsize)::integer AS height,
         ST_XMin(ext) AS upperleftx,ST_YMax(ext) AS upperlefty
    FROM inputs )SELECT 1 AS rid,
         ST_InterpolateRaster( geom,
         algorithm,
         ST_SetSRID(ST_AddBand(ST_MakeEmptyRaster(width,
         height,upperleftx,upperlefty,pixelsize),
         '32BF'), ST_SRID(geom)) ) AS rast FROM sizes, inputs;

其中ST_InterpolateRaster 参见:

https://postgis.net/docs/RT_ST_InterpolateRaster.html

其参数algorithm参见：

https://gdal.org/programs/gdal_grid.html#interpolation-algorithms