PythonでGDALを使ってみる。[Chapter 3]
PythonでOGRをちょろっと使ってみる。[Chapter 2] に引き続き、GDAL/OGR Libraryについての基本を学んでおきます。PythonでGeoSpatialをやってみる。[Chapter 1] でPythonとGeoSpatialに関するlibraryが設定できているもとして進めていきます。
本tutorialは、Python Geospatial Developmentを参考にして、実践的に使えるPythonでGeoSpatialの使い方を学んでいくものです。
また、本Blog中のsource code等に関しては、あくまでも参考としてください。なにがあっても責任とれませんので。
そこんところ、ヨロシク~~!!
GDALとは?
OGR-GroPacific.orgに日本語の情報がありますので、参考にしてください。
本家のHPは、GDAL/OGR HPです。
GDALとは、「ラスター空間情報データを扱うためのlibrary群とApplication群」と言えるでしょうか。
Python Geospatial Developmentによれば、
の図のようなobject構成になっています。
gdalで扱うことのできるformatは、gdal formatを参考にしてください。網羅している数がすごいですね。
GDALのexample code
まず、sample dataのsample rasterをdownloadして、その後解凍しておきます。
私の開発環境では、f:¥raster_dataのfolderを作成して、解凍したfileをcopyしておきます。
Ediorで下記のcodeを入力して、保存します。
保存先はf:¥python_example¥chapter_3¥ folderにgdal-example.pyで保存します。
[python]
from osgeo import gdal, gdalconst
fn = ‘f:¥¥raster_data¥¥example.tif’
ds = gdal.Open(fn, gdalconst.GA_ReadOnly)
if ds is None:
print ‘Could not open ‘ + fn
sys.exit(1)
cols = ds.RasterXSize
rows = ds.RasterYSize
bands = ds.RasterCount
print "cols %d " % (cols)
print "rows %d " % (rows)
print "bands %d " % (bands)
[/python]
注意)上記codeで¥¥は全角ですので、半角にするかもしくは、バックスラッシュに置き換えてください。
次に、cmdターミナル上で、ちまちま入力するのはめんどうなので、GUI Shellを使います。
スタート → プログラム → Python 2.7 → IDLE(Python GUI)でShellを動かすと、
な感じでPython SHellが起動します。
メニューから File → Open で先ほどのgdal-example.pyを読み込みこむと、
な感じでfile openされます。
後は、F5 keyを押すか、メニューからRUN → Run Moduleで実行できます。
な感じでRaster dataのcols,rows,bandが表示されます。
とりあえずは問題ないようですね。
後、libraryに何が含まれているかは、helpコマンドでみることができます。
たとえば、osgeo.gdal.Openの引数はどのように呼び出すのか?
な感じで使うことができます。ただし、最初にLibraryをimportしておく必要があります。
Raster data 図化
せっかくですので、raster dataを可視化してみましょう。
作図するには、pylab libraryをinstallします。このpylabはNumpyに含まれているので、Numpy downloadからnumpy-1.7.1-win32-superpack-python2.7.exeをdownloadして、installします。
次に、Ediorで下記のcodeを入力して、保存します。
保存先はf:¥python_example¥chapter_3¥ folderにgdal-example2.pyで保存します。
[python]
from osgeo import gdal, gdalconst
import pylab
# register drivers
gdal.AllRegister()
# create a data set object
dataset = gdal.Open(‘f:¥¥raster_data¥¥example.tif’, gdalconst.GA_ReadOnly)
# create a band object
band = dataset.GetRasterBand(1)
cols = dataset.RasterXSize
rows = dataset.RasterYSize
bands = dataset.RasterCount
print "cols %d " % (cols)
print "rows %d " % (rows)
print "bands %d " % (bands)
# read the band as matrix
dat_matrix = band.ReadAsArray()
# visualize data matrix
pylab.imshow(dat_matrix)
pylab.show()
[/python]
注意)上記codeで¥¥は全角ですので、半角にするかもしくは、バックスラッシュに置き換えてください。
保存後、Python GUI ShellでF5 keyを押すか、メニューからRUN → Run Moduleで実行すると、
な感じでraster dataを表示することができます。
簡単に表示できますね。
その他の使い方については、Pythonでラスターに詳しく書かれています。非常に参考になります。
今回のまとめ
今回は、Python GUI Shellを使って、超~~簡単なGDALのexampleを作成しました。Numpyもinstallして作図も行いました。
次回は、PythonでのProjection(座標系、測地系)の取り扱いについて学んでいきます。
また、本tutorialは、Python Scriptの基本的なことはある程度理解している前提で今後も話を進めていきます。また、誤字、脱字、spell間違いや勘違いも多々出てくると考えられます。
それは違うじゃん!!とかいろんな意見をいただければと思います。
そこんところ ヨロシク~~!!
最近のコメント