如何修改GDAL命令实现多栅格堆叠且避免Alpha覆盖底层数据？

原命令的核心问题有两个：

1. gdal_rasterize的-init 1参数给所有波段（包括Alpha）设置了初始值1，导致背景Alpha接近透明但并非完全透明，叠加时会干扰底层显示；同时生成的两个栅格都是4波段（RGB+Alpha）。
2. gdalbuildvrt默认将输入栅格的波段按顺序拼接，最终生成8波段的VRT，gdal_translate输出的8波段文件在多数软件中仅显示最后4波段（r2.tiff），并非预期的图层叠加效果。

步骤1：正确生成带Alpha通道的栅格 #

确保背景Alpha为0（完全透明），要素Alpha为255（完全不透明）：

# 生成r1.tiff：黄色要素，透明背景
gdal_rasterize -burn 255 -burn 255 -burn 0 -burn 255 \
  -l korkeusviiva \
  -init 0 -init 0 -init 0 -init 0 \
  -ts 512 512 -ot byte -te x1 y1 x2 y2 \
  dat.gpkg r1.tiff

# 生成r2.tiff：品红色要素，透明背景（修正原命令的dat.pkgp笔误为dat.gpkg）
gdal_rasterize -burn 255 -burn 0 -burn 255 -burn 255 \
  -l tieviiva \
  -init 0 -init 0 -init 0 -init 0 \
  -ts 512 512 -ot byte -te x1 y1 x2 y2 \
  dat.gpkg r2.tiff

步骤2：实现图层叠加（二选一） #

方法A：用gdalwarp完成叠加 #

先以r1为基础，再叠加r2，最终输出4波段叠加结果：

# 先创建带Alpha的基础图层
gdalwarp -te x1 y1 x2 y2 -ts 512 512 -ot byte -dstalpha r1.tiff temp_base.tiff

# 将r2叠加到基础图层上
gdalwarp -te x1 y1 x2 y2 -ts 512 512 -ot byte -overwrite temp_base.tiff r2.tiff merge.tiff

# 清理临时文件
rm temp_base.tiff

方法B：用自定义VRT+gdal_translate #

创建merge.vrt文件（替换<GeoTransform>为r1.tiff的实际地理变换，可通过gdalinfo r1.tiff获取）：

<VRTDataset rasterXSize="512" rasterYSize="512">
  <GeoTransform>x1_res, res_x, 0, y2_res, 0, -res_y</GeoTransform>
  <VRTRasterBand dataType="Byte" band="1">
    <ColorInterp>Red</ColorInterp>
    <ComplexSource>
      <SourceFilename relativeToVRT="1">r2.tiff</SourceFilename>
      <SourceBand>1</SourceBand>
      <AlphaSourceBand>4</AlphaSourceBand>
      <NODATA>0</NODATA>
    </ComplexSource>
    <ComplexSource>
      <SourceFilename relativeToVRT="1">r1.tiff</SourceFilename>
      <SourceBand>1</SourceBand>
      <AlphaSourceBand>4</AlphaSourceBand>
      <NODATA>0</NODATA>
    </ComplexSource>
  </VRTRasterBand>
  <VRTRasterBand dataType="Byte" band="2">
    <ColorInterp>Green</ColorInterp>
    <ComplexSource>
      <SourceFilename relativeToVRT="1">r2.tiff</SourceFilename>
      <SourceBand>2</SourceBand>
      <AlphaSourceBand>4</AlphaSourceBand>
      <NODATA>0</NODATA>
    </ComplexSource>
    <ComplexSource>
      <SourceFilename relativeToVRT="1">r1.tiff</SourceFilename>
      <SourceBand>2</SourceBand>
      <AlphaSourceBand>4</AlphaSourceBand>
      <NODATA>0</NODATA>
    </ComplexSource>
  </VRTRasterBand>
  <VRTRasterBand dataType="Byte" band="3">
    <ColorInterp>Blue</ColorInterp>
    <ComplexSource>
      <SourceFilename relativeToVRT="1">r2.tiff</SourceFilename>
      <SourceBand>3</SourceBand>
      <AlphaSourceBand>4</AlphaSourceBand>
      <NODATA>0</NODATA>
    </ComplexSource>
    <ComplexSource>
      <SourceFilename relativeToVRT="1">r1.tiff</SourceFilename>
      <SourceBand>3</SourceBand>
      <AlphaSourceBand>4</AlphaSourceBand>
      <NODATA>0</NODATA>
    </ComplexSource>
  </VRTRasterBand>
  <VRTRasterBand dataType="Byte" band="4">
    <ColorInterp>Alpha</ColorInterp>
    <ComplexSource>
      <SourceFilename relativeToVRT="1">r2.tiff</SourceFilename>
      <SourceBand>4</SourceBand>
      <NODATA>0</NODATA>
    </ComplexSource>
    <ComplexSource>
      <SourceFilename relativeToVRT="1">r1.tiff</SourceFilename>
      <SourceBand>4</SourceBand>
      <NODATA>0</NODATA>
    </ComplexSource>
  </VRTRasterBand>
</VRTDataset>

然后执行转换：

gdal_translate merge.vrt merge.tiff

核心是利用GDAL的GDALDataset和GDALRasterBand进行逐像素的Alpha合成：

1. 打开r1和r2两个数据集，获取各自的RGB+Alpha波段。
2. 创建输出数据集，设置4波段（Byte类型），地理范围、分辨率与输入一致。
3. 逐块读取r1和r2的像素数据，对每个像素执行Alpha混合：
   • 若r2的Alpha值为255（完全不透明），直接使用r2的RGB值和Alpha值。
   • 若r2的Alpha值为0（完全透明），使用r1的RGB值和Alpha值。
   • 中间Alpha值则按公式计算：out_r = (r2_r * r2_a + r1_r * (255 - r2_a)) / 255，同理计算G、B通道，Alpha通道取max(r1_a, r2_a)或按混合公式计算。
4. 将计算后的像素数据写入输出数据集。

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