PostgreSQL中查询建筑点到指定值最近栅格像素的距离

背景 #

PostgreSQL数据库中存在两张表：

• hs_buildings：存储建筑点位置，核心字段coordinates（点几何，EPSG 2056坐标系）
• hs_raster_values：存储土地利用栅格数据，分辨率100x100，几何数据采用EPSG 2056坐标系

需求：计算建筑点到值为指定值（如170005）的栅格像素的最短距离。

原始低效实现 #

以下查询可实现需求，但因将所有目标栅格像素转为点再计算距离，全表运行时效率极低：

WITH hs_points AS
(
    WITH gv AS 
    (
        SELECT
            (ST_PixelAsCentroids(hs_raster_values.geometry, 1, true)).*
        FROM hs_raster_values
        WHERE hs_raster_values.raster_id = 1
    )

    SELECT
        (gv).x,
        (gv).y,
        (gv).val,
        gv.geom 
    FROM gv 
    WHERE val = 170005
)
SELECT 
    ST_Distance(hs_points.geom, hs_buildings.coordinates)
FROM hs_points, hs_buildings
WHERE hs_buildings.id = 1
ORDER BY hs_buildings.coordinates <-> hs_points.geom
LIMIT 1

使用PostGIS栅格函数ST_MinDist4ma优化 #

错误尝试及问题 #

直接调用ST_MinDist4ma会触发如下报错：

WITH building_pixel AS 
(
    SELECT 
        (ST_WorldToRasterCoordX(hs_raster_values.geometry, hs_buildings.coordinates),
        ST_WorldToRasterCoordY(hs_raster_values.geometry, hs_buildings.coordinates)) AS xy
    FROM hs_raster_values
    INNER JOIN hs_buildings ON ST_Intersects(hs_buildings.coordinates, hs_raster_values.geometry)
    WHERE hs_raster_values.raster_id = 1
    AND hs_buildings.egid = 1
    AND ST_Intersects(hs_buildings.coordinates, hs_raster_values.geometry)
)

SELECT 
    ST_MinDist4ma(170005, building_pixel.xy, hs_raster_values.geometry) AS min_distance
FROM hs_raster_values, building_pixel
WHERE raster_id = 1;

报错信息：

ERROR:  function st_mindist4ma(integer, record, raster) does not exist LINE 14:  ST_MinDist4ma(170005, building_pixel.xy, hs_raster_values.g...

原因：ST_MinDist4ma是移动窗口分析函数，不能直接独立调用，必须配合ST_MapAlgebraFctNgb使用。

正确用法 #

ST_MinDist4ma用于在指定邻域窗口内，计算目标像素到第一个匹配指定值像素的最短像素距离，再结合栅格分辨率转换为实际地理距离。

完整查询代码：

WITH building_raster_pos AS (
    SELECT
        hs_buildings.id AS building_id,
        ST_WorldToRasterCoordX(hs_raster_values.geometry, hs_buildings.coordinates) AS x,
        ST_WorldToRasterCoordY(hs_raster_values.geometry, hs_buildings.coordinates) AS y
    FROM hs_raster_values
    JOIN hs_buildings ON ST_Intersects(hs_buildings.coordinates, hs_raster_values.geometry)
    WHERE hs_raster_values.raster_id = 1
      AND hs_buildings.id = 1 -- 单个建筑计算，批量可移除该条件
)
SELECT
    brp.building_id,
    -- 像素距离转换为实际米数（栅格分辨率100m）
    (ST_MapAlgebraFctNgb(
        hs_raster_values.geometry,
        1, -- 分析的栅格列索引
        brp.x, brp.y, -- 建筑点对应的栅格行列号
        50, 50, -- 搜索窗口大小（可按需调整，如50x50像素范围）
        'ST_MinDist4ma(double precision[][][], integer[][], text[], double precision)',
        ARRAY[170005]::double precision[], -- 目标栅格值
        NULL, NULL,
        'FIRST' -- 匹配到第一个目标值即停止计算
    )) * 100 AS min_distance_meters
FROM hs_raster_values
JOIN building_raster_pos brp ON TRUE
WHERE hs_raster_values.raster_id = 1;

关键参数说明 #

1. ST_MapAlgebraFctNgb核心参数：
   • 第1个参数：目标栅格对象
   • 第2个参数：要分析的栅格列索引（通常为1）
   • 第3、4个参数：建筑点对应的栅格行列号（x,y）
   • 第5、6个参数：搜索窗口的宽高（比如50x50表示以目标像素为中心，搜索周围50像素范围）
   • 第7个参数：移动窗口函数的签名，必须与ST_MinDist4ma的参数匹配
   • 第8个参数：ST_MinDist4ma的输入参数（这里是目标栅格值170005）
2. 距离转换：因栅格分辨率为100m，将像素距离乘以100得到实际米数。

优化建议 #

• 空间索引：确保hs_buildings.coordinates和hs_raster_values.geometry都创建空间索引：

  CREATE INDEX idx_hs_buildings_coords ON hs_buildings USING GIST(coordinates);
  CREATE INDEX idx_hs_raster_geom ON hs_raster_values USING GIST(geometry);
  

• 窗口调整：根据目标栅格值的分布范围，调整搜索窗口大小（范围越小计算越快，需平衡准确性与性能）
• 批量计算：移除hs_buildings.id = 1条件，可一次性计算所有建筑的最短距离

内容的提问来源于stack exchange，提问作者Sybille Roemer