国土数値情報の土地利用細分メッシュを広範囲に結合してラスタデータに変換する方法を紹介します。
0.OSGeo4Wの準備
OSGeo4Wを↓からとってきてインストールします。
http://trac.osgeo.org/osgeo4w/wiki/OSGeo4W_jp
advanced install でgdal-dev（バージョン1.8.0以上）をインストールします。
OSGeo4Wのコンソールを起動して、gdaldevと入力して、gdal（バージョン1.8.0以上）を利用できるようにします。
以下、コマンドの実行はOSGeo4Wコンソールでの作業になります。
1.ダウンロード
土地利用細分メッシュを↓からダウンロードします。
http://nlftp.mlit.go.jp/ksj/jpgis/datalist/KsjTmplt-L03-b.html
2.shpファイルに変換
国土数値情報JPGISをshpファイルに変換するには↓こちらを利用します。
http://nlftp.mlit.go.jp/ksj/jpgis/jpgis_tool.html
デフォルトの設定ではメモリが足りなくてエラーになるので
C:\Program Files\KsjTool\GUI.batの中のメモリ使用量を↓のように増やします。
set option=-Xms1024m -Xmx1024m
3.ラスタに変換
↓のコマンドをshp2rast.batで保存します。
for %%A in (%1) do gdal_rasterize -ot UInt16 -a 土地利用種 -tr 0.00125 0.00083 -l %%~nA %%A %%~nA.tif
※↑ -tr の値を修正しました。（2012/12/3）
2.のshpファイルがある場所で、バッチファイルを実行して解像度0.001秒（約100m）のtifを作成します。
shp2rast.bat *.shp
4.ラスタを結合
以下のコマンドで3.で変換したtifをout.tifに結合します。
gdal_merge -o out.tif *tif
5.ラスタをUTMに変換
以下のコマンドで、緯度経度のラスタをUTM54に変換します。
gdalwarp -s_srs epsg:4612 -t_srs epsg:3100 out.tif outUTM.tif
以上です。
補足
今回、最終的にUTMのラスタに変換したかったのですが、その方法は以下の３パターン考えられました。
1. UTM変換→ラスタ変換→ラスタ結合
2. shp結合→ラスタ変換→UTM変換
3. ラスタ変換→ラスタ結合→UTM変換
しかし、実際には1.、2.は上手くいきませんでした。理由は以下のとおりです。
1.は、ダウンロードしたデータ範囲に重なり部分が無いので、UTM変換→ラスタ→結合の順番だとデータの境界部分で若干のすき間ができてしまいます。
2.は、shp結合→ラスタ変換の部分で、プログラムエラーがでてしまいます。shpの容量が大きすぎて、gdal_rasterizeが処理できないためだと思います。
1.2.を試した際のバッチファイルを一応、メモっときます。
jgd2utm54.bat（shpファイルをutm54のshpに変換。outputフォルダをあらかじめ作っておく。）
for %%A in (%1) do ogr2ogr -s_srs epsg:4612 -t_srs epsg:3100 UTM%%A output\%%A
appendshp.bat（out.shpにshpファイルをくっつける。あらかじめogr2ogrでout.shpを作っておく）
for %%A in (%1) do ogr2ogr -update -append -nln out out.shp %%A
GDALの使い方はこちら↓に載っています。
座標系の詳しいことはこちら↓に載っていします。
