🏕️
GIS × Python Tutorial 6.3 ~ DTMの後処理 ~
この記事は「GIS × Python Tutorial」の関連記事です。
はじめに
Session7 の記事では CSF による点群データ分類と DTM の作成を行いました。しかし作成した DTM は所々穴が開いていたり、その他に地表付近の地物が影響していたのか、ボコボコとした形になっていました。
この記事では、穴（欠損）を埋めしたり、データを滑らかにする為の後処理について解説します。
今回使用するデータ
今回は前回の記事で作成した DTM を使用します。実際に動かしてみたい方は、以下のURLから始めて下さい。
anaconda の環境設定ファイルや、この記事の作成に使用した Notebook は GitHub に保存しています。
インポート
cell_1
import os
from pprint import pprint
import japanize_matplotlib
from matplotlib import pyplot as plt
import matplotlib.patches
from matplotlib.spines import Spine
import numpy as np
import rasterio
import rasterio.fill
import rasterio.mask
import rasterio.plot
import scipy.ndimage
import shapely
from shapely.plotting import plot_polygon
plt.style.use('seaborn-v0_8-whitegrid')
japanize_matplotlib.japanize()
# 読み込みファイル
FILE_PATH = '../datasets/01ID7913_DTM_R05.tif'
FILE_PATH_RGB = '../datasets/01ID7913_RGB.tif'
IMAGE_DIR = '../images/session9'
1. Raster
1-1. Rasterの読み込みとメタデータの確認
Raster データの読み込みには rasterio を使用します。まずは DTM の範囲やサイズを確認してみましょう。
cell_2
def raster_metadata(dst: rasterio.DatasetReader) -> dict:
bounds = dst.bounds
scope = shapely.box(*bounds)
nodata_count = len(dst.read(1)[dst.read(1) == dst.nodata])
return {
'type': type(dst),
'epsg': dst.crs.to_epsg(),
'indexes': dst.indexes,
'data_type': dst.dtypes,
'cells': dst.width * dst.height,
'shape': dst.shape,
'height': dst.height,
'width': dst.width,
'bounds': dst.bounds,
'area': scope.area,
'x_length': abs(bounds.top - (bounds.bottom)),
'y_length': abs(bounds.left - (bounds.right)),
'resolution': dst.meta['transform'][0],
'stats': dst.statistics(1),
'nodata': dst.nodata,
'nodata_cells': nodata_count,
'nodata_ratio': nodata_count / (dst.width * dst.height),
'polygon': shapely.to_wkt(shapely.box(*dst.bounds._asdict().values())),
}
dst = rasterio.open(FILE_PATH)
metadata = raster_metadata(dst)
pprint(metadata, sort_dicts=False, underscore_numbers=True)
output_2
{'type': <class 'rasterio.io.DatasetReader'>,
'epsg': 6_669,
'indexes': (1,),
'data_type': ('float64',),
'cells': 3_003_501,
'shape': (1_501, 2_001),
'height': 1_501,
'width': 2_001,
'bounds': BoundingBox(left=-4000.0, bottom=37500.0, right=-2999.5, top=38250.5),
'area': 750875.25,
'x_length': 750.5,
'y_length': 1000.5,
'resolution': 0.5,
'stats': Statistics(min=143.08, max=516.9, mean=303.36050549932, std=67.100945441007),
'nodata': -9999.0,
'nodata_cells': 45_824,
'nodata_ratio': 0.015256861908819076,
'polygon': 'POLYGON ((-2999.5 37500, -2999.5 38250.5, -4000 38250.5, -4000 ', '37500, -2999.5 37500))'}
1-2. Raster の可視化
Raster は画像と同じく行列のデータなので、matplotlib で可視化する事が出来ます。
rasterio.open で読み込んだデータは rasterio.io.DatasetReader のオブジェクトになり、read メソッドで各配列にアクセスする事が可能です。今回読み込んだデータは DTM であり、地形の高さのデータのみしか保存されていないので配列の数は 1つですが RGB 画像の場合は 3つの配列が保存されています。ここで注意が必要なのが、1枚目の配列にアクセスするには indexes=0 ではなく indexes=1 でアクセスするという事です。感覚的には 0始まりと思いがちですが、上で表示したメタデータでも 'indexes: (1,)'となっているのが確認できます。全て一緒に読み込む場合 read() として indexes を指定しなければ、全ての配列を読み込む事が出来ます。
その他に、matplotlib.pyplot.imread で画像を読み込んだ場合、配列の形状は (高さ, 幅, 次元数) となりますが、rasterio.read で読み込んだ場合は (次元数, 高さ, 幅) となります。その為、普通に読み込んだ画像をそのまま matplotlib で可視化する事は出来ません。可視化するのであれば numpy.dstack で形状を変更してから可視化しましょう。
cell_3
# 'read'では'indexes'で確認した番号を指定します。
dtm = dst.read(indexes=1)
# 描画の範囲をmatplotlibのlimに適用させる
plt.figure(figsize=(7, 5))
extent = (dst.bounds.left, dst.bounds.right, dst.bounds.bottom, dst.bounds.top)
plt.imshow(dtm, extent=extent, cmap='terrain')
plt.colorbar(shrink=0.8)
# plt.savefig(os.path.join(IMAGE_DIR, 'original_dtm.png'), dpi=350)
plt.show()
欠損値（-9999.0）が邪魔して綺麗に可視化されていません。可視化する値の範囲に制限を掛けましょう。ついでに欠損セルだけは赤色で表示させます。
cell_4
fig = plt.figure(figsize=(8, 6))
# 欠損値を赤色で目立たせる
nodata = np.zeros(dtm.shape)
nodata[:, :] = 255
plt.imshow(nodata, cmap='bwr_r', extent=extent)
patches = [matplotlib.patches.Patch(color='red', label="nodata")]
plt.legend(handles=patches, bbox_to_anchor=(1.25, 1.1))
# DTMの可視化
dtm_nan = np.where(dtm < 0, np.nan, dtm)
_cb = plt.imshow(
X=dtm_nan,
cmap='terrain',
extent=extent
)
fig.colorbar(_cb, shrink=0.8)
plt.title('DTMの欠損値を可視化', fontsize=16)
plt.xlabel('X', fontsize=13)
plt.ylabel('Y', fontsize=13)
plt.grid(False)
# plt.savefig(os.path.join(IMAGE_DIR, 'nodata_dtm.png'), dpi=350)
plt.show()
1-3. rasterio による可視化
colorbar が不要な場合は rasterio.plot.show を使用する事で、欠損値の処理や範囲指定をする必要がないので楽です。
cell_5
fig, ax = plt.subplots(figsize=(7, 5))
rasterio.plot.show(dst, cmap='terrain', ax=ax)
ax.grid(False)
# plt.savefig(os.path.join(IMAGE_DIR, 'rasterio_plot_show.png'), dpi=350)
plt.show()
2. 穴埋め（欠損値補完）
rasterio には rasterio.fill.fillnodata というメソッドが用意されており、簡単に欠損値の処理を行う事ができます。これには 'mask'として、同じ形状の配列を渡す必要があります。mask 配列は "欠損値がある場所を 0 それ以外には 255" を入力したものが必要です。ただし、別に自分で作成しなくとも DatasetReader.read_masks で取得する事が可能です。
cell_6
# 欠損値の穴埋め
dtm_fill = rasterio.fill.fillnodata(
# DTMの行列をCopyして渡す
image=dtm.copy(),
mask=dst.read_masks(1),
# 補完する場合にどこまで離れたセルまで検索するか。デフォルトが100です。
max_search_distance=100,
)
# 可視化
fig, ax = plt.subplots(ncols=2, figsize=(14, 6))
arys = [dtm_nan, dtm_fill]
titles = ['処理前のDTM', '欠損値の処理を行ったDTM']
for _ax, ary, title in zip(ax, arys, titles):
_ax.imshow(nodata, cmap='bwr_r')
_ax.imshow(ary, cmap='terrain')
_ax.set_title(title, fontsize=16)
_ax.grid(False)
_ax.axes.xaxis.set_visible(False)
_ax.axes.yaxis.set_visible(False)
plt.subplots_adjust(wspace=0.05)
# plt.savefig(os.path.join(IMAGE_DIR, 'comparison_dtm.png'), dpi=350)
plt.show()
3. データの保存
欠損値を処理する事が出来たので、一度この配列を Raster として保存しましょう。
rasterio では普通に画像として保存する事も当然できますが、今の状況の様に Notebook でデータを見ながら実行している場合、いちいちデータを出力していてはディレクトリの中が汚くなってしまいます。rasterio ではインメモリに保存することも可能なので、今回はその方法を試してみます。
cell_7
# インメモリにRasterを作成する
dst_fill = rasterio.MemoryFile().open(**dst.meta)
# インメモリのRasterに穴埋めしたDTMを書き込む
dst_fill.write(dtm_fill, 1)
pprint(raster_metadata(dst_fill), sort_dicts=False, underscore_numbers=True)
output_7
{'type': <class 'rasterio.io.DatasetWriter'>,
'epsg': 6_669,
'indexes': (1,),
'data_type': ('float64',),
'cells': 3_003_501,
'shape': (1_501, 2_001),
'height': 1_501,
'width': 2_001,
'bounds': BoundingBox(left=-4000.0, bottom=37500.0, right=-2999.5, top=38250.5),
'area': 750875.25,
'x_length': 750.5,
'y_length': 1000.5,
'resolution': 0.5,
'stats': Statistics(min=143.0800018310547, max=516.9000244140625, mean=301.83231935691396, std=67.80678643427129),
'nodata': -9999.0,
'nodata_cells': 0,
'nodata_ratio': 0.0,
'polygon': 'POLYGON ((-2999.5 37500, -2999.5 38250.5, -4000 38250.5, -4000 ', '37500, -2999.5 37500))'}
4. DTMの平滑化
4-1. 拡大して表示する
DTMの平滑化を行う前に、何故わざわざ今ある情報を消すような事（平滑化）をするのか。実際に穴埋めした DTM を拡大して見てみましょう。
cell_8
# RGB画像を読み込む
dst_rgb = rasterio.open(FILE_PATH_RGB)
# 拡大する範囲指定
extended_scope = shapely.box(-3250, 37700, -3100, 37800)
# オリジナルのDTMを可視化
fig = plt.figure(figsize=(10, 7))
shape = (5, 6)
ax1 = plt.subplot2grid(shape, loc=(0, 0), rowspan=3, colspan=4)
ax1.imshow(dtm_fill, cmap='gray', extent=extent)
plot_polygon(
extended_scope,
ax=ax1,
color='red',
facecolor=(0, 0, 0, 0),
linewidth=2,
add_points=False
)
# 拡大したDTMを可視化
ax2 = plt.subplot2grid(shape, loc=(2, 0), rowspan=2, colspan=3)
extended_dtm = rasterio.mask.mask(dst_fill, shapes=[extended_scope], crop=True)[0][0]
ax2.imshow(extended_dtm, cmap='terrain')
# 拡大したRGB画像を可視化
ax3 = plt.subplot2grid(shape, loc=(2, 3), rowspan=2, colspan=3)
extended_rgb = rasterio.mask.mask(dst_rgb, shapes=[extended_scope], crop=True)[0]
ax3.imshow(np.dstack(extended_rgb), cmap='terrain')
# 各種設定
titles = ['Original DTM', 'Extended DTM', 'Extended RGB']
for i, (_ax, title) in enumerate(zip(fig.axes, titles)):
_ax.set_title(title, fontsize=16)
_ax.grid(False)
_ax.axes.xaxis.set_visible(False)
_ax.axes.yaxis.set_visible(False)
if 0 < i:
for children in _ax.get_children():
if isinstance(children, Spine):
children.set_color('red')
# plt.savefig(os.path.join(IMAGE_DIR, 'Extended_DTM.png'), dpi=350)
plt.show()
拡大すると丸い粒のようなものが沢山確認できます。本当にこの様な地形なのでしょうか?
RGB画像を見るとわかる様に、この拡大している場所は森林です。恐らくこれらは地表付近の下層植生（低木やシダ類など）が DTM の計算に含まれてしまったからだと推測できます。出来ればこれらを処理して取り除き、なるべく滑らかに、元の情報を消さずに新たな DTM を作成したいですね。
4-2. サンプルの確認
RasterData から座標を指定してサンプルを取得します。
cell_9
def crop_raster(dst, scope):
return rasterio.mask.mask(dst, shapes=[scope], crop=True)[0][0][0]
# サンプル範囲。
crop_scope = shapely.box(-3650, 38200, -3450, 38200.5)
fig = plt.figure(figsize=(7, 8))
# DTMの可視化
shape = (3, 1)
ax1 = plt.subplot2grid(shape, loc=(0, 0), rowspan=2)
ax1.imshow(dtm_fill, cmap='terrain', extent=extent)
plot_polygon(crop_scope, ax=ax1, color='red', add_points=False, label='sample')
ax1.grid(False)
ax1.set_title('サンプル範囲の確認', fontsize=16)
ax1.legend(bbox_to_anchor=(1, 1))
ax1.axes.xaxis.set_visible(False)
ax1.axes.yaxis.set_visible(False)
# サンプル範囲の可視化
ax2 = plt.subplot2grid(shape, loc=(2, 0), rowspan=1)
ori_cs = crop_raster(dst_fill, crop_scope)
ax2.plot(ori_cs, c='red', label='terrain');
ax2.set_title('断面図', fontsize=16)
# plt.savefig(os.path.join(IMAGE_DIR, 'crop_sample.png'), dpi=350)
plt.show()
4-3. 畳み込みで滑らかにする
DTM に対して畳み込みを行う事で、近傍のセルの平均を計算して滑らかな DTM を作成する事が出来ます。
対象とするセルを含む25個のセル（5 × 5）で畳み込みを行って結果を見てみましょう。
cell_10
def convolution_kernel(distance):
if distance % 2 == 0:
distance += 1
shape = (distance, distance)
cells = distance * distance
return np.ones(shape) / cells
# 25個のセルで平滑化
rows = 5
kernel = convolution_kernel(rows)
dtm_conv = scipy.ndimage.convolve(dtm_fill, kernel, mode='nearest')
# インメモリに書き込み
dst_conv = rasterio.MemoryFile().open(**dst.meta)
dst_conv.write(dtm_conv, 1)
# 可視化範囲の切り抜き
crop_conv = crop_raster(dst_conv, crop_scope)
cell_11
def plot_cross_section(original_ary1d, crop_ary1d, color, title, name, outpath=None):
fig = plt.figure(figsize=(12, 6))
shape = (3, 1)
# 横断図の可視化
ax1 = plt.subplot2grid(shape, loc=(0, 0), rowspan=2)
idx = np.arange(0, len(original_ary1d))
ax1.plot(original_ary1d, lw=1, c='black', zorder=1, label='original')
ax1.plot(crop_ary1d, c=color, lw=1, zorder=3, label=name)
ax1.set_title(title, fontsize=17)
# オリジナルデータと平滑化後の差分を可視化
ax2 = plt.subplot2grid(shape, loc=(2, 0))
ax2.bar(idx, crop_ary1d - original_ary1d, lw=0.3,
facecolor=color, edgecolor='black', label='diff')
ax2.set_ylim(-3, 3)
# 諸々
for _ax, ylabel in zip(fig.axes, ['terrain', 'difference']):
_ax.set_ylabel(ylabel)
_ax.legend()
if outpath:
plt.savefig(outpath, dpi=350)
plt.show()
plot_cross_section(ori_cs, crop_conv, 'red',
'DTMの断面図と畳み込みによる平滑化',
'convolution')
4-4. ガウシアンフィルターを使用した平滑化
ガウシアンフィルターは、画像処理において平滑化を目的として使用されるフィルタの一種です。ガウシアンフィルターは対象とするセルとの距離に応じて重みづけを行う為、自然な平滑化を行う事が出来ます。
cell_12
# ガウシアンフィルターを使用した平滑化
dtm_gauss = scipy.ndimage.gaussian_filter(dtm_fill, sigma=2, mode='nearest')
# インメモリに書き込み
dst_gauss = rasterio.MemoryFile().open(**dst.meta)
dst_gauss.write(dtm_gauss, 1)
# 可視化範囲の切り抜き
crop_gauss = crop_raster(dst_gauss, crop_scope)
plot_cross_section(ori_cs, crop_gauss, 'green',
'DTMの断面図とガウシアンフィルターによる平滑化',
'gaussian')
4-5. 畳み込みの繰り返し
当然ですがただ単に滑らかにしたいだけであれば、畳み込みを繰り返すだけで滑らかなデータを作成できます。しかし滑らかにしすぎると、窪地などの地形の情報が消えてしまうので注意が必要です。
cell_13
def multi_conv(ary_2d, start_distance, count):
for _ in range(count):
kernel = convolution_kernel(start_distance)
ary_2d = scipy.ndimage.convolve(ary_2d, kernel, mode='nearest')
start_distance += 2
return ary_2d
# 複数回の畳み込みによる平滑化
dtm_conv_m = multi_conv(dtm_fill, 3, 5)
# インメモリに書き込み
dst_conv_m = rasterio.MemoryFile().open(**dst.meta)
dst_conv_m.write(dtm_conv_m, 1)
# 可視化範囲の切り抜き
crop_conv_m = crop_raster(dst_conv_m, crop_scope)
plot_cross_section(ori_cs, crop_conv_m, 'deeppink',
'DTMの断面図と複数回の畳み込みによる平滑化',
'multi convolution')
4-6. 最小値フィルターを組み合わせた平滑化
基本的にはガウシアンフィルターを一度通すだけで、データが滑らかになるのでそれで問題ないと思います。しかし下層植生が非常に豊かな場合そう単純ではありません。
豊かな森林であれば、そもそも地表にレーザーパルスが届く事が非常に少ないです。
-
背の高い樹木
1-1. 葉
1-2. 枝
1-3. 幹
-
背の低い樹木
2-1. 葉
2-2. 枝
2-3. 幹
-
シダ類など地表に生えている植物
樹木の密度によっても変わりますが、樹木が密であれば樹冠同士が近く、レーザーパルスが届きづらくなります。また、樹木が疎であれば日光が当たる為に下層植生が豊かになります。落葉広葉樹林であれば、晩秋に測るのが一番簡単に地表を捉えられるでしょうが、必ずしもその様に計測を実行できるとは限りません。
ここではガウシアンフィルターと最小値フィルターを組み合わせて、DTM を滑らかにしつつも、オリジナル DTM の下に潜るような計算を行います。点群データがあまり地表部分に届いていないと判断した場合は、この方法を使用する事で、仮想的に新たな地表を作成する事が出来ます。
cell_14
# ガウシアンフィルターによる平滑化
dtm_gauss = scipy.ndimage.gaussian_filter(dtm_fill, sigma=2, mode='nearest')
# 最小値フィルターで近傍の最小値を取得
rows = 3
dtm_min = scipy.ndimage.minimum_filter(dtm_fill, rows, mode='nearest')
dtm_min = scipy.ndimage.minimum_filter(dtm_min, rows + 2, mode='nearest')
# 上記の情報を使用して
dtm_under = (dtm_gauss + dtm_min) / 2
dtm_under = scipy\
.ndimage\
.gaussian_filter(dtm_under, sigma=2, mode='nearest', truncate=3)
# インメモリに書き込み
dst_under = rasterio.MemoryFile().open(**dst.meta)
dst_under.write(dtm_under, 1)
# 可視化範囲の切り抜き
crop_under = crop_raster(dst_under, crop_scope)
plot_cross_section(ori_cs, crop_under, 'blue',
'DTMの断面図と最小値フィルターを組み合わせた平滑化',
'under convolution')
4-7. 各平滑化結果を同時に可視化する
計算したものを合わせて見てみましょう。
cell_15
idxs = slice(50, 150)
plt.figure(figsize=(12, 5))
plt.plot(ori_cs[idxs], color='black', lw=1, label='original')
plt.plot(crop_conv[idxs], color='red', lw=1, label='convolution')
plt.plot(crop_gauss[idxs], color='green', lw=1, label='gaussian')
plt.plot(crop_conv_m[idxs], color='deeppink', lw=1, label='multi convolution')
plt.plot(crop_under[idxs], color='blue', lw=1, label='under convolution')
plt.legend()
plt.title('DTMの断面図と平滑化手法の違いによる差①', fontsize=17)
plt.ylabel('terrain', fontsize=13)
plt.yticks([])
plt.xticks([])
# plt.savefig(os.path.join(IMAGE_DIR, 'cs_comparison_1.png'), dpi=350)
plt.show()
4-8. 各DTMを表示する
最後に拡大範囲の DTM を可視化してどの様に変化しているかを見てみましょう。
cell_15
datasets = [dst_fill, dst_conv, dst_gauss,
dst_conv_m, dst_under]
titles = ['Original', 'Convolution', 'Gaussian filter',
'Multiple convolution', 'Underground convolution']
fig = plt.figure(figsize=(9, 10))
for i, (_dst, title) in enumerate(zip(datasets, titles)):
ax = fig.add_subplot(2, 3, i + 1)
dst_croped = rasterio.mask.mask(_dst, shapes=[extended_scope], crop=True)[0][0]
ax.imshow(dst_croped, cmap='terrain')
ary1d = dst_croped.flatten()
level = np.arange(ary1d.min(), ary1d.max(), 2)
ax.contour(dst_croped, levels=level, linewidths=0.3, colors='black')
ax.axes.yaxis.set_visible(False)
ax.axes.xaxis.set_visible(False)
ax.grid(False)
ax.set_title(title)
plt.subplots_adjust(wspace=0.05, hspace=0.15)
# plt.savefig(os.path.join(IMAGE_DIR, 'pp_dtm.png'), dpi=350)
plt.show()
cell_16
# データの保存
with rasterio.open(r'../datasets/01ID7913_DTM_R05_SMOOTH.tif', 'w', **dst.meta) as _dst:
_dst.write(np.round(dtm_under, 3), 1)
おわりに
ここまでで、DTM の穴埋めや平滑化について紹介しました。rasterio も結構使用しているので、Raster の扱いも大分慣れてきたかと思います。次回は今回作成した DTM と点群データを使用して、DSMとCHMの作成を行います。
参考
Notebookの置いてあるリポジトリ
畳み込み
rasterio.plot
rasterio.fill
rasterio.mask
ガウシアンフィルター
scipy.ndimage.
Discussion