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

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

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

代码摘要：

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

不足：

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

附代码

import gdalimport ogrimport sysimport cv2import numpy as npimport osimport osrdef 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)