Python で GeoTIFF を使ってみる

【概要】

GeoTIFF ファイルの読み込みと,緯度経度の取得を,Python で行う.

キーワード: GeoTIFF, Python で GeoTIFF ファイルの読み込み,Python で GeoTIFF ファイルからの緯度経度の取得,osr, gdal, 基盤地図情報, 数値標高モデル, 基盤地図情報ダウンロードサービス

【目次】

前準備

Python のインストール

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\""
if not "%ERRORLEVEL%"=="0" ( color 0c & echo Python 3.12 のインストールに失敗しました & ping 127.0.0.1 -n 6 >nul & color )

REM Python と Scripts を PATH 先頭に追加
powershell -NoProfile -Command "$p='C:\Program Files\Python312'; $s=\"$p\Scripts\"; if(Test-Path $p){$k=[Microsoft.Win32.Registry]::LocalMachine.OpenSubKey('SYSTEM\CurrentControlSet\Control\Session Manager\Environment',$true); $c=$k.GetValue('Path','',[Microsoft.Win32.RegistryValueOptions]::DoNotExpandEnvironmentNames); $t=$k.GetValueKind('Path'); $new=$c; if((';'+$new+';') -notlike \"*;$p;*\"){$new=$p+';'+$new}; if((';'+$new+';') -notlike \"*;$s;*\"){$new=$s+';'+$new}; if($new -ne $c){$k.SetValue('Path',$new,$t)}; $k.Close()}"

REM 現在のセッションにも反映(システムPATHを再取得して連結)
for /f "usebackq tokens=2,*" %A in (`reg query "HKLM\SYSTEM\CurrentControlSet\Control\Session Manager\Environment" /v Path`) do set "PATH=%B"

REM pip / wheel の更新
python -m pip install --no-user -U pip wheel
if not "%ERRORLEVEL%"=="0" ( color 0c & echo pip / wheel の更新に失敗しました & ping 127.0.0.1 -n 6 >nul & color )

方法 2:インストーラーによるインストール

  1. Python公式サイト(https://www.python.org/downloads/)にアクセスし、「Download Python 3.x.x」ボタンからWindows用インストーラーをダウンロードする。
  2. ダウンロードしたインストーラーを実行する。
  3. 初期画面の下部に表示される「Add python.exe to PATH」にチェックを入れてから「Customize installation」を選択する。このチェックを入れ忘れると、コマンドプロンプトから python コマンドを実行できない。
  4. 「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 ============================================================
REM Build Tools + Desktop development with C++(VCTools)+ 追加コンポーネント(一括)
REM 未インストール時: winget で新規インストール
REM インストール済み時: setup.exe modify でコンポーネント追加(バージョンは変更しない)
winget list --id Microsoft.VisualStudio.BuildTools 2>nul | findstr /i "BuildTools" >nul 2>&1
if %ERRORLEVEL% EQU 0 (
    for /f "usebackq delims=" %P in (`"C:\Program Files (x86)\Microsoft Visual Studio\Installer\vswhere.exe" -products Microsoft.VisualStudio.Product.BuildTools -property installationPath`) do start /wait "" "C:\Program Files (x86)\Microsoft Visual Studio\Installer\setup.exe" modify --installPath "%P" --add Microsoft.VisualStudio.Workload.VCTools --add Microsoft.VisualStudio.Workload.MSBuildTools --add Microsoft.VisualStudio.Component.VC.CMake.Project --add Microsoft.VisualStudio.Component.VC.Llvm.Clang --add Microsoft.VisualStudio.Component.VC.Llvm.ClangToolset --add Microsoft.VisualStudio.Component.Windows11SDK.26100 --add Microsoft.VisualStudio.Component.VC.v143.x86.x64 --includeRecommended --quiet --norestart --nocache
    if not "%ERRORLEVEL%"=="0" ( color 0c & echo Build Tools のコンポーネント追加に失敗しました & ping 127.0.0.1 -n 6 >nul & color )
) else (
    winget install --scope machine --id Microsoft.VisualStudio.BuildTools -e --silent --disable-interactivity --force --accept-source-agreements --accept-package-agreements --override "--quiet --wait --norestart --nocache --add Microsoft.VisualStudio.Workload.VCTools --includeRecommended --add Microsoft.VisualStudio.Workload.MSBuildTools --add Microsoft.VisualStudio.Component.VC.CMake.Project --add Microsoft.VisualStudio.Component.VC.Llvm.Clang --add Microsoft.VisualStudio.Component.VC.Llvm.ClangToolset --add Microsoft.VisualStudio.Component.Windows11SDK.26100 --add Microsoft.VisualStudio.Component.VC.v143.x86.x64"
    if not "%ERRORLEVEL%"=="0" ( color 0c & echo Build Tools のインストールに失敗しました & ping 127.0.0.1 -n 6 >nul & color )
)

REM 全ユーザー共有の拡張機能フォルダ
if not exist "C:\ProgramData\vscode-extensions" mkdir "C:\ProgramData\vscode-extensions"
icacls "C:\ProgramData\vscode-extensions" /grant "Everyone:(OI)(CI)M" /T

REM スタートメニューのショートカットを --extensions-dir 付きで再作成
if exist "C:\ProgramData\Microsoft\Windows\Start Menu\Programs\Visual Studio Code" rmdir /s /q "C:\ProgramData\Microsoft\Windows\Start Menu\Programs\Visual Studio Code"
if exist "C:\ProgramData\Microsoft\Windows\Start Menu\Programs\Visual Studio Code.lnk" del "C:\ProgramData\Microsoft\Windows\Start Menu\Programs\Visual Studio Code.lnk"
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

REM settings.json を作成(自動更新オフ、Python、Claude Code 設定)
if not exist "%APPDATA%\Code\User" mkdir "%APPDATA%\Code\User"
python -c "import json,os;data={'update.mode':'none','update.enableWindowsBackgroundUpdates':False,'extensions.autoUpdate':False,'python.defaultInterpreterPath':r'C:\Program Files\Python312\python.exe','claudeCode.environmentVariables':[{'name':'ANTHROPIC_API_KEY','value':'not-needed'},{'name':'ANTHROPIC_AUTH_TOKEN','value':'ollama'},{'name':'ANTHROPIC_BASE_URL','value':'http://localhost:11434'},{'name':'ANTHROPIC_MODEL','value':'glm-4.7-flash'},{'name':'CLAUDE_CODE_DISABLE_NONESSENTIAL_TRAFFIC','value':'1'}]};p=os.path.join(os.environ['APPDATA'],'Code','User','settings.json');open(p,'w',encoding='utf-8').write(json.dumps(data,indent=4));print('Done:',p)"

REM 自動更新の抑止ポリシー(settings.json に加えて、レジストリ側でも明示的にオフ)
reg add "HKLM\SOFTWARE\Policies\Microsoft\VSCode" /v "UpdateMode" /t REG_SZ /d "none" /f
echo === セットアップ完了 ===

2. Python インタプリタの選択

同一マシンに複数の Python がインストールされている場合,VS Code で使用する Python 本体(インタプリタ:Python プログラムを解釈・実行するソフトウェア)を選択する必要がある.

  1. コマンドパレット(コマンド名で機能を呼び出す VS Code の入力欄)を開く(Ctrl+Shift+P)
  2. Python: Select Interpreter と入力する
  3. 表示される一覧から,使用する Python(例:C:\Program Files\Python312\python.exe)を選択する.

gdal Windows 版のインストール

OSGeo4W のインストールを行っておくこと

GeoTIFF サンプルデータファイルの準備 (1)

  1. GeoTIFF サンプルデータファイルのダウンロード

    次の Web ページから cea.tif をダウンロード

    http://download.osgeo.org/geotiff/samples/gdal_eg/

  2. ダウンロードした .tif ファイルを,分かりやすいディレクトリ(例えばd:\)に保存する.

GeoTIFF サンプルデータファイルの準備 (基盤地図情報の数値標高モデル)

基盤地図情報・数値標高モデルのダウンロード

  1. 国土地理院の「基盤地図情報ダウンロードサービス」の Web ページを開く

    https://service.gsi.go.jp/kiban/

  2. 画面右上の「ログイン」をクリックし,「新規登録」から利用者登録(無料)を行い,ID とパスワードを取得してログインする.

    * ログイン ID の有効期限は 1 年である.

  3. 「基盤地図情報「基本項目」・「数値標高モデル」のダウンロード」で,「数値標高モデル」を選ぶ.
  4. ダウンロードしたい数値標高モデルの範囲(メッシュ)を選び,検索を行う.
  5. 検索結果から,ダウンロードするデータを選ぶ.

    データの種別は次の通りである.

    • DEM1A(1m メッシュ), DEM5A(5m メッシュ): 航空レーザ測量
    • DEM5B, DEM5C(5m メッシュ): 写真測量
    • DEM10B(10m メッシュ): 地形図の等高線
  6. ダウンロードを行う.

    「まとめてダウンロード」を選んだ場合は,利用者登録したメールアドレスにダウンロード用の URL が送られる.1 ファイルずつ「ダウンロード」を選んだ場合は,Web ブラウザでダウンロードが始まる.

  7. ダウンロードした .zip ファイルを展開(解凍)する.分かりやすいディレクトリに置く.

    ※ 展開(解凍)すると .zip ファイルができるので確認する.

  8. さらに展開(解凍)すると .xml ファイルができるので確認する.

    ファイル名の例: FG-GML-5133-41-00-DEM5A-20161001.xml

    • 5133: 1次メッシュ番号
    • 41: 2次メッシュ番号
    • 00: 3次メッシュ番号
    • DEM5A: データの種別
    • 20161001: 作成年月日
  9. データの座標参照系を確認する.

    令和 7 年 7 月 31 日以降に提供された基盤地図情報は,標高成果の改定に伴い,座標参照系が JGD2011 から JGD2024 に変更されている.

基盤地図情報の数値標高モデルを GeoTIFF に変換

数値標高モデルの .xml ファイルを GeoTIFF に変換するには,QGIS のプラグイン QuickDEM4JP を用いる方法がある.

https://plugins.qgis.org/plugins/QuickDEM4JP/

  1. QGIS をインストールする.

    OSGeo4W のインストールで QGIS をインストールできる.

  2. QGIS を起動し,メニューの「プラグイン」で「プラグインの管理とインストール」を選ぶ.
  3. 「QuickDEM4JP」を検索し,インストールする.
  4. メニューの「プラグイン」で「QuickDEM4JP」を選ぶ.
  5. ダウンロードした .zip ファイル,または .xml ファイルを指定する.
  6. 出力形式として「GeoTIFF」を選び,実行する.

    変換された GeoTIFF ファイルが作られ,QGIS に読み込まれる.

* 別のツールを用いる方法は 基盤地図情報標高DEMデータ変換ツール DEMTOOL の紹介 で説明している.

GeoTIFF ファイルの読み込み

Python で GeoTIFF のデータを読み込む

ここでは,GeoTIFF のファイル d:/cea.tif を numpy 形式のオブジェクト a に読み込む.確認のため print コマンドで,a の中身, a の要素数, GeoTIFF の縦横, 画素値の最大値と最小値を表示する.

  1. Python プログラムの実行
    from osgeo import gdal
    import numpy as np
    
    gdal.UseExceptions()
    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 プログラムの編集と実行

  2. 実行結果を確認する.

    「514 415」は,GeoTIFF の横と縦の画素数である.

  3. 続いて次のプログラムを実行する.

    「255」や「0」は,画素値の最大値と最小値である.

    np.max(a)
    np.min(a)
    
  4. 数値標高モデルから変換した別の GeoTIFF ファイル d:/dem/5133-41-00.tifで,同じことを繰り返す.

    太字のところは,実際のファイル名に読み替える.

    from osgeo import gdal
    import numpy as np
    
    gdal.UseExceptions()
    ds = gdal.Open('d:/dem/5133-41-00.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 では,画素値は標高 [m]である.データが無い画素には,-9999 などの値が入る.

GeoTIFF ファイルからの緯度経度の取得

  1. Python で GeoTIFF の緯度経度を取得

    次のプログラムは,GeoTIFF の座標系から緯度経度 (EPSG:4326) への座標変換を行い,左下と右上の座標を表示する.

    * gdal 3 以降では,座標の順序が座標系の定義に従う.SetAxisMappingStrategy で OAMS_TRADITIONAL_GIS_ORDER を指定すると,「経度, 緯度」の順序になる.

    from osgeo import gdal
    from osgeo import osr
    
    gdal.UseExceptions()
    ds = gdal.Open('d:/cea.tif', gdal.GA_ReadOnly)
    
    # 変換元の座標系(GeoTIFF ファイルの座標系)
    old_cs = osr.SpatialReference()
    old_cs.ImportFromWkt(ds.GetProjectionRef())
    old_cs.SetAxisMappingStrategy(osr.OAMS_TRADITIONAL_GIS_ORDER)
    
    # 変換先の座標系(WGS 84 の緯度経度)
    new_cs = osr.SpatialReference()
    new_cs.ImportFromEPSG(4326)
    new_cs.SetAxisMappingStrategy(osr.OAMS_TRADITIONAL_GIS_ORDER)
    
    # 座標変換のオブジェクト
    transform = osr.CoordinateTransformation(old_cs, new_cs)
    
    # 左下と右上の座標を求める
    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]
    
    # 経度, 緯度, 高さ の順に表示される
    print(transform.TransformPoint(minx, miny))
    print(transform.TransformPoint(maxx, maxy))
    

    Python プログラムの編集と実行

    実行結果の例

  2. 結果を確認するために,gdal に付属の gdalinfo コマンドで .tif ファイルの情報を取得し,四隅の緯度経度を比べる.

    太字のところは,実際のディレクトリに読み替える.

    gdalinfo.exe d:\cea.tif
  3. 数値標高モデルから変換した別の GeoTIFF ファイル d:/dem/5133-41-00.tifで,同じことを繰り返す.

    太字のところは,実際のファイル名に読み替える.

    from osgeo import gdal
    from osgeo import osr
    
    gdal.UseExceptions()
    ds = gdal.Open('d:/dem/5133-41-00.tif', gdal.GA_ReadOnly)
    
    # 変換元の座標系(GeoTIFF ファイルの座標系)
    old_cs = osr.SpatialReference()
    old_cs.ImportFromWkt(ds.GetProjectionRef())
    old_cs.SetAxisMappingStrategy(osr.OAMS_TRADITIONAL_GIS_ORDER)
    
    # 変換先の座標系(WGS 84 の緯度経度)
    new_cs = osr.SpatialReference()
    new_cs.ImportFromEPSG(4326)
    new_cs.SetAxisMappingStrategy(osr.OAMS_TRADITIONAL_GIS_ORDER)
    
    # 座標変換のオブジェクト
    transform = osr.CoordinateTransformation(old_cs, new_cs)
    
    # 左下と右上の座標を求める
    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]
    
    # 経度, 緯度, 高さ の順に表示される
    print(transform.TransformPoint(minx, miny))
    print(transform.TransformPoint(maxx, maxy))