PythonでGeoSpatial Dataの読み書きについて_2[Chapter 9] に引き続き、geospatial raster dataの基本を学んでおきます。

本tutorialは、Python Geospatial Developmentを参考にして、実践的に使えるPythonでGeoSpatialの使い方を学んでいくものです。

また、本Blog中のsource code等に関しては、あくまでも参考としてください。なにがあっても責任とれませんので。
そこんところ、ヨロシク～～！！

PythonでGeoSpatial Dataの読み書きについて_1[Chapter 8] でPythonとGeoSpatialに関するlibraryが設定できているもとして進めていきます。

Rasterデータとして代表的なDEM(Digital Elevation Map)を使ってみましょう。
DEMファイルがは標高データを含んでいるので、与えられたエリアに対する高さを分析することは面白いかもしれません。
例えば、私たちは、ある国の領域にどれくらいの標高があるか示すヒストグラムを引くことができまます。
GLOBEデータセットからあるDEMデータをとり、そのデータを使用して、標高ヒストグラムを計算してみましょう。

今回対象とする国は、ニュージーランドです。

DEMデータをNOAA DEM からdownloadします。

上記のNOAAのHPのGet Data Onlineをclickすると、

のタイルが表示されます。ニュージーランドを含むLをdownloadします。
Downloadするファイルは、l10g.zipです。

DEMファイルを扱うには、georeferenceが必要になります。これは、http://www.ngdc.noaa.gov/mgg/topo/elev/esri/hdrから、先ほどのl10gに対応するl10g.hdrをdownloadします。

そして、l10g.zipを解凍してF:￥raster_data￥l10gand folderにl10g.hdrとともにcopyします。
とりあえずは、DEMデータの作成が終わりました。

とりあえず、Ediorで下記のcodeを入力して、保存します。
保存先はf:￥python_example￥chapter_5￥ folderにhistogram.pyで保存します。
今回は、matplotlib libraryを使うので、matplotlib-1.2.1.win32-py2.7をdownloadして、installしておいてください。

histogram.pyが何を行っているのかを順にみていきましょう。

1.library 読み込み

import sys, struct
from osgeo import gdal
from osgeo import gdalconst
 
import matplotlib.pyplot as plt  # import the Matplotlib module

ここで、matplotlibは作図を行うために必要なlibraryです。

2.ニュージーランドのboubding boxの値

minLat  = -48
maxLat  = -33
minLong = 165
maxLong = 179

ニュージーランドを含む最大最小の経度緯度で、

の図のような感じになります。

3.DEMデータファイルのopen

dataset = gdal.Open("f:￥￥raster_data￥￥l10g￥￥l10g")
band = dataset.GetRasterBand(1)

GetRasterBand(1)でband 1を使います。

4.アフィン逆変換

t = dataset.GetGeoTransform()
success,tInverse = gdal.InvGeoTransform(t)
if not success:
    print "Failed!"
    sys.exit(1)

datasetをアフィン逆変換します。　失敗したらporgramは終了します。

5.最大値と最小値の計算

x1,y1 = gdal.ApplyGeoTransform(tInverse, minLong, minLat)
x2,y2 = gdal.ApplyGeoTransform(tInverse, maxLong, maxLat)
 
minX = int(min(x1, x2))
maxX = int(max(x1, x2))
minY = int(min(y1, y2))
maxY = int(max(y1, y2))
 
width = (maxX - minX) + 1
fmt = "<" + ("h" * width)

ニュージーランドのbounding boxの経度緯度はアフィン逆変換して、DEMデータメッシュ番号の最大最小値を計算します。

6.listの定義

histogram = {} # Maps height to # pixels with that height.
xx = []
yy = []

7.標高値のヒストグラム作成

for y in range(minY, maxY+1):
    scanline = band.ReadRaster( minX, y, width, 1, width, 1, band.DataType)
 
    values = struct.unpack(fmt, scanline)
 
    for value in values:
        if value != band.GetNoDataValue():
            try:
                histogram[value] += 1
            except KeyError:
                histogram[value] = 1

for y in range(minY, maxY+1)で移動方向のDEM呼び出しをloopします。
band.ReadRasterでDEMファイルを必要な部分を呼び出します。
struct.unpack(fmt, scanline)でscanlineから必要な標高値を取り出します。
取り出したvaluesの個数分だけloopして、標高値とその標高値のカウント数を計算して、histogram[value]に代入します。

8.結果のインタプリタ出力

for height in sorted(histogram.keys()):
    print height,histogram[height]
    yy.append(height)
    xx.append(histogram[height])

インタプリタに標高値とカウント数が

の感じに表示されます。

9.作図の出力

# plot
plt.plot(yy,xx)
plt.ylabel('Height')
plt.xlabel('Count of height')
plt.show()

xy plot図を表示します。

この作図画面を閉じるとprogramが終了します。

順番にみていくと理解できると思います。

objectの引数が不明な場合は、PyScriptのインタプリタでhelp(gdal.Band.ReadRaster)などを入力すると、

な感じで表示されます。

今回は、PythonでのGeoSpatial dataのraster dataを使ってヒストグラムを作成してみました。
Exercisesで世界各国のbounding boxを使って、日本の標高値のヒストグラムを出力できるように変更してください。
次回は、測地系と投影法の変換をやってみます。

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

Python GeoSpatial Tutorialの目次に戻る。