PythonでGeoSpatial Dataの読み書きについて_3[Chapter 12] に引き続き、GeoSpatial DataのVector dataの表示(抽出)の基本を学んでおきます。

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

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

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

地理空間情報を扱う場合、多くはShapefileで提供されることが多いですが、Shapefilesは不適当であるか、能率が悪いことがあります。
これらの状況で、効率よくPython libraryで地理空間情報を扱ったり、保存したりできる独立したFormatを持つ必要があります。このformat(地理空間情報データ用の共通語のためのvector format)は、 Well-Known TextあるいはWKTと呼ばれます。

WKTは、point、lineあるいはpolygonのような地理空間情報のobjectのコンパクトなテキストベースの記述です。
例えば、バチカン市国のWKTは

で表されます。
また、WKTは、geometryがどんなに複雑になってもテキストベースです。

WKTは、投影法、データおよび(または)座標系を包含する空間的基準を表わすためにも使用することができます。
例えば、ここに、osgeoがあります。
WKTに変換されたWGS84データを使用して、測地系を表わすosr.SpatialReferenceオブジェクトは、

として表現されます。
WKT formatのspatial referenceはpython libraryを使って、他のspatial referenceに変換するとができます。

ここでは、shapely libraryをつかってタイとミャンマーの国境のintersectionをやってみます。

sample dataのTM_WORLD_BORDERS-0.3からView raw fileをclickしてdownload、その後解凍しておきます。
私の開発環境では、f:￥shapefile_dataのfolderを作成して、解凍したfileをcopyしておきます。

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

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

1.library 読み込み

import os,os.path,shutil
 
import osgeo.ogr
import shapely.wkt

今回は、shapely libraryを使います。
shapelyについては、shapely manualを参考にしてください。

2.入力Shapefileのopen

# Load the thai and myanmar polygons from the world borders
# dataset.
 
shapefile = osgeo.ogr.Open("f:￥￥shapefile_data￥￥TM_WORLD_BORDERS-0.3.shp")
layer = shapefile.GetLayer(0)

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

3.タイとミャンマーのObjectの定義とgeometry作成

thailand = None
myanmar = None
 
for i in range(layer.GetFeatureCount()):
    feature = layer.GetFeature(i)
    if feature.GetField("ISO2") == "TH":
        geometry = feature.GetGeometryRef()
        thailand = shapely.wkt.loads(geometry.ExportToWkt())
    elif feature.GetField("ISO2") == "MM":
        geometry = feature.GetGeometryRef()
        myanmar = shapely.wkt.loads(geometry.ExportToWkt())

タイとミャンマーのfeatureをそれぞれthailandとmyanmarにgeometry.ExportToWkt()で追加します。

4.国境の抽出

# Calculate the common border.
 
commonBorder = thailand.intersection(myanmar)

intersectionでcommonBorderを抽出します。

5.出力先folderの作成

# Save the common border into a new shapefile.
 
outPath = "f:￥￥shapefile_data￥￥common-border"
if os.path.exists(outPath):
    shutil.rmtree(outPath)
os.mkdir(outPath)

6.出力Shapefileの作成

spatialReference = osgeo.osr.SpatialReference()
spatialReference.SetWellKnownGeogCS('WGS84')
 
driver = osgeo.ogr.GetDriverByName("ESRI Shapefile")
dstPath = os.path.join(outPath, "border.shp")
dstFile = driver.CreateDataSource(dstPath)
dstLayer = dstFile.CreateLayer("layer", spatialReference)

spatialReferenceを”WGS84”として出力Shapefileを定義します。
PythonでGeoSpatial Dataの読み書き_2[Chapter 9] を参考にして下さい。

7.出力Shapefileにfeatureの作成

wkt = shapely.wkt.dumps(commonBorder)
 
feature = osgeo.ogr.Feature(dstLayer.GetLayerDefn())
feature.SetGeometry(osgeo.ogr.CreateGeometryFromWkt(wkt))
dstLayer.CreateFeature(feature)

commonBorderをshapely.wkt.dumpsでwtkにdumpします。
featureを定義して、wktからGeometryを作成してfeatureにsetします。
dstLayerの出力Shqpefileにfeatureを作成します。

8.Objectの破棄

# All done
 
feature.Destroy()
dstFile.Destroy()

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

今回は、PythonでのGeoSpatial dataの表示について、WKTの基本を学びました。
Exercisesでshapely manualのshapelyのexampleを実行して、shapelyでできることを確認すると良いと思います。
次回は、PythonでGeoSpatial Dataの保存をやってみます。

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

Python GeoSpatial Tutorialの目次に戻る。