PythonでGDALをちょろっと使ってみる。[Chapter 3]

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によれば、
blog.godo-tys.jp_wp-content_gallery_python_03_gdal_image.jpg

の図のような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で保存します。

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)

注意)上記codeで¥¥は全角ですので、半角にするかもしくは、バックスラッシュに置き換えてください。

次に、cmdターミナル上で、ちまちま入力するのはめんどうなので、GUI Shellを使います。
スタート → プログラム → Python 2.7 → IDLE(Python GUI)でShellを動かすと、
blog.godo-tys.jp_wp-content_gallery_python_03_image01.jpg

な感じでPython SHellが起動します。

メニューから File → Open で先ほどのgdal-example.pyを読み込みこむと、
blog.godo-tys.jp_wp-content_gallery_python_03_image02.jpg

な感じでfile openされます。

後は、F5 keyを押すか、メニューからRUN → Run Moduleで実行できます。

では、実行してみると、
blog.godo-tys.jp_wp-content_gallery_python_03_image03.jpg

な感じでRaster dataのcols,rows,bandが表示されます。

とりあえずは問題ないようですね。

後、libraryに何が含まれているかは、helpコマンドでみることができます。
たとえば、osgeo.gdal.Openの引数はどのように呼び出すのか?
blog.godo-tys.jp_wp-content_gallery_python_03_image04.jpg

な感じで使うことができます。ただし、最初にLibraryをimportしておく必要があります。

Raster data 図化

せっかくですので、raster dataを可視化してみましょう。
作図するには、pylab libraryをinstallします。このpylabNumpyに含まれているので、Numpy downloadからnumpy-1.7.1-win32-superpack-python2.7.exeをdownloadして、installします。

次に、Ediorで下記のcodeを入力して、保存します。
保存先はf:¥python_example¥chapter_3¥ folderにgdal-example2.pyで保存します。

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()

注意)上記codeで¥¥は全角ですので、半角にするかもしくは、バックスラッシュに置き換えてください。

保存後、Python GUI ShellでF5 keyを押すか、メニューからRUN → Run Moduleで実行すると、
blog.godo-tys.jp_wp-content_gallery_python_03_image06.jpg

な感じでraster dataを表示することができます。
簡単に表示できますね。
その他の使い方については、Pythonでラスターに詳しく書かれています。非常に参考になります。

今回のまとめ

今回は、Python GUI Shellを使って、超~~簡単なGDALのexampleを作成しました。Numpyもinstallして作図も行いました。
次回は、PythonでのProjection(座標系、測地系)の取り扱いについて学んでいきます。

また、本tutorialは、Python Scriptの基本的なことはある程度理解している前提で今後も話を進めていきます。また、誤字、脱字、spell間違いや勘違いも多々出てくると考えられます。
それは違うじゃん!!とかいろんな意見をいただければと思います。
そこんところ ヨロシク~~!!

Python GeoSpatial Tutorialの目次に戻る。

Leave a Comment


NOTE - You can use these HTML tags and attributes:
<a href="" title=""> <abbr title=""> <acronym title=""> <b> <blockquote cite=""> <cite> <code> <del datetime=""> <em> <i> <q cite=""> <strike> <strong>

 

WP-SpamFree by Pole Position Marketing

Social Widgets powered by AB-WebLog.com.