パイソンドライグッズ|リモートセンシング画像ステッチ

序文

「秋初のミルクティー」を飲まなかったので、投稿を更新します。

前回のツイートでは、「Python」と「Landsat」を組み合わせてリモートセンシング画像を作成する方法を紹介しました（[Python干货|リモートセンシング画像を作成]（http://mp.weixin.qq.com/s?__biz=MzUxODUxODI2OQ==&mid=2247483805&idx=1&sn=03bd9134b8db283f0d1e84bfbfc294a8&chksm=f986e50fcef16c19029bf8923326c34c38fa4b1ca61e2c5485750d36096f5118b5ea3c12122d&scene=21#wechat_redirect））。

Landsatデータの場合、特定のエリアの再訪問期間は16日であり、各場所はグローバル参照システム（** WRS ）を使用してインデックス付けされます。つまり、各場所はパス**に対応し、行、隣接する画像間にいくつかの重複領域があります。

Fig.1 World Reference System

一部のリモートセンシング画像アプリケーションのシナリオでは、下の図の翔山港のように、焦点を当てる領域が2つの画像の接合部にある場合、使用する前に画像をつなぎ合わせる必要があります。

一つの画像はこんな感じです。

これは、この記事がマージされた後のケースです。

1. 準備オーケー＃＃

前回の投稿と比較して、リモートセンシング画像のモザイクスプライシングを実現するために、 rasterioと gdalの2つのライブラリを使用しました。

import rasterio as rio
import gdal

まず、2セットのリモートセンシング画像スプライシングを実現するというアイデアを紹介します。まず、2つの隣接する画像を選択してそれぞれの空間範囲を取得し、次に2つのシーンを組み合わせた空間範囲を取得し、「gdal」を使用して新しい「tif」を作成します。ファイル（データ転送用）で、新しく作成した「tif」ファイルで元の2つのシーン画像の開始位置をそれぞれ取得し、対応するデータを新しい「tif」ファイルに書き込んでモザイクスプライシングを実現します。

上記は2シーン画像のスプライシングです。より多くの画像スプライシングであれば適用可能ですが、現在の方法ではより多くの画像をスプライシングする場合、多くのメモリスペースが必要になり、メモリオーバーフローが発生しやすくなります。興味のある友人はそれについて考えることができます。マルチシーン画像のスティッチングを効率的に実現する方法。

2つの重要な処理があります。1つは重複領域の無効な情報を削除する方法であり、もう1つは重複領域のデータを選択する方法です。コードから答えが見つかるといいのですが。

2. 始めましょう##

入力画像の4つのコーナーポイントを取得します。

def tiffileList2filename(tiffileList):
 filename =[]
 prefix =[]for ifile in tiffileList:
  file0 = ifile.split("\\")[-1]
  prefix.append(os.path.join(ifile, file0))
  filename.append(os.path.join(ifile, file0)+"_B1.TIF")return filename, prefix
def get_extent(tiffileList):
 filename, prefix =tiffileList2filename(tiffileList)
 rioData = rio.open(filename[0])
 left = rioData.bounds[0]
 bottom = rioData.bounds[1]
 right = rioData.bounds[2]
 top = rioData.bounds[3]for ifile in filename[1:]:
  rioData = rio.open(ifile)
  left =min(left, rioData.bounds[0])
  bottom =min(bottom, rioData.bounds[1])
  right =max(right, rioData.bounds[2])
  top =max(top, rioData.bounds[3])return left, bottom, right, top, filename, prefix

新しく作成したtifファイルのサイズを取得します。ここでは、 Landsatの空間解像度が30mであることがわかっています。他のリモートセンシングデータの場合は、それに応じて変更する必要があります。

def getRowCol(left, bottom, right, top):
 cols =int((right - left)/30.0)
 rows =int((top - bottom)/30.0)return cols, rows

メインプログラム。plot_rgbは前の記事で使用した関数です。

if __name__ =='__main__':
 tiffileList =[r'PathofLandsat8\LC08_L1TP_118039_20160126_20170330_01_T1',
     r'PathofLandsat8\LC08_L1TP_118040_20160126_20170330_01_T1']
 left, bottom, right, top, filename, prefix =get_extent(tiffileList)
 cols, rows=getRowCol(left, bottom, right, top)
 bands =['B7','B5','B3']
 n_bands =len(bands)
 arr = np.zeros((n_bands, rows, cols), dtype=np.float)
 # tifファイルを開く
 in_ds = gdal.Open(filename[0])for i inrange(len(bands)):
  ibands = bands[i]
  # 新しいtifファイルを作成します
  driver = gdal.GetDriverByName('gtiff')
  out_ds = driver.Create(ibands +'mosaic.tif', cols, rows)
  # tifファイルの投影を設定します
  out_ds.SetProjection(in_ds.GetProjection())
  out_band = out_ds.GetRasterBand(1)
  # 新しいtifファイルの地理的変換を設定します
  gt =list(in_ds.GetGeoTransform())
  gt[0], gt[3]= left, top
  out_ds.SetGeoTransform(gt)
  # ステッチする画像のループ読み取り
  for ifile in prefix:
   in_ds = gdal.Open(ifile +'_'+ ibands +'.tif')
   # 新しく作成されたtifファイルと今回開いたtifファイルの間の座標ドリフトを計算します
   trans = gdal.Transformer(in_ds, out_ds,[])
   # オフセット開始点を取得
   success, xyz = trans.TransformPoint(False,0,0)
   x, y, z =map(int, xyz)
   # バンド情報を読む
   fnBand = in_ds.GetRasterBand(1)
   data = fnBand.ReadAsArray()
   # tifファイルを書き込む前に、最大値を255に設定します。この手順は非常に重要です。
   data = data /65535*255
   data[np.where(data ==255)]=0
   # 画像オーバーラップ部分処理、オーバーラップ部分が最大値を取ります
   xSize = fnBand.XSize
   ySize = fnBand.YSize
   outData = out_band.ReadAsArray(x, y, xSize, ySize)
   data = np.maximum(data, outData)
   out_band.WriteArray(data, x, y)
  del out_band, out_ds
  file2read = ibands +'mosaic.tif'
  arr[i,:,:]= tiff.imread(file2read)
  os.remove(file2read)plot_rgb(arr, rgb=(0,1,2))