サイトの全体構成

Python で GeoTIFF を使ってみる

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

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

前準備

Python のインストール（Windows 上）

Windows での Python 3.10，関連パッケージ，Python 開発環境のインストール: 別ページ »で説明

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

Python のまとめ: 別ページ »にまとめ

【関連する外部ページ】

Python の公式ページ: https://www.python.org/

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

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

Windows では，コマンドプロンプトを管理者として実行する．
インストールしたいので、次のコマンドを実行
※ 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
```
本サイトは金子邦彦研究室のWebページである．
資料等の公開では，原則，「クリエイティブコモンズ BY NC SA」として公開するようにしている． PDFファイル，パワーポイントファイルなどには，「クリエイティブコモンズ BY NC SA」を明記するとともに，ロゴを記載するようにしている（作業が間に合っていない分もあるのでご容赦ください）．
公開している資料をご利用になる場合の，再配布の条件，剽窃の防止などについて，別ページ »で説明再配布や資料改変の際には，そのページをご確認ください．

サイトマップは，サイトマップのページをご覧下さい．本サイト内の検索は，サイト内検索のページをご利用下さい．

問い合わせ先：金子邦彦（かねこくにひこ）

Python で GeoTIFF を使ってみる

前準備

Python のインストール（Windows 上）

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

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

GeoTIFF ファイルの読み込み

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