python3+osgeo处理高分影像初探

之前用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 tmp
from osgeo import gdal, osr
from 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 odir

def 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 = None

def 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, pan

def pansharpen(mss, pan, ofn):
  ds = gdal.Open(mss, GA_ReadOnly)
  nb = ds.RasterCount
  ds = None
  vrt = """<VRTDataset subClass="VRTPansharpenedDataset">
  <PansharpeningOptions>
    <PanchroBand>
      <SourceFilename relativeToVRT="0">%s</SourceFilename>
      <SourceBand>1</SourceBand>
    </PanchroBand>\n""" % pan
  for i in range(nb):
    vrt += """    <SpectralBand dstBand="%s">
      <SourceFilename relativeToVRT="0">%s</SourceFilename>
      <SourceBand>%s</SourceBand>
    </SpectralBand>\n""" % (i+1, mss, i+1)
  vrt += """  </PansharpeningOptions>
</VRTDataset>\n"""
  pansharpends = gdal.Open(vrt)
  newds = gdal.GetDriverByName('GTiff').CreateCopy(ofn, pansharpends)
  pansharpends = newds = None

def 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 0

def main():
    return snowyDove(sys.argv)

if __name__ == '__main__':
    sys.exit(snowyDove(sys.argv))

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

后来分析才发现，其实产生这种错位差异的原因可能在于正射校正环节：GDAL的正射校正保持了原汁原味的锯齿状，而ENVI的正射校正则使用了一定的平滑，所以才可能导致全色和多光谱的错位。