本文共 2916 字，大约阅读时间需要 9 分钟。

之前用IDL写高分预处理的时候，就有想过可不可以用python+GDAL写，可是一直卡在了第一步的正射校正，gdal.Warp()函数始终找不到放DEM的位置，最近终于找到了。我尝试了一景1.3G的GF1/WFV，采用ENVI/IDL的脚本运行每次都需要500s以上，而python3+osgeo则稳定在惊人的15s以内！就速度而言，python3+osgeo远远快于ENVI接口。

以下是今天写的简单的代码，包括解压函数，正射校正函数和融合函数（GDAL的融合方法只有默认的加权brovey变换）。运行了一景GF6/PMS，并和IDL版白鸽的结果作了简单对比。

import sys, os, tarfile, tempfile as tmpfrom osgeo import gdal, osrfrom gdalconst import *def unpackage(fn):  dirname, basename = os.path.split(fn)  odir = os.path.join(dirname, 'snowyDove_' + basename[0:-7])  with tarfile.open(fn) as file:    file.extractall(path = odir)  return odirdef ortho(ifn, demfn, res, ofn):  ds = gdal.Open(ifn, GA_ReadOnly)  isNorth = 1 if os.path.basename(ifn).split('_')[3][0] == 'N' else 0  zone = str(int(float(os.path.basename(ifn).split('_')[2][1:])/6) + 31)  zone = int('326' + zone) if isNorth else int('327' + zone)  dstSRS = osr.SpatialReference()  dstSRS.ImportFromEPSG(zone)  tds = gdal.Warp(ofn, ds, format = 'GTiff', xRes = res, yRes = res, dstSRS = dstSRS, rpc = True, transformerOptions = demfn)  ds = tds = Nonedef findImage(dir):  fn = []  with os.scandir(dir) as it:    for entry in it:      if entry.name.endswith('.tiff'):        fn.append(os.path.join(dir, entry.name))  if len(fn) == 1:    return fn[0]  elif len(fn) == 2:    if fn[0].find('PAN') != -1:      return fn[1], fn[0]    else:      return fn[0], fn[1]  elif len(fn) == 3:    return fn  elif len(fn) == 6:    mss = pan = None    for element in fn:      if element.find('MUX.') != -1:        mss = element        continue      if element.find('PAN.') != -1:        pan = element        continue      if mss != None and pan != None:        break    return mss, pandef pansharpen(mss, pan, ofn):  ds = gdal.Open(mss, GA_ReadOnly)  nb = ds.RasterCount  ds = None  vrt = """
     
        
           
      
       %s
            
      
       1
          
     \n""" % pan  for i in range(nb):    vrt += """    
           
      
       %s
            
      
       %s
          
     \n""" % (i+1, mss, i+1)  vrt += """  
    
   \n"""  pansharpends = gdal.Open(vrt)  newds = gdal.GetDriverByName('GTiff').CreateCopy(ofn, pansharpends)  pansharpends = newds = Nonedef snowyDove(argv):  imgdir = unpackage(sys.argv[1])  mss, pan = findImage(imgdir)  odir = os.path.join(os.path.dirname(sys.argv[1]), 'snowyDove_tempFiles')  os.mkdir(odir)  ortho_mss_fn = tmp.mkstemp('.tiff', 'snowyDove_', odir, False)[1]  ortho(mss, sys.argv[2], 8, ortho_mss_fn)  print('mss done')  ortho_pan_fn = tmp.mkstemp('.tiff', 'snowyDove_', odir, False)[1]  ortho(pan, sys.argv[2], 2, ortho_pan_fn)  print('pan done')  sharpen_fn = tmp.mkstemp('.tiff', 'snowyDove_', odir, False)[1]  pansharpen(ortho_mss_fn, ortho_pan_fn, sharpen_fn)  return 0def main():    return snowyDove(sys.argv)if __name__ == '__main__':    sys.exit(snowyDove(sys.argv))

结果对比（上层为gdal结果，下层为IDL结果）：