超大影像的矢量裁剪

最近工作遇到了一个小需求，用一个矢量文件去裁剪一个超级大的栅格文件（单文件超过140G）。

这个需求在arcgis很容易实现，但是arcgis有时候会卡顿。综上，用代码的方式去实现这个功能或是更好的选择。

代码

from osgeo import gdal
import os

#  根据shp 裁剪影像
def clip_image_by_shp(input_file, output_file, clip_shape, nodata=0):
    '''
    @todo 根据shp 裁剪影像
    @param shp_file: shp file
    @param img_file: image file
    @param outpath: output path
    @return:
    '''
    outpath = os.path.dirname(output_file)
    if os.path.exists(outpath):
        pass
    else:
        os.makedirs(outpath)


    # 根据已知矩形矢量裁剪图像
    input_file_raster = gdal.Open(input_file)  # 打开待裁剪图像

    input_file_raster.GetProjection()

    # 获取栅格影像的地理参考信息
    geo_transform = input_file_raster.GetGeoTransform()
    projection = input_file_raster.GetProjection()

    # 创建裁剪后的影像文件
    driver = gdal.GetDriverByName('GTiff')
    out_raster = driver.Create(output_file, input_file_raster.RasterXSize, input_file_raster.RasterYSize, input_file_raster.RasterCount, gdal.GDT_Byte)
    out_raster.SetGeoTransform(geo_transform)
    out_raster.SetProjection(projection)


    #shp_ds.GetProjection()
    
    '''
    裁剪影像： 使用 gdal.Warp 函数来进行裁剪操作，
    其中 cutlineDSName 参数指定了Shapefile路径，
    cropToCutline=True 表示按照Shapefile的边界裁剪，
    dstNodata=0 设置了裁剪区域外的像素值为0。'''

    gdal.Warp(output_file, input_file_raster, format='GTiff', cutlineDSName=clip_shape, cropToCutline=True, dstNodata=nodata)

if __name__ == '__main__':
    import time
    a = time.time()
    img_file = r'D:\CGCS2000_202310_Dom.tif'
    shp_file = r'D:\temp\fish.shp'
    outfile = r'D:\temp\\out\\' + os.path.splitext(os.path.basename(img_file))[0] + '.tif'

    clip_image_by_shp(img_file, outfile, shp_file)

    b = time.time()
    print(b-a)

用时2598.5940265655518秒，合计约43分钟。

有时候会感叹，gdal的wrap函数真好强大！