GDAL修復Landsat ETM+影像條帶

Landsat7 ETM+衛星影像由於衛星感測器故障，導致此後獲取的影像出現了條帶。如下圖所示,影像中均勻的佈滿條帶。

使用GDAL修復影像條帶的程式碼如下：

def gdal_repair(tif_name,out_name,bands):
  """
    tif_name(string): 源影像名
    out_name(string): 輸出影像名
    bands(integer): 影像波段數
  """
  # 開啟影像檔案
  tif = gdal.Open(tif_name)
  
  # 根據檔案型別獲取對應的驅動程式
  driver = gdal.GetDriverByName('GTiff')
  
  # 根據指定檔案的驅動程式，使用現有資料集建立新的可寫資料集
  # 所有支援建立新檔案的驅動程式都支援該`CreateCopy()`方法，   # 但僅`Create()`部分支援該方法
  # CreateCopy的第一個引數為目標檔名，第二個引數為源資料集
  # 第三個引數的值是`0`或`1`，值是`0`。即使無法將原始資料準確地轉換為目標資料，程式仍將執行
  new_img = driver.CreateCopy(out_name,tif,0)
 
  for i in tqdm(range(1,bands)):
    # 分別對每個波段處理
    band = new_img.GetRasterBand(i)
    
    # 使用FillNodata對條帶部分進行插值
    gdal.FillNodata(targetBand = band,maskBand = band,maxSearchDist = 15,smoothingIterations=0)
    
    # 將修復好的波段寫入新資料集中
    new_img.GetRasterBand(i).WriteArray(band.ReadAsArray())

修復之後的效果圖如下所示：

Opencv修復Landsat ETM+影像條帶

使用opencv修復影像的程式碼如下：

def cv2_repair(tif_name):
  # 讀取tif影像
  tif_data = gdal_array.LoadFile(tif_name).astype('float32')

  # 獲取掩膜
  mask = tif_data.sum(axis=0)
  mask = (mask == 0).astype(np.uint8)
  
  bands = tif_data.shape[0]

  res = []
  for i in tqdm(range(bands)):
    # cv.Inpaint(src,inpaintMask,dst,inpaintRadius,flags)
    # src：源影象，可以是8位、16位無符號整型和32位浮點型1通道或者8位無符號3通道
    # inpaintMask：掩膜，8位無符號整型
    # dst：和源影象具有一樣大小的輸出
    # inpaintRadius：演算法考慮的每個已修復點的圓形鄰域的半徑     # flags：修復演算法型別，可選cv2.INPAINT_NS和cv2.INPAINT_TELEA
    
    repaired = cv2.inpaint(tif_data[i],mask,3,flags=cv2.INPAINT_TELEA)
    res.append(repaired)

  return np.array(res)

修復之後的結果圖：

使用opencv修復影像，速度要比Gdal慢許多，但修復質量更好。

Reference

https://www.bogotobogo.com/python/OpenCV_Python/python_opencv3_Image_reconstruction_Inpainting_Interpolation.php

https://gis.stackexchange.com/questions/151020/how-to-use-gdal-fillnodata-in-python

到此這篇關於詳解Python修復遙感影像條帶的兩種方式的文章就介紹到這了,更多相關Python修復遙感影像條帶內容請搜尋我們以前的文章或繼續瀏覽下面的相關文章希望大家以後多多支援我們！