已知一个GeoTiff，得到其边界矢量的方法

在制作GeoTiff的样本时，经常需要知道它的边界矢量，原因如下：

1. 没有原始的裁切矢量；
2. 裁切它的矢量栅格化后与影像相差一个像素；

代码摘要：

1. 制作一个与原影像同坐标系但是值全部为0的中间影像A；
2. 将A矢量化

不足：

1. 没有坐标系，需要手动添加

附代码

import gdal
import ogr
import sys
import cv2
import numpy as np
import os
import osr


def main(argv):
    if len(argv) < 2:
        print('[target image] [shpfile name]')
        sys.exit(0)

    # this allows GDAL to throw Python Exceptions
    gdal.UseExceptions()

    #  get raster datasource
    src_DS = gdal.Open(argv[1])
    if src_DS is None:
        print('Unable to open %s' % src_filename)
        sys.exit(1)

    # create a temp tif file to simplify polygonize
    xsize = src_DS.RasterXSize
    ysize = src_DS.RasterYSize
    arr = np.zeros([ysize, xsize])
    # in case of 'tif' or 'tiff'
    if argv[1][-5] == '.':
        temp_tif_fn = argv[1][0:-5] + '_temp.tif'
    else:
        temp_tif_fn = argv[1][0:-4] + '_temp.tif'
    cv2.imwrite(temp_tif_fn, arr)
    # add geo info
    geoTrans = src_DS.GetGeoTransform()
    srcPro = src_DS.GetProjection()
    target_ds = gdal.Open(temp_tif_fn,gdal.GA_Update)
    target_ds.SetGeoTransform(geoTrans)
    target_ds.SetProjection(srcPro)
    target_ds.FlushCache()
    target_ds = None
    
    src_DS = gdal.Open(temp_tif_fn)

    srcband = src_DS.GetRasterBand(1)

    #  create output datasource
    # output datasource file name
    if argv[1][-5] == '.':
        dst_layername = argv[1][0:-5]
    else:
        dst_layername = argv[1][0:-4]
    drv = ogr.GetDriverByName("ESRI Shapefile")
    dst_ds = drv.CreateDataSource(dst_layername + ".shp")
    dst_layer = dst_ds.CreateLayer(dst_layername, srs = None )

    gdal.Polygonize( srcband, None, dst_layer, -1, [], callback=None )
    # TODO: add geo info to shpfile
    
    del src_DS
    os.remove(temp_tif_fn)


if __name__ == "__main__":
    main(sys.argv)