Python で GeoTIFF を使ってみる

GeoTIFF ファイルの読み込み, GeoTIFF ファイルからの緯度経書の取得

キーワード: GeoTIFF, Python で GeoTIFF ファイルの読み込み,Python で GeoTIFF ファイルからの緯度経書の取得,osr, gdal

前準備

Python のインストール(Windows上)

注:既にPython(バージョン3.12を推奨)がインストール済みの場合は,この手順は不要である.

winget(Windowsパッケージマネージャー)を使用してインストールを行う

  1. Windowsで,コマンドプロンプト管理者権限で起動する(例:Windowsキーを押し,「cmd」と入力し,「管理者として実行」を選択)
  2. winget(Windowsパッケージマネージャー)が利用可能か確認する:
    winget --version
    
  3. Pythonのインストール(下のコマンドにより Python 3.12 がインストールされる).
    winget install --scope machine Python.Launcher
    winget install --scope machine Python.Python.3.12
    
  4. 【関連する外部サイト】

    【サイト内の関連ページ】

    gdal パッケージと前提となる Python パッケージのインストール手順

    gdal前提パッケージをインストールする.

    1. Windows では,コマンドプロンプト管理者として実行する.
    2. インストールしたいので、次のコマンドを実行

      * Ubuntu では「sudo python3 setup.py install」

      conda install -c conda-forge gdal
      

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

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

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

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

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

    GeoTIFF ファイルの読み込み

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

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

    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)
    

    実行結果の例

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

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

      https://stackoverflow.com/questions/2922532/obtain-latitude-and-longitude-from-a-geotiff-file に記載のプログラムを次のように書き換えて使用.(書き換えた部分は太字で示す)。

      
      import gdal
      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)
      

      実行結果の例

    2. 確認のためgdal に付属の gdalinfo コマンドで, .tif ファイルの情報を取得. 先ほどの結果が正しいか確認できる.
      gdalinfo.exe d:\cea.tif