Daily Archives: 07/13/2013

PythonでGeoSpatial Dataの読み書きについて_3[Chapter 10]

PythonでGeoSpatial Dataの読み書きについて_3[Chapter 10]

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が設定できているもとして進めていきます。

DEMを使った標高値のヒストグラム

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

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

DEM data download

DEMデータをNOAA DEM からdownloadします。
blog.godo-tys.jp_wp-content_gallery_python_10_image01.jpg

上記のNOAAのHPのGet Data Onlineをclickすると、
blog.godo-tys.jp_wp-content_gallery_python_10_image02.jpg
のタイルが表示されます。ニュージーランドを含む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データの作成が終わりました。

標高値のヒストグラフ作成 code

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

[python]
import sys, struct
from osgeo import gdal
from osgeo import gdalconst

import matplotlib.pyplot as plt # import the Matplotlib module

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

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

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

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)

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

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 height in sorted(histogram.keys()):
print height,histogram[height]
yy.append(height)
xx.append(histogram[height])

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

[/python]

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

PyScriptでrunすると、標高値とその標高値のカウント数がグラフで、
blog.godo-tys.jp_wp-content_gallery_python_10_image03.jpg

な感じで、表示されます。
このmatplotlibは非常に強力なlibraryです。詳しくは、matplotlibを参照してください。これ使いこなせれば、本当にすごいです。

codeの説明

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

ニュージーランドを含む最大最小の経度緯度で、
blog.godo-tys.jp_wp-content_gallery_python_10_image04.jpg

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

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

インタプリタに標高値とカウント数が
blog.godo-tys.jp_wp-content_gallery_python_10_image05.jpg

の感じに表示されます。

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)などを入力すると、
blog.godo-tys.jp_wp-content_gallery_python_10_image06.jpg

な感じで表示されます。

今回のまとめ

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

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

Python GeoSpatial Tutorialの目次に戻る。

1 / 11

Social Widgets powered by AB-WebLog.com.