Python で GeoTIFF を使ってみる
GeoTIFF ファイルの読み込み,
GeoTIFF ファイルからの緯度経書の取得
キーワード: GeoTIFF, Python で GeoTIFF ファイルの読み込み,Python で GeoTIFF ファイルからの緯度経書の取得,osr, gdal
前準備
Python のインストール(Windows上)
注:既にPython(バージョン3.12を推奨)がインストール済みの場合は,この手順は不要である.
winget(Windowsパッケージマネージャー)を使用してインストールを行う
- Windowsで,コマンドプロンプトを管理者権限で起動する(例:Windowsキーを押し,「cmd」と入力し,「管理者として実行」を選択)
- winget(Windowsパッケージマネージャー)が利用可能か確認する:
winget --version
- Pythonのインストール(下のコマンドにより Python 3.12 がインストールされる).
- Python詳細ガイド:Pythonまとめ »
- Windows では,コマンドプロンプトを管理者として実行する.
- インストールしたいので、次のコマンドを実行
* Ubuntu では「sudo python3 setup.py install」
conda install -c conda-forge gdal
- GeoTIFF サンプルデータファイルのダウンロード
次の Web ページから cea.tif をダウンロード
http://download.osgeo.org/geotiff/samples/gdal_eg/
- ダウンロードした .tif ファイルを,分かりやすいディレクトリ(例えばd:\)に保存
- 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)
実行結果の例
- 確認のためgdal に付属の gdalinfo コマンドで, .tif ファイルの情報を取得. 先ほどの結果が正しいか確認できる.
gdalinfo.exe d:\cea.tif
【関連する外部サイト】
【サイト内の関連ページ】
gdal パッケージと前提となる Python パッケージのインストール手順
gdal前提パッケージをインストールする.
GeoTIFF サンプルデータファイルの準備
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)

実行結果の例
