Python で GeoTIFF を使ってみる
【概要】
GeoTIFF ファイルの読み込みと,緯度経度の取得を,Python で行う.キーワード: GeoTIFF, Python で GeoTIFF ファイルの読み込み,Python で GeoTIFF ファイルからの緯度経度の取得,osr, gdal, 基盤地図情報, 数値標高モデル, 基盤地図情報ダウンロード
【目次】
前準備
Python 3.12 のインストール
Pythonのインストールを行い、Pythonのプログラムを実行する環境を整える。扱う環境は、Windows搭載パソコンである。金子研究室では、Python 3.12.10を推奨する。
[Windows での Python 3.12 のインストール手順を見るには、ここをクリック]
Windows での Python 3.12 のインストール
以下のいずれかの方法でPython 3.12をインストールする。Pythonがインストール済みの場合、この手順は不要である。
方法 1:winget によるインストール
【インストールコマンドの実行方法】
管理者権限でコマンドプロンプトを起動する(手順:Windowsキーまたはスタートメニュー → cmd と入力 → 右クリック → 「管理者として実行」)。そして、コマンド全体をコマンドプロンプトにコピー&ペーストする。
--scope machine を指定することで、システム全体(全ユーザー向け)にインストールされる。このオプションの実行には管理者権限が必要である。インストール完了後、コマンドプロンプトを再起動するとPATHが反映される。
REM Python 3.12 をシステム領域にインストール
winget install --id Python.Python.3.12 -e --scope machine --silent --accept-source-agreements --accept-package-agreements --override "/quiet InstallAllUsers=1 PrependPath=1 Include_test=0 Include_pip=1 Include_launcher=1 InstallLauncherAllUsers=1 TargetDir=\"C:\Program Files\Python312\""
REM Python と Scripts を PATH 先頭に追加
powershell -NoProfile -Command "$p='C:\Program Files\Python312'; $s=\"$p\Scripts\"; $c=[Environment]::GetEnvironmentVariable('Path','Machine'); if((Test-Path $p) -and (';'+$c+';' -notlike \"*;$p;*\") -and (';'+$c+';' -notlike \"*;$s;*\")){[Environment]::SetEnvironmentVariable('Path',\"$p;$s;$c\",'Machine')}"
方法 2:インストーラーによるインストール
- Python公式サイト(https://www.python.org/downloads/)にアクセスし、「Download Python 3.x.x」ボタンからWindows用インストーラーをダウンロードする。
- ダウンロードしたインストーラーを実行する。
- 初期画面の下部に表示される「Add python.exe to PATH」にチェックを入れてから「Customize installation」を選択する。このチェックを入れ忘れると、コマンドプロンプトから
pythonコマンドを実行できない。 - 「Install Python 3.xx for all users」にチェックを入れ、「Install」をクリックする。
インストールの確認
コマンドプロンプトで以下を実行する。
python --version
バージョン番号(例:Python 3.12.x)が表示されればインストール成功である。「'python' は、内部コマンドまたは外部コマンドとして認識されていません。」と表示される場合は、インストールが正常に完了していない。
Python の開発環境 Visual Studio Code のインストールと Python 用の設定
Python の開発環境Visual Studio Code(プログラムを編集するソフトウェア。以下、VS Code)を整える。
[Windows での Visual Studio Code のインストールと Python 用の設定手順を見るには、ここをクリック]
Windows での Visual Studio Code のインストールと Python 用の設定手順
1. VS Code と拡張機能のインストール
以下のコマンドにより,既存の VS Code を削除し,全ユーザー共有の設定で再インストールしたうえで,拡張機能(VS Code に機能を追加するソフトウェア)をまとめて導入する.
【インストールコマンドの実行方法】
管理者権限でコマンドプロンプトを起動する(手順:Windows キーまたはスタートメニュー → cmd と入力 → 右クリック → 「管理者として実行」)。そして,コマンド全体をコマンドプロンプトにコピー&ペーストする。
インストールコマンド
REM ============================================================
REM Microsoft Visual Studio Code
REM ============================================================
winget uninstall -e --id Microsoft.VisualStudioCode --silent --disable-interactivity --accept-source-agreements
rmdir /s /q C:\ProgramData\vscode-extensions 2>nul
rmdir /s /q "%APPDATA%\Code" 2>nul
rmdir /s /q "%USERPROFILE%\.vscode" 2>nul
rmdir /s /q "%LOCALAPPDATA%\Microsoft\vscode-update" 2>nul
REM VS Code をシステム領域に新規インストール
winget install --scope machine --id Microsoft.VisualStudioCode -e --silent --accept-source-agreements --accept-package-agreements
REM 全ユーザー共有の拡張機能フォルダ
mkdir C:\ProgramData\vscode-extensions 2>nul
icacls "C:\ProgramData\vscode-extensions" /grant "Everyone:(OI)(CI)M" /T
REM スタートメニューのショートカットを --extensions-dir 付きで再作成
rmdir /s /q "C:\ProgramData\Microsoft\Windows\Start Menu\Programs\Visual Studio Code" 2>nul
del "C:\ProgramData\Microsoft\Windows\Start Menu\Programs\Visual Studio Code.lnk" 2>nul
powershell -NoProfile -Command "$s=New-Object -ComObject WScript.Shell; $lnk=$s.CreateShortcut('C:\ProgramData\Microsoft\Windows\Start Menu\Programs\Visual Studio Code.lnk'); $lnk.TargetPath='C:\Program Files\Microsoft VS Code\Code.exe'; $lnk.Arguments='--extensions-dir \"C:\ProgramData\vscode-extensions\"'; $lnk.Save()"
REM ショートカットの検証
powershell -NoProfile -Command "$s=New-Object -ComObject WScript.Shell; $lnk=$s.CreateShortcut('C:\ProgramData\Microsoft\Windows\Start Menu\Programs\Visual Studio Code.lnk'); Write-Host 'TargetPath:' $lnk.TargetPath; Write-Host 'Arguments:' $lnk.Arguments"
REM ファイル / フォルダ右クリックの「Code で開く」を登録
reg add "HKLM\SOFTWARE\Classes\*\shell\VSCode\command" /ve /d "\"C:\Program Files\Microsoft VS Code\Code.exe\" --extensions-dir \"C:\ProgramData\vscode-extensions\" \"%1\"" /f
reg add "HKLM\SOFTWARE\Classes\Directory\shell\VSCode\command" /ve /d "\"C:\Program Files\Microsoft VS Code\Code.exe\" --extensions-dir \"C:\ProgramData\vscode-extensions\" \"%1\"" /f
reg add "HKLM\SOFTWARE\Classes\Directory\Background\shell\VSCode\command" /ve /d "\"C:\Program Files\Microsoft VS Code\Code.exe\" --extensions-dir \"C:\ProgramData\vscode-extensions\" \"%V\"" /f
REM --extensions-dir 付きで起動する code.cmd ラッパを作成
REM (%* を echo で書くと対話的 cmd で失われるため、PowerShell で [char]37+'*' を書き出す)
powershell -NoProfile -Command "$pct=[char]37; $q=[char]34; $c='@echo off'+[char]13+[char]10+$q+'C:\Program Files\Microsoft VS Code\bin\code.cmd'+$q+' --extensions-dir '+$q+'C:\ProgramData\vscode-extensions'+$q+' '+$pct+'*'+[char]13+[char]10; [IO.File]::WriteAllText('C:\ProgramData\vscode-extensions\vscode.cmd',$c,[Text.Encoding]::ASCII)"
REM 拡張機能のインストール
set "CODE=C:\Program Files\Microsoft VS Code\bin\code.cmd"
"%CODE%" --extensions-dir "C:\ProgramData\vscode-extensions" --uninstall-extension GitHub.copilot
"%CODE%" --extensions-dir "C:\ProgramData\vscode-extensions" --uninstall-extension GitHub.copilot-chat
"%CODE%" --extensions-dir "C:\ProgramData\vscode-extensions" --install-extension ms-python.python
"%CODE%" --extensions-dir "C:\ProgramData\vscode-extensions" --install-extension ms-python.vscode-pylance
"%CODE%" --extensions-dir "C:\ProgramData\vscode-extensions" --install-extension ms-python.debugpy
"%CODE%" --extensions-dir "C:\ProgramData\vscode-extensions" --install-extension MS-CEINTL.vscode-language-pack-ja
"%CODE%" --extensions-dir "C:\ProgramData\vscode-extensions" --install-extension saoudrizwan.claude-dev
"%CODE%" --extensions-dir "C:\ProgramData\vscode-extensions" --install-extension rust-lang.rust-analyzer
"%CODE%" --extensions-dir "C:\ProgramData\vscode-extensions" --install-extension tamasfe.even-better-toml
"%CODE%" --extensions-dir "C:\ProgramData\vscode-extensions" --install-extension anthropic.claude-code
"%CODE%" --extensions-dir "C:\ProgramData\vscode-extensions" --install-extension almenon.arepl
"%CODE%" --extensions-dir "C:\ProgramData\vscode-extensions" --list-extensions --show-versions
echo === セットアップ完了 ===
2. Python インタプリタの選択
同一マシンに複数の Python がインストールされている場合,VS Code で使用する Python 本体(インタプリタ:Python プログラムを解釈・実行するソフトウェア)を選択する必要がある.
- コマンドパレット(コマンド名で機能を呼び出す VS Code の入力欄)を開く(
Ctrl+Shift+P) Python: Select Interpreterと入力する
- 表示される一覧から,使用する Python(例:
C:\Program Files\Python312\python.exe)を選択する.
gdal Windows 版のインストール
OSGeo4W のインストールを行っておくこと
GeoTIFF サンプルデータファイルの準備 (1)
- GeoTIFF サンプルデータファイルのダウンロード
次の Web ページから cea.tif をダウンロード
http://download.osgeo.org/geotiff/samples/gdal_eg/
- ダウンロードした .tif ファイルを,分かりやすいディレクトリ(例えばd:\)に保存する.
GeoTIFF サンプルデータファイルの準備 (基盤地図情報の数値標高モデル)
【関連する外部ページ】 http://sanvarie.hatenablog.com/entry/2016/01/10/163027
ありがとうございます.
基盤地図情報・数値標高モデルのダウンロード
- 国土地理院の「基盤地図情報ダウンロードサービス」の Web ページを開く
- 「ログイン画面はこちら」をクリックする.
- 「基盤地図情報・数値標高モデル」の下の
「ファイル選択へ」をクリックする.
- ダウンロードしたい数値標高モデルの範囲を選ぶ.
- 選び終わったら,「ダウンロードファイル確認へ」をクリックする.
- 「すべてチェック」をクリックする.
* 「5A」は航空レーザー測量,「5B」は写真測量である.
- 「まとめてダウンロード」をクリックする.
- ログインする.
- アンケートに協力する(正しく回答する).
- .zip ファイルのダウンロードが始まるので確認する.
- ダウンロードした .zip ファイルを展開(解凍)する.分かりやすいディレクトリに置く.
※ 展開(解凍)すると .zip ファイルができるので確認する.
- さらに展開(解凍)すると .xml ファイルができるので確認する.
ファイル名: FG_GML-5133-41-00-DEM5A-201601001.xml
- 5133: 1次メッシュ番号
- 41: 2次メッシュ番号
- 00: 3次メッシュ番号
- 5A: 「5A」は航空レーザー測量,「5B」は写真測量である.
- 20161001: 年月日
基盤地図情報の数値標高モデルを GeoTIFF に変換
- 「GeoTIFFを格納するフォルダ」を1つ作り,そこに,すべての .xml ファイル を集める.
- http://sanvarie.hatenablog.com/entry/2016/01/10/163027 に記載のプログラムを実行してみる.
優れたソフトウェアの公開に感謝を表明します.
まず,プログラム内のXMLを格納するフォルダ,GeoTIFFを格納するフォルダは適切に設定する必要がある(下図のように).
実行は簡単でした.
* Windows で実行するとき,次のようなエラーが出ることがある.UnicodeDecodeError: 'cp932' codec can't decode byte 0x85 in position 395: illegal multibype sequence.
* 次のように,プログラムファイルを書き換えて回避しました(書き換え箇所2か所).
- Windows でファイルを見てみると,つぎのようになる.
- Windows でファイルを見てみると,つぎのようになる.
GeoTIFF ファイルの読み込み
Python で GeoTIFF のデータを読み込む
ここでは,GeoTIFF のファイル d:/cea.tif を numpy 形式のオブジェクト a に読み込む.確認のため print コマンドで,a の中身, a の要素数, GeoTIFF の縦横, 画素値の最大値と最小値を表示している.
- Python プログラムの実行
from osgeo import gdal import numpy as np ds = gdal.Open('d:/cea.tif', gdal.GA_ReadOnly) a = np.array([ds.GetRasterBand(i + 1).ReadAsArray() for i in range(ds.RasterCount)]) print(a) print(a.shape) print(ds.RasterXSize, ds.RasterYSize)Python プログラムの編集と実行
- 実行結果を確認する.
「514 415」は、GeoTIFF の縦横を表示している.
- 続いて次のプログラムを実行してみる.
「255」や「0」は、画素値の最大値と最小値である.
np.max(a) np.min(a)
- 念のため別の GeoTIFF ファイル E:/FG-GML-5133-41-DEM5A/51334100.tifで,同じことを繰り返してみる.
from osgeo import gdal import numpy as np ds = gdal.Open('E:/FG-GML-5133-41-DEM5A/51334100.tif', gdal.GA_ReadOnly) a = np.array([ds.GetRasterBand(i + 1).ReadAsArray() for i in range(ds.RasterCount)]) print(a) print(a.shape) print(ds.RasterXSize, ds.RasterYSize) print( np.max(a) ) print( np.min(a) )そして、実行結果を確認する.
GeoTIFF ファイルからの緯度経度の取得
- Python で GeoTIFF の緯度経度を取得
https://stackoverflow.com/questions/2922532/obtain-latitude-and-longitude-from-a-geotiff-file に記載のプログラムを次のように書き換えて使用する.(書き換えた部分は太字で示す)。
from osgeo import gdal from osgeo import osr ds = gdal.Open('d:/cea.tif', gdal.GA_ReadOnly) old_cs = osr.SpatialReference() old_cs.ImportFromWkt(ds.GetProjectionRef()) # create the new coordinate system wgs84_wkt = """ GEOGCS["WGS 84", DATUM["WGS_1984", SPHEROID["WGS 84",6378137,298.257223563, AUTHORITY["EPSG","7030"]], AUTHORITY["EPSG","6326"]], PRIMEM["Greenwich",0, AUTHORITY["EPSG","8901"]], UNIT["degree",0.01745329251994328, AUTHORITY["EPSG","9122"]], AUTHORITY["EPSG","4326"]]""" new_cs = osr.SpatialReference() new_cs.ImportFromWkt(wgs84_wkt) # create a transform object to convert between coordinate systems transform = osr.CoordinateTransformation(old_cs, new_cs) # get the point to transform, pixel (0,0) in this case width = ds.RasterXSize height = ds.RasterYSize gt = ds.GetGeoTransform() minx = gt[0] miny = gt[3] + width*gt[4] + height*gt[5] maxx = gt[0] + width*gt[1] + height*gt[2] maxy = gt[3] # get the coordinates in lat long latlong = transform.TransformPoint(minx, miny) print(latlong) latlong = transform.TransformPoint(maxx, maxy) print(latlong)Python プログラムの編集と実行
実行結果の例
- 正しい値なのか確認のため,gdal に付属の gdalinfo コマンドで, .tif ファイルの情報を取得する.先ほどの結果が正しいか確認できる.
太字のところは、実際のディレクトリを調べて読み替えてください.
gdalinfo.exe d:\cea.tif
- 念のため別の GeoTIFF ファイル E:/FG-GML-5133-41-DEM5A/51334100.tifで,同じことを繰り返してみる.
from osgeo import gdal from osgeo import osr ds = gdal.Open('E:/FG-GML-5133-41-DEM5A/51334100.tif', gdal.GA_ReadOnly) old_cs = osr.SpatialReference() old_cs.ImportFromWkt(ds.GetProjectionRef()) # create the new coordinate system wgs84_wkt = """ GEOGCS["WGS 84", DATUM["WGS_1984", SPHEROID["WGS 84",6378137,298.257223563, AUTHORITY["EPSG","7030"]], AUTHORITY["EPSG","6326"]], PRIMEM["Greenwich",0, AUTHORITY["EPSG","8901"]], UNIT["degree",0.01745329251994328, AUTHORITY["EPSG","9122"]], AUTHORITY["EPSG","4326"]]""" new_cs = osr.SpatialReference() new_cs.ImportFromWkt(wgs84_wkt) # create a transform object to convert between coordinate systems transform = osr.CoordinateTransformation(old_cs, new_cs) # get the point to transform, pixel (0,0) in this case width = ds.RasterXSize height = ds.RasterYSize gt = ds.GetGeoTransform() minx = gt[0] miny = gt[3] + width*gt[4] + height*gt[5] maxx = gt[0] + width*gt[1] + height*gt[2] maxy = gt[3] # get the coordinates in lat long latlong = transform.TransformPoint(minx, miny) print(latlong) latlong = transform.TransformPoint(maxx, maxy) print(latlong)