虽然python和GDAL配合可以很容易实现影像剪切，据我了解其对于剪切后输出范围有三个：按照原影像范围输出、按照剪切数据（如一个polygon shapefile)的范围输出和按照自定义的范围输出。这看来是比较完美的，能够满足绝大多数对输出范围的要求。但我在实际工作中遇到一种情况，貌似这三种方式不能直接或者完美的解决。先举个例子（图1）：我要剪切出某一现状地物（如河流）的影像，但是这个现状地物很长，超过多幅影像的覆盖范围，另如果要处理的影像非常多，这个时候怎么办。

                              图1. 剪切文件和被剪切影像的叠加示意。

如图1所示，如果我们直接去剪切：若按照影像范围输出，那么输出文件的大小同影像的大小（占据的磁盘空间），当影像非常多的时候，对电脑的硬件要求就很高，而且不利于后续进一步处理；若按照剪切文件的范围输出的话，那么输出影像更大，更不利后续处理；若按照自定义范围去输出的话，问题就是如何自动获取overlap部分的范围呢？

获取利用arcpy可以方便的解决，我也想啊。问题是我的python是64位的，而arcpy是32位的，不能调用，也不能重装符合要求的版本（各种原因），总之只能用64位python去实现。下面是我折腾出来的解决办法，不一定高效，但能用。

思路是分三步解决：

1、首先按照影像范围输出一个剪切文件

2、只提取该剪切文件中有效值部分，无效值(nodata，白色部分)删除，形成掩膜文件

3、将掩膜文件转换成polygon。

# -*- coding: utf-8 -*-
"""
Created on Sun Mar 25 14:21:46 2018
@author: wenzf
"""
import subprocess
from osgeo import ogr
def clip_rst(inraster, inshape, outraster, crop = 0):
    if crop == 0:
        cmd = ['gdalwarp', '-dstnodata', '0', inraster, outraster, '-cutline', inshape, '-overwrite']
    else:
        cmd = ['gdalwarp', '-dstnodata', '0', '-crop_to_cutline', inraster, outraster, '-cutline', inshape, '-overwrite']
    subprocess.call(cmd)
    
#   or:    
#    from osgeo import gdal
#    gdal.Warp(outraster,
#              inraster,
#              cutlineDSName = inshape,
#              dstNodata = 0)
    
def reclass_rst(inraster, outraster):
    exe = 'C:\Users\wenzf\Anaconda2\Lib\site-packages\osgeo\scripts\gdal_calc.py'
    cmd = ['python', exe, '-A', inraster, '--outfile='+outraster, '--calc=(A>0)', \
           '--NoDataValue=0','--overwrite']
    subprocess.call(cmd)
    
def rst2polyg(inraster, outshape):
    driver = ogr.GetDriverByName("ESRI Shapefile")
    driver.DeleteDataSource(outshape)
    exe = 'C:\Users\wenzf\Anaconda2\Scripts\gdal_polygonize.py'
    cmd = ['python', exe, inraster, '-f', 'ESRI Shapefile', outshape]
    subprocess.call(cmd)

#---- clipping raster -----
inshape = folder+'river_buffer.shp'
# ----------生成最小剪切范围 --------------
tmp = fns[0]    # fns是包含需要进行剪切的影像文件路径的数组
inraster = tmp
outraster = inraster.replace('.TIF', '_clip.TIF')
core.clip_rst(inraster, inshape, outraster)
outr_tmp = inraster.replace('.TIF', '_tmp.TIF')
core.reclass_rst(outraster, outr_tmp)
driver = gdal.GetDriverByName('GTiff')
driver.Delete(outraster)
tmp_shape = folder+f_name.replace('.tar.gz','.shp')
core.rst2polyg(outr_tmp, tmp_shape)
driver.Delete(outr_tmp)
# ----------生成最小剪切范围------------
# ----------根据最小范围剪切----------
outfns = []
for fn in fns:
    inraster = fn
    outraster = inraster.replace('.TIF', '_clip.TIF')
    outfns = np.append(outfns, outraster)
    core.clip_rst(inraster, tmp_shape, outraster, 1)

补充说明：

上面的方法也有一些缺陷，就是当得到的新shapefile存在多个feature时，就不能成功运行。对此，有两种方案来解决：一是直接获得新生成的shapefile文件的extent，然后将该extent输入gdalwarp参数中进行剪切；一是直接绕过前面所列三个步骤中的二、三步，即不通过生成新shapefile。在获得初始剪切的raster后，直接获得该raster有效值部分的最小extent，然后将该extent输入gdalwarp参数中进行剪切。

def rst_extent(inraster, nodata):
    inds = gdal.Open(inraster)
    trans = inds.GetGeoTransform()
    upl_x = trans[0]
    upl_y = trans[3]
    rsx = trans[1]
    rsy = trans[5]
    img = inds.ReadAsArray()
    msk = np.where(img != nodata)
    row_min = msk[0].min()
    col_min = msk[1].min()
    row_max = msk[0].max()
    col_max = msk[1].max()
    x_min = upl_x + (col_min+1)*rsx
    y_max = upl_y + (row_min+1)*rsy
    x_max = upl_x + (col_max+1)*rsx
    y_min = upl_y + (row_max+1)*rsy
    
    inds = img = None    
    return (x_min, y_min, x_max, y_max)
#---- clipping raster -----
# -----生成最小剪切范围 -----
inshape = folder+'river_buffer.shp'
tmp = fns[0]
inraster = tmp
outraster = inraster.replace('.TIF', '_clip.TIF')
cmd = ['gdalwarp', '-dstnodata', '0', inraster, outraster, '-cutline', inshape, '-overwrite']
subprocess.call(cmd)
min_extent = rst_extent(outraster, 0)
driver = gdal.GetDriverByName('GTiff')
driver.Delete(outraster)
# ----------生成最小剪切范围------------
# ----------根据最小范围剪切----------
outfns = []
extent = '-te '+ str(min_extent[0]) +' '+str(min_extent[1]) +' '+str(min_extent[2]) +' '+str(min_extent[3])
for fn in fns:
    inraster = fn
    outraster = inraster.replace('.TIF', '_clip.TIF')
    outfns = np.append(outfns, outraster)
    cmd = ['gdalwarp', '-dstnodata', '0', '-te', str(int(min_extent[0])), str(int(min_extent[1])), \
           str(int(min_extent[2])), str(int(min_extent[3])), '-cutline', inshape, inraster, outraster, '-overwrite']
    subprocess.call(cmd)