PythonとGDAL/OGRによる地理座標の投影座標への変換

/ Python GDAL

GDALとはGeospatial Data Abstraction Libraryの略であり、もともとはC/C++のライブラリであるがPythonやJavaなどで使えるバインディングも用意されている。これを用いて地理座標を投影座標に変換する方法を書いておく。当然逆も可能。

準備

  • Python 3.8.1
  • GDAL 3.0.1

とう状況のもとで行った。Pythonが使えることは前提としてGDALについてはpipかもしくはWindowsであればUnofficial Windows Binaryからダウンロードしてインストールできる。

EPSGコードを使って変換

from osgeo import ogr, osr, gdal

src_srs, dst_srs = osr.SpatialReference(), osr.SpatialReference()
src_srs.ImportFromEPSG(4612) # from EPSG6668
dst_srs.ImportFromEPSG(6674) # to   EPSG6674
trans = osr.CoordinateTransformation(src_srs, dst_srs)

print(trans.TransformPoint(35, 135))
# (-110481.75032585638, -91280.6397650397, 0.0)

EPSG6668はJGD2011の地理座標、EPSG6674はJGD2011の平面直角座標に投影した座標の6系。これで北緯35度、東経135度という地理座標系を投影座標系に変換すると-110481.75032585638, -91280.6397650397という結果になる。

Shapefileから座標系を取得して変換

ファイルの指定する座標系に合わせたいという場合がある。

from osgeo import ogr, osr, gdal

shapefile_name = 'yourshapefile.shp'

shp = ogr.GetDriverByName('ESRI Shapefile').Open(shapefile_name, 0) # '0' means 'read'
src_srs, dst_srs = shp.GetLayer().GetSpatialRef(), osr.SpatialReference()
dst_srs.ImportFromEPSG(6674)
trans = osr.CoordinateTransformation(src_srs, dst_srs)

print(trans.TransformPoint(35, 135))
# (-110481.75032585638, -91280.6397650397, 0.0)

Shapefileといっても実際に座標系の情報が含まれているのは同名ファイルの拡張子.prjであるので、.prjを直接読み込んでImportFromWktで読み込むことも可能だがShapefileありきであればこの方がスマートそうである。少しだけ問題があって、TransformPointの引数の順番が緯度、経度の順番なのかその逆なのかがEPSGコードだけでは定まらない様子である。よって、場合によってはEPSG6668に合わせた順番である緯度、経度の順で引数を指定すると読み込んでいるShapefileによってはinfが返される。その出力側の座標系の範囲外の座標が渡されるとこうなるようだ。

この例では入力側の座標系をShapefileから読んでいるが出力側も同じようにしてShapefileから読み込める。ちなみにz座標はデフォルトで0になっているが指定することもできる。

複数の点を変換

TransformPointsを使うとできる。ただ、こうなってくるとShapefileごと座標系を変換して一時ファイルに書き出すとか、あるいはogr2ogrであらかじめ所望の座標系に変えておくかすればよいと思われる。引数はおそらくタプルの配列だろうけどよくわかっていないので放置しておく。気が向いたら（いずれ使うことがあれば）メモがてら載せたいと思う。

まとめ

ウェブメルカトル座標のバウンディングボックスの4点を変換したり、とか少数の点を変換する際には使えるのではないでしょうか。