状況

JAXAのGCOM-Cという衛星データを触ります。GCOM-Cという衛星から日本領域のデータをダウンロードしてきてもろもろの前処理を行うときに困っていた問題が最近解決したので解決策を残しておきます。

問題

このような日本領域のグレースケール画像が欲しいのですが、gdal.Warp()で投影変換すると。

このような横長の画像になってしまいます。
gdal_translateではこのような挙動は見られないので、それぞれのデフォルトセッティングに影響がありそう。

解決策

geo referenceから、output Boundsを正しく定義する必要がありました。

例えば日本領域を統合したtif画像のメタデータはこのように表示されます。

Metadata:
  AREA_OR_POINT=Area
Image Structure Metadata:
  COMPRESSION=LZW
  INTERLEAVE=PIXEL
Corner Coordinates:
Upper Left  (10007554.678, 5559752.599) (140d 0'55.68"E, 50d 0' 0.20"N)
Lower Left  (10007554.678, 2223901.040) ( 95d46'34.04"E, 20d 0' 0.08"N)
Upper Right (14455356.757, 5559752.599) (157d45'19.57"W, 50d 0' 0.20"N)
Lower Right (14455356.757, 2223901.040) (138d20'35.83"E, 20d 0' 0.08"N)
Center      (12231455.717, 3891826.819) (134d17' 7.52"E, 35d 0' 0.14"N)

これを見ると、Upper Rightの座標のみ西経が使われていることがわかります。日付変更線を超えてしまっているからです。

gdal.Warpでは投影変換を行うときに、東経180度の続きに西経180度が繋がっていないのだと思います。これを解決するためには、まず東経と西経の表記を統一し、そこから変換対象の画像のminX, minY, maxX, maxYの値を正しく取得する必要があります。

options={"outputBounds": (minX, minY, maxX, maxY)}
gdal.Warp(output_file, input_file, options=gdal.WarpOptions(**options))

ここでは、Lower Leftの座標とUpper Rightの座標がそれぞれのmin maxになります。