PythonでGeoSpatial Dataの読み書きについて_2[Chapter 9]
PythonでGeoSpatial Dataの読み書きについて_1[Chapter 8] に引き続き、geospatial dataの書き出しの基本を学んでおきます。
本tutorialは、Python Geospatial Developmentを参考にして、実践的に使えるPythonでGeoSpatialの使い方を学んでいくものです。
また、本Blog中のsource code等に関しては、あくまでも参考としてください。なにがあっても責任とれませんので。
そこんところ、ヨロシク~~!!
準備
PythonでGeoSpatial Dataの読み書きについて_1[Chapter 8] でPythonとGeoSpatialに関するlibraryが設定できているもとして進めていきます。
Writing GeoSpatial data
PythonでGeoSpatial Dataの読み書きについて_1[Chapter 8] のcodeに、計算した「bounding box」をShapefileを書き出すことをやってみます。
世界の国のbounding boxをShapefileに書き出す
PythonでGeoSpatial Dataの読み書きについて_1[Chapter 8] では、単純にインタプリタ出力でしたが、ここでは、その「bounding box」の値を使って、polygonを作成して、Shapefileを作成してみます。
このexample codeはShapefile作成時の基本となります。
Shapefile writing code
とりあえず、Ediorで下記のcodeを入力して、保存します。
保存先はf:¥python_example¥chapter_5¥ folderにboundingBoxesToShapefile.pyで保存します。
[python]
import os, os.path, shutil
from osgeo import ogr
from osgeo import osr
# Open the source shapefile.
srcFile = osgeo.ogr.Open("f:¥¥shapefile_data¥¥TM_WORLD_BORDERS-0.3.shp")
srcLayer = srcFile.GetLayer(0)
# Open the output shapefile.
if os.path.exists("bounding-boxes"):
shutil.rmtree("bounding-boxes")
os.mkdir("bounding-boxes")
spatialReference = osgeo.osr.SpatialReference()
spatialReference.SetWellKnownGeogCS(‘WGS84′)
driver = osgeo.ogr.GetDriverByName("ESRI Shapefile")
dstPath = os.path.join("bounding-boxes", "boundingBoxes.shp")
dstFile = driver.CreateDataSource(dstPath)
dstLayer = dstFile.CreateLayer("layer", spatialReference)
fieldDef = osgeo.ogr.FieldDefn("COUNTRY", osgeo.ogr.OFTString)
fieldDef.SetWidth(50)
dstLayer.CreateField(fieldDef)
fieldDef = osgeo.ogr.FieldDefn("CODE", osgeo.ogr.OFTString)
fieldDef.SetWidth(3)
dstLayer.CreateField(fieldDef)
# Read the country features from the source shapefile.
for i in range(srcLayer.GetFeatureCount()):
feature = srcLayer.GetFeature(i)
countryCode = feature.GetField("ISO3")
countryName = feature.GetField("NAME")
geometry = feature.GetGeometryRef()
minLong,maxLong,minLat,maxLat = geometry.GetEnvelope()
# Save the bounding box as a feature in the output
# shapefile.
linearRing = osgeo.ogr.Geometry(osgeo.ogr.wkbLinearRing)
linearRing.AddPoint(minLong, minLat)
linearRing.AddPoint(maxLong, minLat)
linearRing.AddPoint(maxLong, maxLat)
linearRing.AddPoint(minLong, maxLat)
linearRing.AddPoint(minLong, minLat)
polygon = osgeo.ogr.Geometry(osgeo.ogr.wkbPolygon)
polygon.AddGeometry(linearRing)
feature = osgeo.ogr.Feature(dstLayer.GetLayerDefn())
feature.SetGeometry(polygon)
feature.SetField("COUNTRY", countryName)
feature.SetField("CODE", countryCode)
dstLayer.CreateFeature(feature)
feature.Destroy()
# All done.
srcFile.Destroy()
dstFile.Destroy()
[/python]
注意)上記codeで¥は全角ですので、半角にするかもしくは、バックスラッシュに置き換えてください。
PyScriptでrunすると、何も起きませんが、f:¥python_example¥chapter_8¥のfolderにbounding-boxes folderが作成されて、Shapefileが作成されています。
たとえば、日本だけを抜き出すと、
な感じで、「bounding box」がpolygonとして作成されてることがわかります。
Shapefileのviewとしては、QGISを使用しています。
codeの説明
boundingBoxesToShapefile.pyが何を行っているのかを順にみていきましょう。
1.library 読み込み
import os, os.path, shutil from osgeo import ogr from osgeo import osr
ここで、osrは空間参照を行うために必要なlibraryです。
2.入力ShapefileのOpenとLayerの定義
# Open the source shapefile. srcFile = osgeo.ogr.Open("f:¥¥shapefile_data¥¥TM_WORLD_BORDERS-0.3.shp") srcLayer = srcFile.GetLayer(0)
Shapefileのlayerはただ一つのlayerですから、GetLayer(0)となります。
3.出力先folderの作成
# Open the output shapefile. if os.path.exists("bounding-boxes"): shutil.rmtree("bounding-boxes") os.mkdir("bounding-boxes")
4.出力Shapefileの作成
spatialReference = osgeo.osr.SpatialReference() spatialReference.SetWellKnownGeogCS('WGS84') driver = osgeo.ogr.GetDriverByName("ESRI Shapefile") dstPath = os.path.join("bounding-boxes", "boundingBoxes.shp") dstFile = driver.CreateDataSource(dstPath) dstLayer = dstFile.CreateLayer("layer", spatialReference)
spatialReferenceで空間参照を作成します。ここではWGS84として定義します。
driverで出力形式をShapfileとします。
dstFileでShapefileを作成して、dstLayerでlayerを定義します。
5.属性dbfのfield作成
fieldDef = osgeo.ogr.FieldDefn("COUNTRY", osgeo.ogr.OFTString) fieldDef.SetWidth(50) dstLayer.CreateField(fieldDef) fieldDef = osgeo.ogr.FieldDefn("CODE", osgeo.ogr.OFTString) fieldDef.SetWidth(3) dstLayer.CreateField(fieldDef)
属性データの作成をfieldDefで行います。ここでは、countryとcodeのfieldを作成します。
6.countriesのbounding boxの計算
# Read the country features from the source shapefile. for i in range(srcLayer.GetFeatureCount()): feature = srcLayer.GetFeature(i) countryCode = feature.GetField("ISO3") countryName = feature.GetField("NAME") geometry = feature.GetGeometryRef() minLong,maxLong,minLat,maxLat = geometry.GetEnvelope()
srcLayer.GetFeatureCount()でデータ数分for文でloopします。ここでは国の数になります。
featureはここでは各国となります。Shapefileはfeatureの集合と考えると良いでしょう。
GetFieldで属性値を取り出します。ここでは、countryCodeとcountryNameを取り出しています。
feature.GetGeometryRef()でgeometryを取得して、 geometry.GetEnvelope()で最大最小の経度緯度を取り出します。
7.出力ShaefileのFeature作成
# Save the bounding box as a feature in the output # shapefile. linearRing = osgeo.ogr.Geometry(osgeo.ogr.wkbLinearRing) linearRing.AddPoint(minLong, minLat) linearRing.AddPoint(maxLong, minLat) linearRing.AddPoint(maxLong, maxLat) linearRing.AddPoint(minLong, maxLat) linearRing.AddPoint(minLong, minLat) polygon = osgeo.ogr.Geometry(osgeo.ogr.wkbPolygon) polygon.AddGeometry(linearRing)
linearRingでpointを追加します。追加する順番は、
になります。
point追加後に、各pointを結ぶpolygonを作成して、AddGeometryで追加していきます。
pointは4点ですが、polygonは閉じてなければならないので、5点目に最初のpointを追加します。
8.属性dbfの追加
feature = osgeo.ogr.Feature(dstLayer.GetLayerDefn()) feature.SetGeometry(polygon) feature.SetField("COUNTRY", countryName) feature.SetField("CODE", countryCode) dstLayer.CreateFeature(feature) feature.Destroy()
出力Shapefileに属性dbfを書き込んでいきます。
9.Objectの削除
# All done. srcFile.Destroy() dstFile.Destroy()
となります。順番にみていくと理解できると思います。
今回のまとめ
今回は、PythonでのGeoSpatial dataを読み込んで「bounding box」を取得してShapefileを作成しました。
Exercisesで世界各国のbounding boxを出力しましたが、日本だけとか特定の国を出力できるように変更してください。
次回は、raster dataを扱ってみます。
また、本tutorialは、Python Scriptの基本的なことはある程度理解している前提で今後も話を進めていきます。また、誤字、脱字、spell間違いや勘違いや翻訳間違いなど多々出てくると考えられます。
それは違うじゃん!!とかいろんな意見をいただければと思います。
そこんところ ヨロシク~~!!
最近のコメント