Pythonでgeometryと距離の変換と標準化_2[Chapter 17]

Pythonでgeometryと距離の変換と標準化_2[Chapter 17]

Pythonでgeometryと距離の変換と標準化_2[Chapter 16] に引き続き、pyproj libraryを使用した距離から緯度経度を求める計算の基本を学んでおきます。

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

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

準備

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

距離から緯度経度計算

前回http://blog.godo-tys.jp/2013/08/02/3525/|PythonでShapely geometriesを使ってみる[Chapter 15]で使ったCA_Features_XXX.txtのdataからCaliforniaのLas Vegasの小さな町のShoshoneから、西へ132.7kmの場所の緯度経度を計算してみます。

Sample data

PythonでShapely geometriesを使ってみる[Chapter 15] でdownloadしたhttp://www2.census.gov/geo/tiger/2013CBSA/CBSAs_Feb2013_delineation.zipを使います。

私の開発環境では、解凍後、f:¥CA_Features_20130602¥のfolderを作成して、解凍したfileをcopyしておきます。

距離から緯度経度計算example code

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

# findShoshone.py

import pyproj

distance = 132.7 * 1000
angle    = 270.0

f = file("f:¥¥CA_Features_20130602¥¥CA_Features_20130602.txt", "r")
for line in f.readlines():
    chunks = line.rstrip().split("|")
    if chunks[1] == "Shoshone" and ¥
       chunks[2] == "Populated Place":
        latitude = float(chunks[9])
        longitude = float(chunks[10])

        geod = pyproj.Geod(ellps='WGS84')
        newLong,newLat,invAngle = geod.fwd(longitude, latitude, angle, distance)

        print "Shoshone is at %0.4f,%0.4f" % (latitude, longitude)
        print "The point %0.2f km west of Shoshone " % (distance/1000.0) ¥
            + "is at %0.4f, %0.4f" % (newLat, newLong)

f.close()

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

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

で、Shoshoneから西へ132.7kmの場所の緯度経度が表示されます。

codeの説明

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

1.library 読み込み

import pyproj

pyproj libraryを使います。

2.変数の定義

distance = 132.7 * 1000
angle    = 270.0

距離は132.7km、西は北を基準に270度を設定する。

3.緯度経度計算

f = file("f:¥¥CA_Features_20130602¥¥CA_Features_20130602.txt", "r")
for line in f.readlines():
    chunks = line.rstrip().split("|")
    if chunks[1] == "Shoshone" and ¥
       chunks[2] == "Populated Place":
        latitude = float(chunks[9])
        longitude = float(chunks[10])
 
        geod = pyproj.Geod(ellps='WGS84')
        newLong,newLat,invAngle = geod.fwd(longitude, latitude, angle, distance)
 
        print "Shoshone is at %0.4f,%0.4f" % (latitude, longitude)
        print "The point %0.2f km west of Shoshone " % (distance/1000.0) ¥
            + "is at %0.4f, %0.4f" % (newLat, newLong)

Test fileを開いて、ShoshoneとPopulated Placeの行を取得して、緯度経度を取得します。
pyproj.Geod(ellps='WGS84')でgeodを定義して、geod.fwd(…)で緯度経度を計算します。
そして、画面上に結果を出力します。
距離と角度から緯度、経度をgeod.fwd(…)関数で求めます。

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

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

4.fileを閉じる

f.close()

読み込んだtext fileを閉じる。

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

今回のまとめ

今回は、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.