import numpy as np
import pandas as pd
# 図化用のモジュール
import holoviews as hv
hv.extension('bokeh', logo=False)
とある地形データです。Lが平面距離(m)、Zが標高(m)です。Lは0.5m間隔に正規化されています。
df = pd.read_csv('test_bedElevation.csv', index_col=0)
df
L | Z | |
---|---|---|
0 | 0.0 | -5.06 |
1 | 0.5 | -5.07 |
2 | 1.0 | -4.94 |
3 | 1.5 | -4.94 |
4 | 2.0 | -4.89 |
... | ... | ... |
4359 | 2179.5 | -5.49 |
4360 | 2180.0 | -5.49 |
4361 | 2180.5 | -5.68 |
4362 | 2181.0 | -5.68 |
4363 | 2181.5 | -5.68 |
4364 rows × 2 columns
L, Z = df.L.values, df.Z.values
g = hv.Curve((df.L,df.Z)).options(xlabel='L[m]', ylabel='Z[m]',width=500)
g
複雑なデータですが、地形の表面に出現する細かい波長のデータ、つまり、高周波成分の特徴を分析します。
numpyのFFTを使います。データの個数とかも気にせず上手に変換してくれます。細かい中身は見てません。
dx = 0.5
N = len(Z)
freq = np.fft.fftfreq(N, d=dx)
F = np.fft.fft(Z)
g = hv.Curve((freq[1:int(len(F)/2)], np.abs(F/N*2)[1:int(len(F)/2)])).options(xlabel='波長の逆数[1/m]', ylabel='振幅[m]',width=500)
g.redim.range(x=(0,0.2))
高周波成分のみを使うため、ハイパスフィルターを使います。上記のFFTの結果の低周波成分を0として、フーリエ逆変換を行います。低周波成分の分類は波長が50m以上(逆数で0.02)としました。
F2 = np.copy(F)
ind = (freq < 0.02) & (freq > -0.02)
F2[ind] = 0.0
Z2 = np.fft.ifft(F2)
g = hv.Curve((freq[1:int(len(F2)/2)], np.abs(F2/N*2)[1:int(len(F2)/2)])).options(xlabel='波長の逆数[1/m]', ylabel='振幅[m]',width=500)
g.redim.range(x=(0,0.2))
g = hv.Curve((1/freq[1:int(len(F2)/2)], np.abs(F2/N*2)[1:int(len(F2)/2)])).options(xlabel='波長[m]', ylabel='振幅[m]',width=500)
g.redim.range(x=(0,50))
g = hv.Curve((L,Z2.real)).options(xlabel='L[m]', ylabel='Z[m]',width=500)
g
いい感じ高周波成分のみ抽出できました。
こんな感じでみると、0~900m間は、0.2-0.4m程度の波高の成分、900m以降は、0.1-0.2m程度の成分が主となっています。
場所によって卓越する成分が異なっていそうな感じです。
上記を参照に2区間に分けて、同じ処理をしてみます。
def highpass(Z, freq):
F = np.fft.fft(Z)
F2 = np.copy(F)
ind = (freq < 0.02) & (freq > -0.02)
F2[ind] = 0.0
Z2 = np.fft.ifft(F2)
return Z2, F2
Z1 = Z[:1800]
Z2 = Z[1800:]
dx = 0.5
freq1 = np.fft.fftfreq(len(Z1), d=dx)
Z1h, F1h = highpass(Z1, freq1)
freq2 = np.fft.fftfreq(len(Z2), d=dx)
Z2h, F2h = highpass(Z2, freq2)
卓越する波長の差は無さそうです。
g1 = hv.Curve((1/freq1[1:int(len(Z1)/2)], np.abs(F1h/len(Z1)*2)[1:int(len(Z1)/2)]), label = 's1')
g2 = hv.Curve((1/freq2[1:int(len(Z2)/2)], np.abs(F2h/len(Z2)*2)[1:int(len(Z2)/2)]), label = 's2')
g1*g2.redim.range(x=(0,50)).options(xlabel='波長[m]', ylabel='振幅[m]',width=500)
g1 = hv.Curve((np.array(range(len(Z1h)))*dx,Z1h.real), label = 's1').options(xlabel='L[m]', ylabel='Z[m]',width=500)
g2 = hv.Curve((np.array(range(len(Z2h)))*dx,Z2h.real), label = 's2').options(xlabel='L[m]', ylabel='Z[m]',width=500)
g1 * g2