gdalでパンシャープンする方法 - 自然環境保全のための周辺技術

landsatの画像をgdal使ってパンシャープンとかの処理ができる便利なプログラムがあります。
https://github.com/gina-alaska/dans-gdal-scripts/wiki

Windowsで使いたかったら自分でコンパイルしてね！とあるのでその方法をメモしておきます。

※プログラムの修正部分で上手くいってない気がするので、現在検証中。

プログラムのダウンロード

githubから取ってきます。
git clone https://github.com/gina-alaska/dans-gdal-scripts.git

visual studio 2010のインストール

下の方にあります。
https://www.visualstudio.com/downloads/download-visual-studio-vs

boostのダウンロード＆ビルド

boostを利用してるので取ってきてビルドします。
http://www.boost.org/users/download/

C:\に展開してvisual studioのコマンドプロンプトで

cd c:\boost_1_58_0
boostrap.bat
.\b2

gdalのダウンロード＆ビルド

http://trac.osgeo.org/gdal/wiki/DownloadSource

展開してvisual studioのコマンドプロンプトで

cd gdal-1.10.1
nmake /f makefile.vc
nmake /f makefile.vc install
nmake /f makefile.vc devinstall

コンパイル

ソースに移動

cd c:\Users\mizutani\Documents\GitHub\dans-gdal-scripts\src

gdal_landsat_pansharpのコンパイル

gdal_landsat_pansharp.ccに下のroundを追加。（vc++にはroundがないらしい。）

double round(double r) {
    return (r > 0.0) ? floor(r + 0.5) : ceil(r - 0.5);
}

コンパイル

cl /EHsc /Ic:\boost_1_58_0 /Ic:\warmerda\bld\include gdal_landsat_pansharp.cc common.cc c:\warmerda\bld
\lib\gdal_i.lib /link /LIBPATH:"c:\boost_1_58_0\stage\lib"

gdal_merge_simpleのコンパイル

131行目を変更

-std::vector<uint8_t> buf(w * chunk_size);
+std::vector<int> buf(w * chunk_size);

コンパイル

cl /EHsc /Ic:\boost_1_58_0 /Ic:\warmerda\bld\include gdal_merge_simple.cc common.cc c:\warmerda\bld
\lib\gdal_i.lib /link /LIBPATH:"c:\boost_1_58_0\stage\lib"

gdal_contrast_stretchのコンパイル

datatype_conversion.ccのstd::isnanを_isnanに変更
gdal_contrast_stretch.ccのstd::isnanを_isnanに変更、std::_isinfを!_finiteに変更
gdal_contrast_stretch.ccに下のroundを追加。（vc++にはroundがないらしい。）

double round(double r) {
    return (r > 0.0) ? floor(r + 0.5) : ceil(r - 0.5);
}

コンパイル

cl /EHsc /Ic:\boost_1_58_0 /Ic:\warmerda\bld\include gdal_contrast_stretch.cc common.cc ndv.cc datatype
_conversion.cc c:\warmerda\bld\lib\gdal_i.lib /link /LIBPATH:"c:\boost_1_58_0\st
age\lib"

lndsat8の画像をダウンロード

https://console.developers.google.com/storage/browser/earthengine-public/landsat/
たとえば　LC81070332014151LGN00

実行

gdal_landsat_pansharp.exe、gdal_contrast_stretch.exe、gdal_merge_simple.exeをLC81070332014151LGN00と同じ場所にコピー

C:\warmerda\bld\bin\の中の　gdal110.dllを、LC81070332014151LGN00と同じ場所にコピー

以下のサイトからスクリプトをprocess_landsatの名前で保存
https://gist.github.com/andymason/6780828

中身を以下のように変更

-gdal_merge_simple
+./gdal_merge_simple
-gdal_landsat_pansharp -rgb ./tmp/rgb.tif -lum ./tmp/rgb.tif 0.25 0.23 0.52 -pan ./tmp/b3-8bit.tif -ndv 0 -o ./tmp/pan.tif
+./gdal_landsat_pansharp -rgb ./tmp/rgb.tif -lum ./tmp/rgb.tif 0.25 0.23 0.52 -pan ./tmp/b8-8bit.tif -ndv 0 -o ./tmp/pan.tif
-gdal_contrast_stretch
+./gdal_contrast_stretch
-convert -verbose -channel B -gamma 0.8 -quality 95 ./tmp/pan.tif final-pan-rgb-corrected.jpg

OSgeo4Wのコンソールで
bash process_landsat LC81070332014151LGN00

その他

他の選択肢も色々あります。結局どれが良いのかそのうち試します。

orfeoツールボックスでもできます。なんか遅いような。
otbcli_BundleToPerfectSensor

geosageを使えばできますが、オープンソースじゃないので。。。
http://www.geosage.com/highview/features_landsat8.html

grassでもできますが、なんか面倒。
http://grass.osgeo.org/grass70/manuals/i.pansharpen.html

landsat-utilでもできます。
https://github.com/developmentseed/landsat-util