Pythonでgeometryと距離の変換と標準化_1[Chapter 16]

Pythonでgeometryと距離の変換と標準化_1[Chapter 16]

PythonでShapely geometriesを使ってみる[Chapter 15] に引き続き、pyproj libraryを使用した距離計算の基本を学んでおきます。

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

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

準備

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

geometryと距離の変換

今、下図のような2点間の距離を求めることを考えます。
blog.godo-tys.jp_wp-content_gallery_python_16_image01.jpg

緯度経度で2点間の距離が表されている場合は、直線距離は角距離を求めることになります。平面直角座標系であれば、2点間の距離は単純に求めることができます。
ここでは、緯度経度(WGS84)の距離を求めてみます。

まずは、PythonでGeoSpatial Dataの表示と保存_1[Chapter 13] でタイとミャンマーの国境を抽出したshapefileのlineデータから国境の総延長を求めるexample codeを試してみましょう。

Sample data

PythonでGeoSpatial Dataの表示と保存_1[Chapter 13] で作成したborder.shpを使います。

国境の総延長計算 example code

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

import os.path
import osgeo.ogr
import pyproj

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

inPath = "f:¥¥shapefile_data¥¥common-border"
filename = os.path.join(inPath, "border.shp")
shapefile = osgeo.ogr.Open(filename)
layer = shapefile.GetLayer(0)
feature = layer.GetFeature(0)
geometry = feature.GetGeometryRef()
segments = getLineSegmentsFromGeometry(geometry)

geod = pyproj.Geod(ellps='WGS84')

totLength = 0.0
for segment in segments:
    for i in range(len(segment)-1):
        pt1 = segment[i]
        pt2 = segment[i+1]

        long1,lat1 = pt1
        long2,lat2 = pt2

        angle1,angle2,distance = geod.inv(long1, lat1,
                                          long2, lat2)
        totLength += distance

print "Total border length = %0.2f km" % (totLength/1000)

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

PyScriptでrunすると、python インタプリタに、
blog.godo-tys.jp_wp-content_gallery_python_16_image02.jpg

でタイとミャンマーの国境の総延長の計算結果が表示されます。

codeの説明

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

1.library 読み込み

import os.path
import osgeo.ogr
import pyproj

pyprojlibraryを使います。

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.入力Shapefileのopenとsegment作成

inPath = "f:¥¥shapefile_data¥¥common-border"
filename = os.path.join(inPath, "border.shp")
shapefile = osgeo.ogr.Open(filename)
layer = shapefile.GetLayer(0)
feature = layer.GetFeature(0)
geometry = feature.GetGeometryRef()
segments = getLineSegmentsFromGeometry(geometry)

Shapefileのlayerはただ一つのlayerですから、GetLayer(0)となります。
featureからgeometryを作って、getLineSegmentsFromGeometry関数でsegmentを作成します。

4.geodの定義

geod = pyproj.Geod(ellps='WGS84')

読み込んだshapefileのprojectionはWGS84です。

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

totLength = 0.0
for segment in segments:
    for i in range(len(segment)-1):
        pt1 = segment[i]
        pt2 = segment[i+1]
 
        long1,lat1 = pt1
        long2,lat2 = pt2
 
        angle1,angle2,distance = geod.inv(long1, lat1,
                                          long2, lat2)
        totLength += distance

segmentの数でloopします。
そして、segmentごとの緯度、経度から角度と距離をgeod.inv(…)関数で求めます。

geod.inv関数とは、PyScripterのインタプリタ上でhelp(geod.inv)
blog.godo-tys.jp_wp-content_gallery_python_16_image03.jpg

のような引数と関数であることがわかります。

7.結果の表示

print "Total border length = %0.2f km" % (totLength/1000)

タイとミャンマーの国境の総延長距離の出力です。

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

ここで初めてdef関数objectが出てきましたが、fortranではsubroutineのようなものですね。

今回のまとめ

今回は、Pythonでpyprojを使って距離の計算について学びました。pyproj libraryは機能が非常に豊富ですので、ぜひ、pyproj documentationで自習することをお勧めします。
次回は、もう一つpyproj libraryを使ったexampleで距離から緯度経度を計算する方法について学びます。

また、本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.