PythonでGeoSpatial Dataの読み書きについて_3[Chapter 11] に引き続き、GeoSpatial Dataの測地系と投影法の変更の基本を学んでおきます。

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

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

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

GISのデータは多くの測地系や投影法があり、いつも自分の使いたい座標系ではないことが多いですが、この測地系と投影法の座標変換を学んでおくことは非常に有用です。

ここでは、古い測地系を新しい測地系に変換する例を学びます。
ただし、例としてはUSAの測地系ですので、はっきり言ってよくわかりませんが、とりあえずは書籍の通り実行してみましょう。

測地系NAD27は1927年に設定されたもので、それがNAD83と取り替えられた1980年代までの北アメリカの地理空間分析に一般に使用されました。ウェブサイトに記述されるように、ここでは、NAD27データを使用します。

ESRIのCensus 2000 TIGER/Line Dataからdownloadします。

上記のサイトへ行って、 Preview and Download linkをclickして、次に、ドロップダウン・メニューから Alaskaを選んでください。Line Features – Roads layerを選択して、Submit buttonをclickします。

at_tigeresri484022272.zipデータを最終的にはdownloadします。

そして、at_tigeresri577572118.zipを解凍して、その中のlkA02020.zipを解凍するして、F:￥shapefile_data￥lkA02020 folderに解凍したfileをcopyします。
copyするfileは

1. tgr02020lkA.shp
2. tgr02020lkA.shx
3. tgr02020lkA.dbf

の3つのfileです。

ウェブサイトに記述されるように、NAD27データを使用します。
もしこのShapefileがWSG84データを使用したなら、すべての特徴は間違った場所にあるでしょう:

のように位置のずれが生じています。

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

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

1.library 読み込み

import os, os.path, shutil
from osgeo import ogr
from osgeo import osr

2.入出力Shapefileのprojection定義

# Define the source and destination datums, and a
# transformation object to convert from one to the other.
 
srcDatum = osr.SpatialReference()
srcDatum.SetWellKnownGeogCS('NAD27')
 
dstDatum = osr.SpatialReference()
dstDatum.SetWellKnownGeogCS('WGS84')
 
transform = osr.CoordinateTransformation(srcDatum, dstDatum)

入力ShapefileのprojectionをsrcDatumで定義します。そしてSetWellKnownGeogCS('NAD27')とします。
出力ShapefileのdstDatumをSpatialReferenceで定義します。そしてSetWellKnownGeogCS('WGS84')で世界測地系とします。
transformでCoordinateTransformationを定義します。

3.入力Shapefileのopen

# Open the source shapefile.
 
inPath = "f:￥￥shapefile_data￥￥lkA02020"
srcFile = ogr.Open(inPath + "￥￥" +"tgr02020lkA.shp")
srcLayer = srcFile.GetLayer(0)

4.出力先folderの作成

# Create the dest shapefile, and give it the new projection.
 
outPath = "f:￥￥shapefile_data￥￥lkA-reprojected"
if os.path.exists(outPath):
    shutil.rmtree(outPath)
os.mkdir(outPath)

5.出力Shapefileの作成

driver = ogr.GetDriverByName("ESRI Shapefile")
dstPath = os.path.join(outPath, "lkA02020.shp")
dstFile = driver.CreateDataSource(dstPath)
dstLayer = dstFile.CreateLayer("layer", dstDatum)

PythonでGeoSpatial Dataの読み書き_2[Chapter 9]を参考にして下さい。

6.Projectionの変換

for i in range(srcLayer.GetFeatureCount()):
    feature = srcLayer.GetFeature(i)
    geometry = feature.GetGeometryRef()
 
    newGeometry = geometry.Clone()
    newGeometry.Transform(transform)
 
    feature = ogr.Feature(dstLayer.GetLayerDefn())
    feature.SetGeometry(newGeometry)
    dstLayer.CreateFeature(feature)
    feature.Destroy()

srcLayerのfeature数でloopします。
featureとgeometryに入力Shapefileの情報をcopyして、newGeometryのcloneを作成して、Transform(transform)で投影法の変換を行います。
属性値については変更はありませんので、そのまま出力Shapefileに追加していきます。
feature.SetGeometryで変換後のgeomertyを追加して、dstLayerに新しくfeatureを作成して追加していきます。

7.Objectの破棄

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

実行後、QGISで確認してみてください。

でWGS84に変更されています。

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

objectの引数が不明な場合は、PyScriptのインタプリタでhelp(osr)などを入力すると、使用方法がインタプリタ画面に出力されます。

今回は、PythonでのGeoSpatial dataのShapefileの測地系の変換を行いました。
ExercisesでJGD2000の測地系を使って、Shapefileの変換をtyrしてください。
次回は、PythonでGeoSpatial Dataの表示と保管をやってみます。

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

Python GeoSpatial Tutorialの目次に戻る。