Pythonでgeometryと距離の変換と標準化_2[Chapter 17] に引き続き、pyproj libraryを使用した汎用性を考えたcodeを学んでおきます。

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

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

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

今までは、Python program fileにfileの存在場所をハードコーディングして、汎用性がないcodeでした。また、PyScripterのrunボタンを使ってインタープリタ画面に結果を出力していました。

ここで、汎用性を持たせるために、cmdターミナルからPython programを実行し、任意の場所のfileが読み込めるように修正してみましょう。

まずは、sys.argvを使って、

if len(sys.argv) != 2:
    print "Usage: calcFeatureLengths.py <shapefile>"
    sys.exit(1)

のようにします。
これは、cmdターミナルで

python calcFeatureLengths.py test_file.shp

のように実行時にコマンド　パラメータを渡すことを意味します。
この場合、

• sys.argv[0]には、calcFeatureLengths.py
• sys.argv[1]には、test_file.shp

がsys.argvにimportされます。

このsys.argvを使ってexample codeを作成してみましょう。

とりあえず、Ediorで下記のcodeを入力して、保存します。
保存先はf:￥python_example￥chapter_5￥ folderにcalcFeatureLengths.pyで保存します。

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

1.library 読み込み

import sys
from osgeo import ogr, osr
import pyproj

osgeoのogrとosrを使います。
pyproj libraryを使います。

2.def関数objectの定義

def getLineSegmentsFromGeometry(geometry):
    segments = []
    if geometry.GetPointCount() > 0:
        segment = []
        for i in range(geometry.GetPointCount()):
            segment.append(geometry.GetPoint_2D(i))
        segments.append(segment)
    for i in range(geometry.GetGeometryCount()):
        subGeometry = geometry.GetGeometryRef(i)
        segments.extend(
            getLineSegmentsFromGeometry(subGeometry))
    return segments

これは、line geometryからsegmentに分割する関数です。戻り値はsegmentsとなります。
segmentとは、「始点と終点をもった途中に座標点を持たない直線」と考えることができます。
したがって、lineは複数のsegmentで構成されます。

3.coammnd parameterのcheck

if len(sys.argv) != 2:
    print "Usage: calcFeatureLengths.py <shapefile>"
    sys.exit(1)

4.入力Shapefileのopen

filename = sys.argv[1]
shapefile = ogr.Open(filename)
layer = shapefile.GetLayer(0)

Shapefileのlayerはただ一つのlayerですから、GetLayer(0)となります。

5.spatialRefのcheck

spatialRef = layer.GetSpatialRef()
if spatialRef == None:
    print "Shapefile lacks a spatial reference, using WGS84."
    spatialRef = osr.SpatialReference()
    spatialRef.SetWellKnownGeogCS('WGS84')
 
if spatialRef.IsProjected():
    srcProj = pyproj.Proj(spatialRef.ExportToProj4())
    dstProj = pyproj.Proj(proj='longlat', ellps='WGS84',
                          datum='WGS84')

読み込んだshapefileのspatialRefをcheckちます。defaultはWGS84です。
spatialRef.IsProjected()でprojectionの定義をします。

6.国境の総延長距離計算

for i in range(layer.GetFeatureCount()):
    feature = layer.GetFeature(i)
    geometry = feature.GetGeometryRef()
    segments = getLineSegmentsFromGeometry(geometry)
 
    geod = pyproj.Geod(ellps='WGS84')
 
    totLength = 0.0
    for segment in segments:
        for j in range(len(segment)-1):
            pt1 = segment[j]
            pt2 = segment[j+1]
 
            long1,lat1 = pt1
            long2,lat2 = pt2
 
            if spatialRef.IsProjected():
                long1,lat1 = pyproj.transform(srcProj,
                                              dstProj,
                                              long1, lat1)
                long2,lat2 = pyproj.transform(srcProj,
                                              dstProj,
                                              long2, lat2)
 
            try:
                angle1,angle2,distance = geod.inv(long1, lat1,
                                                  long2, lat2)
            except ValueError:
                print "Unable to calculate distance from " ￥
                    + "%0.4f,%0.4f to %0.4f,%0.4f" ￥
                    % (long1, lat1, long2, lat2)
                distance = 0.0
            totLength += distance
 
    print "Total length of feature %d is %0.2f km" ￥
        % (i, totLength/1000)

layer.GetFeatureCount()数でloopします。
そして、getLineSegmentsFromGeometryでsegmentに分割します。
spatialRef.IsProjected()がtrueの場合、pyproj.transformでWGS84の緯度経度に変換します。
分割したsegmentの数だけloopしてgeod.inv(…)で距離を求めます。
そして、国土の総延長距離の計算結果を表示します。

geod.inv関数とは、PyScripterのインタプリタ上でhelp(geod.inv)で

のような引数と関数であることがわかります。
うまくhelpを使うことで、それぞれのobjectのmethodなどを確認することができます。　

7.try…except : statementについて

try:
    angle1,angle2,distance = geod.inv(long1, lat1, long2, lat2)
except ValueError:
    print "Unable to calculate distance from " ￥
          + "%0.4f,%0.4f to %0.4f,%0.4f" ￥
          % (long1, lat1, long2, lat2)
    distance = 0.0

これは、エラー回避のためのstatementです。
もし、「緯度経度変換にエラーがあった場合は、programが異常終了しないように、segmentの距離を0とします。」
と言う意味になります。
汎用性を持たす場合は、I/O checkやerror checkを行い、programの異常終了等がないようにしなければなりません。

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

今回は、pyproj libraryを使った汎用性のあるexample codeについて学びました。これで、5章が終了しました。
次回はいよいよGeoSpatialでDatabaseを使う際の準備を行います。

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

Python GeoSpatial Tutorialの目次に戻る。