複数枚に分割されているオルソ化空中写真を結合して、その後に、2次メッシュで再分割する方法を紹介します。
さらに、その縦×横サイズを10000,10000に正規化する方法を紹介します。なんのために？って思うかもしれませんが、河川水辺の国勢調査のシステムに取り込むためです。。。
1. gdalを使える環境を整えます。
ここからダウンロード http://trac.osgeo.org/osgeo4w/
advanced installから gdalとgdal-pythonをインストールします。
2. 結合します。※以下、OSGeo4Wのコンソールでの作業です。
平面直角座標系にオルソ化された空中写真（GeoTiff）img1.tif img2.tif img3.tifを解像度3mのmerge.tifとして結合します。
gdal_merge -ps 3 3 -o merge.tif img1.tif img2.tif img3.tif
結合するファイルが多い場合は、以下のバッチファイルを使うと便利です。
gdal_wildmerge.bat img*.tif merge.tif "-ps 3 3"
gdal_wildmerge.batという名前で保存してください。
@echo off
:: initial version, 2005-Nov-17, matt wilkie
:: this script is public domain
if [%1]==[] (
echo.
echo -={ gdal_wildmerge }=- Allow gdal_merge to accept wildcards for input files
echo .
echo . gdal_wildmerge.bat [input-file-wildcard] [output-file] {options}
echo .
echo . gdal_wildmerge d:\src\*.tif mosaick.tif "-of GTiff -init -9999"
echo.
goto :EOF
)
:: %1 is the input wildcard mask
:: This is inverted from normal gdal_merge usage in order to allow other
:: options be passed to gdal_merge
:: delayed expansion requires windows NT or later command shell
:: (it won't work on windows 9x)
verify other 2>nul
setlocal ENABLEDELAYEDEXPANSION
if errorlevel 1 goto :wrong_version
:: gdal_merge won't overwrite an existing file
if exist %2 goto :file_exists
:: Build the input file list
set infiles=
for %%a in (%1) do set infiles=!infiles! "%%a"
:: Build the option list
:: set options=
:: for %%b in (%*) do (
:: if %%b != %1 && %%b != %2 (
:: set options=!options! "%%b"
:: )
:: )
:: Here is were the actual work happens
set gdalcmd=gdal_merge %~3 -o %2 %infiles%
echo %gdalcmd%
:: IMPORTANT, if gdal_merge is not call'ed control doesn't return
:: to the batch file
::cd %~dp0
call %gdalcmd%
goto :EOF
:file_exists
echo *** the output file "%2" already exists.
goto :EOF
:wrong_version
echo.
echo Sorry, this batchfile requires DELAYED EXPANSION which is
echo not available in this shell. See
echo http://www.google.com/search?q=enabledelayedexpansion
goto :EOF
3. 平面直角座標系のmerge.tifを緯度経度のmergeLL.tifに変換します。
※epsgコードは、手元のデータに従い変更してくだい。
gdalwarp -r bilinear -s_srs epsg:2452 -t_srs epsg:4612 merge.tif mergeLL.tif
4. 2次メッシュ 574172の範囲で切り取って、tmpA574172.tifに保存します。
gdalwarp.exe -te 141.25 38.5833333333333 141.375 38.6666666666667 mergeLL.tif tmpA574172.tif
5. 正規化座標に変換します。
'gdalinfoでピクセルサイズを読み取ります。 gdalinfo tmpA574172.tif 'コントロールポイントを設定します。gcpは、ピクセルX,ピクセルY,変換後X座標,変換後Y座標,標高値（分からなければ0）の順番です。ピクセルX,ピクセルYは、左上が(0,0)なので注意。 gdal_translate -gcp 0 0 0 10000 0 -gcp 3989 0 10000 10000 0 -gcp 3989 2659 10000 0 0 -gcp 0 2659 0 0 0 tmpA574172.tif tmpB574172.tif 'コントロールポイントに従って、変換します。 gdalwarp -tps tmpB574172.tif 574172.tif
4.〜5.の作業を自動化するには、以下を実行します。
cut.vbs 574172 mergeLL.tif
以下のスクリプトをcut.vbsとして保存してください。
Set objWshShell = WScript.CreateObject("WScript.Shell") Set args = WScript.Arguments '2次メッシュの緯度経度の計算 mesh = args.item(0) ido0 = Left(mesh, 2) / 1.5 + Mid(mesh, 5, 1) * 300 / 3600.0 ido1 = ido0 + 300 / 3600.0 keido0 = Mid(mesh, 3, 2) + 100 + Mid(mesh, 6, 1) * 450 / 3600.0 keido1 = keido0 + 450/3600.0 '2次メッシュの範囲で切り取り strCmdLine = "gdalwarp.exe -te " & keido0 & " " & ido0 & " " & keido1 & " " & ido1 & " " & args.item(1) & " " & "tmpA" & args.item(0) & ".tif" flag = objWshShell.Run(strCmdLine,,True) 'ピクセルサイズをテキストで取得 '%comspec% /c を付けないと、リダイレクトできなかった。/k にすれば、コマンドプロンプトを閉じないので、エラーを確認できる。 strCmdLine = "%comspec% /c gdalinfo.exe" & " " & "tmpA" & args.item(0) & ".tif | findstr /c:""Size is"" " & " > """ & objWshShell.CurrentDirectory & """\tmpA" & args.item(0) & ".txt" flag = objWshShell.Run(strCmdLine,,True) 'テキストから縦、横を取得 Set objFSO = WScript.CreateObject("Scripting.FileSystemObject") Set objFile = objFSO.OpenTextFile("tmpA" & args.item(0) & ".txt") sizeStr = objFile.ReadLine objFile.Close Set objFile = Nothing Set objFSO = Nothing x = Mid(sizeStr,9,InStr(sizeStr,",")-9) y = Mid(sizeStr,InStr(sizeStr,",")+1) '左下(0,0) 右上(10000,10000)の正規化座標となるように、グラントコントロールポイントを設定 strCmdLine = "gdal_translate -gcp 0 0 0 10000 0 -gcp " & x & " 0 10000 10000 0 -gcp " & x & " " & y & " " & "10000 0 0 -gcp 0 " & y & " 0 0 0 " & "tmpA" & args.item(0) & ".tif" & " " & "tmpB" & args.item(0) & ".tif" flag = objWshShell.Run(strCmdLine,,True) 'グランドコントロールポイントに従って変換 strCmdLine = "gdalwarp -tps " & "tmpB" & args.item(0) & ".tif" & " " & args.item(0) & ".tif" flag = objWshShell.Run(strCmdLine,,True)
GDALの使い方が書いてあります。参考下さい。
※自分のためのメモ
- gdalwarp はver1.7以降だと ワイルドカードが使える。ver1.6以前であれば、http://trac.osgeo.org/gdal/wiki/UserDocs/GdalMerge ここのgdal_wildwarp.batを使うと良い。けど、オプションを引数にする部分が間違っているので、↓に修正プログラムを掲載
gdal_mergeでは、ワイルドカードが使えない。- gdalwarpで結合するより、gdal_mergeで結合する方が断然早い。（ファイルが多い場合）
@echo off
:: initial version, 2005-Nov-17, matt wilkie
:: this script is public domain
if [%1]==[] (
echo.
echo -={ gdal_wildwarp }=- Allow gdal_merge to accept wildcards for input files
echo .
echo . gdal_wildwarp.bat [input-file-wildcard] [output-file] {options}
echo .
echo . gdal_wildwarp d:\src\*.tif mosaick.tif -init -9999
echo.
goto :EOF
)
:: %1 is the input wildcard mask
:: This is inverted from normal gdal_merge usage in order to allow other
:: options be passed to gdal_merge
:: delayed expansion requires windows NT or later command shell
:: (it won't work on windows 9x)
verify other 2>nul
setlocal ENABLEDELAYEDEXPANSION
if errorlevel 1 goto :wrong_version
:: gdal_merge won't overwrite an existing file
if exist %2 goto :file_exists
:: Build the input file list
set infiles=
for %%a in (%1) do set infiles=!infiles! "%%a"
:: Build the option list
:: set options=
:: for %%b in (%*) do (
:: if %%b != %1 && %%b != %2 (
:: set options=!options! "%%b"
:: )
:: )
:: Here is were the actual work happens
:: set gdalcmd=gdalwarp -tr 10 10 -r bilinear -s_srs "+proj=tmerc +lat_0=40 +lon_0=140.8333333333333 +k=0.999900 +x_0=0 +y_0=0 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs <>" -t_srs "+proj=longlat +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +no_defs <>" %infiles% %2
set gdalcmd=gdalwarp %infiles% %2 %~3
echo %gdalcmd%
:: IMPORTANT, if gdal_merge is not call'ed control doesn't return
:: to the batch file
::cd %~dp0
call %gdalcmd%
goto :EOF
:file_exists
echo *** the output file "%2" already exists.
goto :EOF
:wrong_version
echo.
echo Sorry, this batchfile requires DELAYED EXPANSION which is
echo not available in this shell. See
echo http://www.google.com/search?q=enabledelayedexpansion
goto :EOF
CADデータの場合の覚書
CADをラスターデータにする方法の覚書
- dxfをcadのフリーソフトで読み込んでpdfに書き出し
- photoshopやgimpでtifに変換。この時に、2階調→白黒にすると軽くなるかも。
- dxfをogrでshpに変換
- shpをillustratorに読み込み。秘密のツールで
- 変換したtifをillustratorに配置して、shpを使って位置合わせ
- 秘密のツールでwldファイルを作成
