如何在PostGIS中将30m多边形重采样为1km分辨率多边形？

不需要转栅格，直接用PostGIS的矢量空间聚合就能实现需求，步骤如下：

1. 核心思路 #

因为SRID 4326是经纬度坐标系，无法直接用米单位计算格网大小，所以先将数据转换为米单位的投影坐标系（如UTM分带投影），生成1km菱形格网后，再关联原30m格网计算均值，最后转回4326坐标系。

2. 完整SQL代码 #

-- 替换为你数据所在区域的UTM投影EPSG代码，比如北半球33带用32633，根据实际调整
WITH projected_bbox AS (
    -- 将原数据的边界转换为米单位投影
    SELECT ST_Transform(ST_Extent(geom), 32633) AS bbox FROM my_data
),
grid_params AS (
    -- 定义1km菱形参数：这里取对角线长度为1000m，可根据需求改为边长1000m
    SELECT
        1000 AS diamond_diagonal,
        ST_XMin(bbox) AS min_x,
        ST_YMin(bbox) AS min_y,
        ST_XMax(bbox) AS max_x,
        ST_YMax(bbox) AS max_y
    FROM projected_bbox
),
grid_points AS (
    -- 生成1km菱形的网格原点，间隔为菱形对角线长度
    SELECT
        ST_MakePoint(
            min_x + (diamond_diagonal * x),
            min_y + (diamond_diagonal * y)
        ) AS point,
        x, y
    FROM grid_params,
         generate_series(0, ceil((max_x - min_x)/diamond_diagonal)::int) x,
         generate_series(0, ceil((max_y - min_y)/diamond_diagonal)::int) y
),
diamond_grid AS (
    -- 从每个原点生成菱形多边形，再转回4326坐标系
    SELECT
        ST_Transform(
            ST_MakePolygon(
                ST_MakeLine(
                    ARRAY[
                        ST_Translate(point, 0, diamond_diagonal/2),
                        ST_Translate(point, diamond_diagonal/2, 0),
                        ST_Translate(point, 0, -diamond_diagonal/2),
                        ST_Translate(point, -diamond_diagonal/2, 0),
                        ST_Translate(point, 0, diamond_diagonal/2)
                    ]
                )
            ),
            4326
        ) AS geom,
        x, y
    FROM grid_points, grid_params
),
aggregated_result AS (
    -- 关联原30m格网，计算每个1km菱形内的val均值
    SELECT
        dg.geom,
        ROUND(AVG(md.val)::numeric, 2) AS avg_val,
        dg.x, dg.y
    FROM diamond_grid dg
    JOIN my_data md ON ST_Intersects(dg.geom, md.geom)
    GROUP BY dg.geom, dg.x, dg.y
)
-- 将结果存入新表，也可直接查询使用
SELECT * INTO my_1km_diamond_grid FROM aggregated_result;

3. 关键细节说明 #

• 投影坐标系选择：必须替换代码中的32633为你数据所在区域的UTM分带EPSG代码（可通过查询原数据中心点的经纬度确定，比如东经10-15度用32631）。如果不需要高精度，也可以用Web墨卡托（EPSG:3857），但高纬度区域会有变形。
• 菱形参数调整：如果需要的是边长为1km的菱形，将diamond_diagonal改为1000*sqrt(2)（因为菱形对角线长度=边长*√2），或者直接调整ST_Translate的偏移值为500（半边长）。
• 对齐原菱形格网：如果原30m菱形有特定旋转角度，可先计算原菱形的旋转角，再用ST_Rotate对生成的1km菱形进行旋转，确保格网对齐：

  -- 计算原菱形的旋转角度（取第一个多边形的边方向）
  WITH first_diamond AS (SELECT geom FROM my_data LIMIT 1),
  rotation AS (SELECT degrees(ST_Azimuth(ST_PointN(geom,1), ST_PointN(geom,2))) AS theta FROM first_diamond)
  -- 在diamond_grid中添加旋转：ST_Rotate(..., radians(theta))
  

4. 关于你尝试的栅格方法补充 #

你之前的栅格创建错误在于：SRID 4326的单位是度，ST_MakeEmptyRaster的像素大小参数填30代表30度，完全不符合30m的需求。如果一定要用栅格方法，必须先转投影到米单位坐标系，再创建30m分辨率栅格，用ST_AsRaster填充值后，用ST_Resample重采样到1km，最后再转成矢量格网，但步骤比矢量聚合复杂。

内容的提问来源于stack exchange，提问作者metal bar