在指定坐标上将PNG叠加到GeoTIFF上

Question

我有一个带地理信息的tiff文件，gdalinfo的输出如下：

Driver: GTiff/GeoTIFF
Files: generated.tiff
       generated.tiff.aux.xml
Size is 6941, 4886
Coordinate System is `'
GCP Projection = 
GEOGCS["WGS 84",
    DATUM["WGS_1984",
        SPHEROID["WGS 84",6378137,298.257223563,
            AUTHORITY["EPSG","7030"]],
        AUTHORITY["EPSG","6326"]],
    PRIMEM["Greenwich",0],
    UNIT["degree",0.0174532925199433],
    AUTHORITY["EPSG","4326"]]
GCP[  0]: Id=1, Info=
          (0,0) -> (0.01,0.05886,0)
GCP[  1]: Id=2, Info=
          (6941,0) -> (0.07941,0.05886,0)
GCP[  2]: Id=3, Info=
          (6941,4886) -> (0.07941,0.01,0)
GCP[  3]: Id=4, Info=
          (0,4886) -> (0.01,0.01,0)
Metadata:
  AREA_OR_POINT=Area
  Software=paint.net 4.0
Image Structure Metadata:
  INTERLEAVE=BAND
Corner Coordinates:
Upper Left  (    0.0,    0.0)
Lower Left  (    0.0, 4886.0)
Upper Right ( 6941.0,    0.0)
Lower Right ( 6941.0, 4886.0)
Center      ( 3470.5, 2443.0)

还有一个文件是地图标记的图片，叫做marker1.png（大小是36x60像素）。

我想把marker1.png叠加到上面生成的.tiff文件上，让它的左上角位于这个地理tiff文件的坐标0.037,0.025。这样看起来就像在谷歌地图上放了一个标记。

我该怎么做才能实现这个呢？

我已经部分实现了这个功能，但不确定这是不是正确的方向。

import gdal

gdal.UseExceptions()
s = gdal.Open('generated.tiff')

drv = gdal.GetDriverByName("VRT")
vrt = drv.CreateCopy('test.vrt', s, 0)
band = vrt.GetRasterBand(1)

source_path = 'marker1.png'
source_band = 1
x_size = 36
y_size = 60
x_block = 36
y_block = 1
x_offset = 0
y_offset = 0
x_source_size = 36
y_source_size = 60
dest_x_offset = 2000
dest_y_offset = 2000
x_dest_size = 36
y_dest_size = 60

simple_source = '<SimpleSource><SourceFilename relativeToVRT="1">%s</SourceFilename>' % source_path + \
    '<SourceBand>%i</SourceBand>' % source_band + \
    '<SourceProperties RasterXSize="%i" RasterYSize="%i" DataType="Byte" BlockXSize="%i" BlockYSize="%i"/>' % (x_size, y_size, x_block, y_block) + \
    '<SrcRect xOff="%i" yOff="%i" xSize="%i" ySize="%i"/>' % (x_offset, y_offset, x_source_size, y_source_size) + \
    '<DstRect xOff="%i" yOff="%i" xSize="%i" ySize="%i"/></SimpleSource>' % (dest_x_offset, dest_y_offset, x_dest_size, y_dest_size)
band.SetMetadata({'source_0': simple_source}, "new_vrt_sources")
band.SetMetadataItem("NoDataValue", '1')

p = gdal.GetDriverByName("PNG")
p.CreateCopy('result.png', vrt, 0)

vrt = None

这个方法使用的是像素坐标，而不是地理坐标（不过转换很简单），但是标记图片显示成了黑色的块（虽然尺寸是对的）——看起来调色板可能有问题？

Answer 1

我尝试了很多不同的方法，但都没有成功——颜色不对、透明度不对或者其他错误。

最后我用PIL这个库解决了问题，下面的代码就是我的做法。代码很简短，容易理解（比起我用gdal想出来的那些方法），最重要的是——它能正常工作。

当然，这个方法还可以进一步改进。

from PIL import Image, ImageFont, ImageDraw
from osgeo import gdal,ogr

image = 'generated.tiff'
src_ds = gdal.Open(image)
gt = src_ds.GetGeoTransform() # used to convert geographical coordinates to pixel coordinates

font = ImageFont.truetype("sans-serif.ttf", 16)

img = Image.open(image)

def add_marker (gt, watermark, font, img, mx, my, text):
    px = int((mx - gt[0]) / gt[1]) #x pixel
    py = int((my - gt[3]) / gt[5]) #y pixel


    wmark = Image.open(watermark)
    draw = ImageDraw.Draw(wmark)
    draw.text((12, 10), text, (0, 0, 0), font=font)

    img.paste(wmark, (px, py), wmark)

add_marker(gt, 'marker1.png', font, img, 0.012, 0.0132, "1")

img.save("result.png", "PNG")