【PythonとGIS】GDALを使ってShapeファイルを作成してみる

皆さん、お久りぶりです。ずっとブログ更新をさぼってしまっていましたが、またちょくちょく書いていこうと思います。さて、今回はGDALですね！GDALは以下記事のように当ブログで紹介はしているのですが、今回はShapeファイルの作成にチャレンジしてみたいと思います。

www.gis-py.com

sanvarie.hatenablog.com

また、今までのブログではArcPyでオープンストリートマップのデータ（.osm)をShapeに変換してみた - GIS奮闘記のようにArcpyを使用してのShapeファイルの作成方法は紹介したことがあるのですが、これだとArcGISに依存した作成方法なので、ArcGISを持っていない方でもShapeファイルの作成が手軽にできるようにGDALのようなフリーのライブラリでの作成方法を紹介します。

GDALとは

GDAL - Wikipedia
色々小難しいことが書いてありますが、要はフリーのGISライブラリです。

インストール

以下のようなサイトを参考にインストールをお願いします。
http://blog.godo-tys.jp/2013/06/19/2465/
https://github.com/domlysz/BlenderGIS/wiki/How-to-install-GDAL

実行環境

Windows7 64bit
ArcGIS10.4.1
ArcMap10.4.1
Python2.7.3

サンプルコード

ポイントShapeを作成するサンプルです。

# -*- coding: utf-8 -*-
import osgeo.ogr as ogr
import osgeo.osr as osr

# 今回作成するポイントの緯度経度
lat_lon = [[35.679933,139.714465],[35.1231433,139.2533465]]

# shapefileドライバ
driver = ogr.GetDriverByName("ESRI Shapefile")

# Shape出力先を設定
data_source = driver.CreateDataSource(ur"D:\gis-py\gdal\sample.shp")

# 投影法を設定（projectionファイルがうまく作成されないので、コメントアウト）
#srs = osr.SpatialReference()
#srs.ImportFromEPSG(4612)

# レイヤ作成
#layer = data_source.CreateLayer("sample", srs, ogr.wkbPoint) # projetionファイル作るときに使う
layer = data_source.CreateLayer("sample",  geom_type= ogr.wkbPoint)

# フィールド追加
field_name = ogr.FieldDefn("column", ogr.OFTString)
field_name.SetWidth(24)
layer.CreateField(field_name)

layer.CreateField(ogr.FieldDefn("Latitude", ogr.OFTReal))
layer.CreateField(ogr.FieldDefn("Longitude", ogr.OFTReal))

# Shape作成
for l in lat_lon:
  # フィーチャ作成
  feature = ogr.Feature(layer.GetLayerDefn())

  # 属性を設定
  feature.SetField("column", "test")
  feature.SetField("Latitude", l[0])
  feature.SetField("Longitude", l[1])

  # 作成するポイントに緯度経度をセット
  wkt = "POINT(%f %f)" %  (float(l[1]) , float(l[0]))

  # ポイント作成
  point = ogr.CreateGeometryFromWkt(wkt)
  feature.SetGeometry(point)
  layer.CreateFeature(feature)
  feature = None

data_source = None