トップページ -> 実践知識 -> 3次元地図データベース応用 -> Python で GeoTIFF を使ってみる
[サイトマップへ], [サイト内検索へ],

Python で GeoTIFF を使ってみる

サイト構成 連絡先,業績など 実践知識 データの扱い コンピュータ 教材 サポートページ
GeoTIFF ファイルの読み込み, GeoTIFF ファイルからの緯度経書の取得

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


前準備

Anaconda のインストール

前準備として,Python 開発環境のAnaconda のインストールが終わっていること.

Windows での Anaconda のインストール手順は、 別のページで説明している

Ubuntu での Anaconda のインストール手順は、 別のページで説明している

以下,Windows に Anaconda をインストール済みであるものとして説明を続ける.

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

次の手順で,Anaconda の Python 環境に gdal前提パッケージをインストールします.

  1. Anaconda プロンプト (Anaconda Prompt)」を管理者として実行

  2. Anaconda プロンプトで,次のコマンドを実行

    ※ (※ Anaconda や Miniconda を使っていないときは conda がないので pip で代用),

    conda install -c conda-forge gdal
    

    ※ 「Proceed ([y]/n)?」のように表示されたときは y, Enter キー


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 ファイルの情報を取得. 先ほどの結果が正しいか確認できる.
    C:\ProgramData\Anaconda3\pkgs\libgdal-2.2.1-vc14_2\Library\bin\gdalinfo.exe d:\cea.tif