PythonでGDALを使い倒そう！

さて、本日はGDALについて書いてみようと思います。GDLについては以下のように以前に何個か記事を書いているのですが、非常に奥が深いライブラリですので紹介していない機能がまだまだたくさん残っています。今回はラスターの扱いについて紹介します。

www.gis-py.com

GDALとは

公式サイトによると、GDALはフリーのGISライブラリでラスターとベクターを扱うことができます。よくGDALはラスター用のライブラリと書かれているサイトを見るのですが、ベクターも扱うことができます。GDAL — GDAL documentation

サンプルデータ

今回はラスターの読み込みとラスターの作成を行います。以下二つのデータを用意してください。

ラスターの読み込み

以下サイトで航空写真がダウンロードできますので、今回はそのうちの一つを使ってみようと思います。
空中写真オープンデータ

航空写真とは？という方は以下エントリーを参考にしてみてください。
www.gis-py.com

ラスターの作成

ShapeファイルをGEOTIFFに変換してみます。以下エントリーで使用した日本地図をGEOTIFFに変換します。

www.gis-py.com

サンプルコード

ラスターの読み込み

from osgeo import gdal
raster = gdal.Open( r'D:\gis-py\library\gdal\data\19640721c-wide-niigata-eq-0-1.tif')

print raster.GetProjection() #投影法
print raster.RasterXSize      #Xサイズ
print raster.RasterYSize      #Yサイズ
print raster.RasterCount    #バンド数
print raster.GetMetadata() #メタデータ
print raster.GetDriver().LongName

geotransform= raster.GetGeoTransform()
print geotransform[0] #左上X
print geotransform[3] #左上Y
print geotransform[1] #セルサイズX
print geotransform[5] #セルサイズY

このようにラスターのプロパティを取得することができます。
f:id:sanvarie:20181029161812p:plain

ラスターの作成

from osgeo import gdal, ogr

# ピクセルサイズを設定
pixel_size = 0.01
NoData_value = -9999

# Shapeファイルを指定
vector_fn = ur'D:\gis-py\library\gdal\data\japan_ver80.shp'

# 作成するTIFFを指定
raster_fn = r'D:\gis-py\library\gdal\outcome\japan_geo.tif'

# Shapeをオープン
source_ds = ogr.Open(vector_fn)
source_layer = source_ds.GetLayer()
x_min, x_max, y_min, y_max = source_layer.GetExtent()

# TIFF作成のための設定
x_res = int((x_max - x_min) / pixel_size)
y_res = int((y_max - y_min) / pixel_size)
target_ds = gdal.GetDriverByName('GTiff').Create(raster_fn, x_res, y_res, 1, gdal.GDT_Int16)
target_ds.SetGeoTransform((x_min, pixel_size, 0, y_max, 0, -pixel_size))
band = target_ds.GetRasterBand(1)
band.SetNoDataValue(NoData_value)

# ラスタ化
gdal.RasterizeLayer(target_ds, [1], source_layer, burn_values=[0])