2022年7月31日日曜日

斜面の向きと地震波

地形量と崩壊の関係を整理してみようかと考えていました。

DEM が 1枚あれば、多くの地形量を作ることができます。それと相関性を見てみるだけです。GIS 仕事です。

文献を探してみると、地震時でトライされた内容がありました。誘因として雨の要因を考慮する必要がなくなるため、地震時の方が良いのでしょう。

安田ほか(2006)動的振動解析による地震時の加速度応答および斜面変位と地形効果に関する考察
https://doi.org/10.11475/sabo1973.59.4_3

統計的に整理するに数が少ないと思われますが、方針はこれと同様です。

地震時だと揺れの方向と斜面の方向が関連するので、数値実験の方が整理しやすいでしょう。斜面のモデルは少なくても、揺れの方向を自由に変えてその応答を見れば良いだけですので理屈は簡単です。時間はかかりそうですが。

実際の地震波では、斜面の法線方向から下向きの場合が最も危険になるでしょう。以下では、ある程度考慮されているようです。下向きの角度はよくわかりませんでしたが。

篠田ほか(2018)広域的な地震時斜面崩壊危険度図の作成方法
https://doi.org/10.2208/jscejge.74.177

どなたか整理されていないでしょうか。単純な話なのでありそうですが。


2022年7月30日土曜日

六甲式の改良

地震時崩壊を調べていると、国土地理院の方々が書かれた以下の文献を見つけました。

神谷ほか「地震による斜面崩壊危険度評価判別式「六甲式」の改良と実時間運用」写真測量とリモートセンシング 51 (6), 381-386, 2013

F=0.075s-8.92c+0.006a-0.3228 (六甲式)
G=4.38log(s-119c)+3.93loga-15.27(修正六甲式)

例えば曲率が0の場合(地表面の凹凸がない場合),傾斜が43度以上では  PGA が 0gal であっても崩壊と予想する。また,PGA が583gal(計測震度5.8に相当)以上では,傾斜が0度であっても崩壊すると予想する。本システムでは,このような場合でも,妥当な結果を得る必要がある。
言われてみると、納得。閾値を調整すれば六甲式でもカバーできそうですが。
六甲式,修正六甲式は,PGAに関係する項と,関係しない項の和である。PGAに関係しない項を,それぞれ「部分六甲式」,「部分修正六甲式」と呼ぶ。あらかじめ,小グリッドの全てのセルにおける部分六甲式の値を計算し,大グリッドの全てのセルにおいて,「大グリッドのセルに含まれる小グリッドの部分六甲式の値」のヒストグラムを作成する。地震発生時には,大グリッドの部分六甲式のヒストグラムと各階級の代表値を用いて,六甲式等のヒストグラムを計算する。この計算は,大グリッド単位で実行するため,小グリッド単位の計算と比較し,大幅に計算量を削減することができる。
あらためて10mメッシュの傾斜角や曲率を計算し表示してみましたが、2億メッシュでも数十秒~数分程度でした。うまい工夫だと思いますが、今では10mメッシュのままでも計算できるでしょう。日本全国の10mメッシュで s, c の項を計算しておけば、地震時に PGA や震度情報を取り込むことですぐに予測結果を得られます。

今後、高度化を目指すなら、まずは5mDEMの併用でしょうか。角度や曲率を計算する際に海岸線沿いの崖などが考慮できない点をカバーできると思われます。
シンプルかつある程度検証されてきた式ですので、地震のたびに検証・更新される方がいらっしゃっても良いように思います。一番手間の書かかるポリゴン制作も衛星画像等から精度よく自動作成できるようになれば、一気に実現しそうです。が、まだ先ですかね。

 

2022年7月29日金曜日

Pandas DataFrame の不一致

Pandas で EXCELから読み込んだ内容と、csv に変換してから読み込んだ内容が一致せず、どこが原因かを探っていました。

最初はエンコードを疑いましたが、違いました。
原因は数字の表記でした。

 

EXCEL上で 100.0 と表記していると、csv 経由では float、excel からは int として読み込まれる仕様でした。

 

天然ダムの構成

先日の引用文献の Fanさん、地質の方でしょうか。天然ダムの構成に関する話題も報告されていました。

Fan et al. (2017.) Characteristics and classification of landslide dams associated with the 2008 Wenchuan earthquake.

Sub-type III: landslide dams with partly intact rock strata at the base topped by large boulders and blocks or soil with rock fragments, showing a two-layered (e.g. the Xiaojiaqiao dam) or three layered (the Tangjiashan dam, Fan et al. 2012b) internal structure. Such dams are mainly formed by deep-seated rock slides or avalanches that fail as a largely intact mass and run into relatively narrow valleys. Due to the topographic constraints on the landslide movement, the sliding rock mass cannot be disintegrated completely, part of which still keeps the original geological
structure.


2022年7月26日火曜日

Seismic vs hydrologic triggering of landslide dams

USGS の Landslide Hazards Seminar を発見。
このような動画を公開されているとは。さすがです。

面白そうな動画を一つ見てみました。
Seismic vs hydrologic triggering of landslide dams in the Oregon Coast Range
https://www.usgs.gov/media/videos/seismic-vs-hydrologic-triggering-landslide-dams-oregon-coast-range


文献としては以下です。
Struble et al. (2020) Dendrochronological dating of landslides in western Oregon: Searching for signals of the Cascadia A.D. 1700 earthquake
https://doi.org/10.1130/B35269.1
LaHusen et al. (2020) Rainfall triggers more deep-seated landslides than Cascadia earthquakes in the Oregon Coast Range, USA
https://www.science.org/doi/10.1126/sciadv.aba6790

Surface morphology: 地形の起伏と年代の関係は、見たことがあります。直接年代を求めることができませんので、14Cなどの結果との対比が必須になります。実務では見かけませんね。

Dendrochronology:  年輪年代学。木の年輪で降水量の多寡を把握できることは知りませんでした(13:50頃~)。年輪の個々の幅を計るのは地道な作業だと思いますが、文献では専用のソフトを利用したと書かれています。検索すると動画がありました。色の急変箇所でわかるのでしょうね。面白い。
長期成分をスプラインで除去する方法は文献には書かれていませんが、見た目以上に繊細な作業です。スプラインの頂点の間隔で答えが変わるのですが、どの木でも同じ間隔にしておけば問題ないのでしょうか。あとは相互相関を取ってマッチングするだけですので理屈だけは簡単なのですが、実装は簡単ではなかったでしょう。地震と同じ手法がこのようなところでも使えるのは驚き。ま、ベクトルにしてしまえば同じことですね。

天然ダムの安定性に関するお話も、シンプルで個人的には好きです(31:50頃~)。日本だと、異なるのでしょうか。関連文献は以下です。
Fan et al. (2020) The formation and impact of landslide dams – State of the art https://doi.org/10.1016/j.earscirev.2020.103116
Oliver Korup (2004) Geomorphometric characteristics of New Zealand landslide dams https://doi.org/10.1016/j.enggeo.2003.11.003

日本を調べている方は、またもや国外の方。怖いですね。

結論としては、地震よりも雨がトリガーになっていたようです。地道な研究を積み重ね、このような壮大な結論を導かれた方々にはリスペクトしかありません。

2022年7月23日土曜日

GeoPandas の投影変換

GeoPandas にて CRS の機能しないPCがありました。
CRSError: Invalid projection: EPSG:4326 #1887

思い当たる点はありました。PostgreSQL です。
CRS 以外は機能するので当面それを避けて凌ぐことに。

今日、別の PC で仮想環境を作りましたが、こちらも失敗。
Error when importing geopandas #1756

午前中で解決せず、気を取り直して午後も作業。最終的には既存の環境に GeoPandas 等の不足しているライブラリを足して解決です。15時までかかりました。

で、実施したかった投影変換。
gdf2000=gdf2011.to_crs(epsg=4612)
gdf2011['geometry'].equals(gdf2000['geometry'])
#True

なぜ?東北地方も含まれているのですが。

epsg =4326 (WGS84) → True
epsg =4301 (TOKYO) →False

大まかにあっている測地系では変更しない、ということでしょうか。微妙な挙動です。

結局、ArcGIS で処理することにしました。ま、このような挙動を知ることができたので、良しとしましょう。

 

2022年7月21日木曜日

GPUと単精度

プロとお話をしていて、GPUの話題になりました。

プロ曰く、陽解法は単精度でもOKとのこと。

調べてみると、ありました。

GPUを用いた粒子法計算の高速化と大規模化
https://www.jstage.jst.go.jp/article/conf/21/0/21_95/_article/-char/ja/

倍精度と単精度で波形が全く異なる経験をしていましたが、よく考えてみると微動でした。強震だと概ね同程度の挙動を示すのでしょうか?

頭においておきましょう。

2022年7月18日月曜日

Arias Intensity

読み飛ばしてい地震動の指標について整理。Iaを知りませんでした。

地震動強さの指標 - PukiWiki (tmu.ac.jp)

\begin{eqnarray}PGA=\sqrt[]{a_X^2(t)+a_Y^2(t)+a_Z^2(t)}_{max}  (m/s^2)\end{eqnarray}

PGA はある時点の加速度を表しているに過ぎず、継続時間を考慮した地震エネルギー規模を表現できません。そこでアリアス強度。

Arias Intensity, IA (Arias, A. (1970) A measure of earthquake intensity., Seismic 692 Design for Nuclear Power Plants., pp. 438-483. 

\begin{eqnarray}I_a=\frac{\pi}{2g}\int_{t_1}^{t_2}[a(t)]^2 dt (m/s)\end{eqnarray}

加速度二乗平均値(加速度RMS)も Ia と同様に、時間方向の積分値を使用しています。

\begin{eqnarray}a_{rms}=\sqrt[]{\frac{1}{t_2-t_1}\int_{t_1}^{t_2}[a(t)]^2 dt} (m/s^2)\end{eqnarray}

2022年7月17日日曜日

コルモゴロフ・スミルノフ検定

2群のデータを比較するためにヒストグラムを作成しました。
それを眺めながら、分布形状が似ている、似ていないを判断したいのですが、主観だと説明性に欠けます。ということで統計的検定の出番です。

以前、2標本(対応ナシ)の条件で Mann-Whitney's U test を使用しました。
https://phreeqc.blogspot.com/2020/12/mann-whitneys-u-test.html

今回は2群のサンプル数が異なるとともに、サンプルが非常に多く何らかの確率分布を仮定できそうでした。
調べてみると、ありますね。
・2標本のコルモゴロフ・スミルノフ検定
two-sample Kolmogorov-Smirnov(K-S) test

手順は、①2群を同じ階級幅で区分して累積度数を求める。②累積度数を累積確率に変換し、その差の絶対値を計算。③絶対値の最大値を使って検定統計量を求め、④棄却限界値と比較し判定する。棄却限界値は自由度f=2のχ2分布から求めます。

これは、好都合。ヒストグラムを全面積で正規化し、縦軸を確率密度で記載していたので処理は半分終わっています。あとは累積確率に直して差の最大値を求めるだけのようなもの。χ2分布の適用も、おかしくはない分布ですから、これで決定です。というか、2群のヒストグラムを書く際に、ルーチンワークとして求めておけば良い内容です。Python コードを残しておきましょう。

def hist_plot(col,x1,x2,label,bins):
    fig = plt.figure(figsize=(15,5))
    ax = fig.add_subplot(1,1,1)

    freq,bin_list,patch=ax.hist([x1,x2], density=True, alpha=0.5 ,edgecolor='black', bins=bins, label=label)
    ax.set_xlabel(col,fontsize=14)
    ax.set_ylabel('確率密度',fontsize=14)
    ax.grid(which='major',color='gray',linestyle='-')
    ax.grid(which='minor',color='gray',linestyle='--')
    ax.legend(borderpad = 0.8,framealpha=1.0,shadow=True)
   
    #2標本コルモゴロフ・スミルノフ検定(χ2分布、各標本数40以上、対応の無い2標本)
    df_p=pd.DataFrame({col:bin_list[1:],label[0]+'_確率密度':freq[0],label[1]+'_確率密度':freq[1]})
    df_p_sum=df_p.cumsum()*bin_list[1]
    df_p_sum=df_p_sum.drop(col,axis=1)
    df_p_sum.columns  = [label[0]+'_累積確率', label[1]+'_累積確率']
    df_p_sum['D']=abs(df_p_sum[label[0]+'_累積確率']-df_p_sum[label[1]+'_累積確率'])
    df_p=pd.concat([df_p, df_p_sum], axis=1)

    n1=len(x1)
    n2=len(x2)
    Dmax=df_p['D'].max()
    Dmax2=Dmax*Dmax
    T=4*Dmax2*(n1*n2)/(n1+n2)
    if T<5.99146:
        comment='有意水準5%で、差があるとは言えない。'
    else:
        comment='有意水準5%で、差があると言える。'       
    display(df_freq,'Dmax',Dmax,'T',T,comment)

Dが小さいと、差があるとは言えない(H0を棄却しない)

H0:帰無仮説(差がない)
H1:対立仮説(差がある)

2022年7月16日土曜日

Civil3D2022 サーフェス作成時のポイント数

他支店の方から、Civil3Dでの等高線の作り方を教えてほしいと連絡がありました。

いつもはV-nas Clairで作成されているようですが、非常に遅いので困っているとのこと。
すぐにできると思いましたが、沼にハマりました。

取り込むポイントは4800万点。これから等高線を作ってV-nasで表示したいとのこと。聞けば、その方も社内で偉い人?から投げられた作業のようでした。技術者なら必要な範囲を判断し、指示するはずですが。センスのない話です。

最初にハマったのは、Civil3D でサーフェスを保存できない現象。原因はメモリかと思いましたが余裕があるので異なりました。少しずつ点数を減らして限界を把握しました。結果、今回は4200万点まで保存できました。
調べてみると同様の報告がありました。サポートさんも御存じないようでしたが、以前から存在する現象のようです。現行の2023でも解決していないのかな。

次にハマったのが V-nas Clair の制約。Civil3D で作成した等高線を LandKit でインポートすると、「要素の制限値を超えた【ポリライン】があります。分解もしくは削除されました」というダイアログが出ます。サポートさんに伺うと、1本のポリラインに2^15以上の要素が含まれるとポリラインは分解、スプラインは削除する仕様とのことでした。分解すると起終点を閉じたポリラインは開くので、その間の要素が消えたように見えます。が、Civil の等高線は開いたポリラインのため、影響なし。よかった。その前に TrendPoint のスプライン等高線を読んでいたのですが、一部削除された原因がコレでした。

単純作業だったのですが、ハマりポイントが複数あったため時間がかかりました。ま、各ソフトウェアの限界を知ることができたので、良しとしましょう。

 

Software ポイントファイル数 ポイント総数 サーフェス作成 サーフェス表示 等高線作成 save-reopen V-nasClairでの等高線Opren 結果
Civil3D
2022.1.3
RCS:1 48,000,000 ×:フィルタなしで15分程度。
作成完了した時点でフリーズ。
https://knowledge.autodesk.com
境界のみ - ※ reopen時に再度サーフェスを作成しているような挙動。
数分かかる。
kome  ×サーフェスも等高線も作成できない。
TXT:1 48,000,000 〇:Civil15分、メモリ20GB コンター 〇サーフェスから取り出し ×:サーフェスを保存できない。
等高線のみなら可。
https://forums.autodesk.com
〇等高線のみ保存可
TXT:1 48,000,000 境界のみ ×  
TXT:2 48,000,000 境界のみ ×  
TXT:1 42,000,000 境界のみ
TXT:1 36,000,000 境界のみ  
TXT:1 24,000,000 境界のみ  
V-nasClair
2021.4
TXT:1 48,000,000 〇点群読み込み2分、TIN作成6分:メモリ18GB
※動かす度に再描画するため遅い。
TIN ×:等高線間隔設定まで6分。
その後一晩以上でも作成されず中止。
メモり100GB
〇:サーフェスreopen8分 - ×遅い。メモリーを使う。
等高線を作成できない
TrendPoint
Ver.9(9006)
TXT:1 48,000,000 ※時間を要するため中止 - 〇:作成に一晩。メモリ33GB 〇点群、等高線のみ。Save20分、reopen10分以上かかる。 ×:スプラインに完全対応していない。AutoCAD等で少しづつポリラインに変換する必要あり。 △遅い。

2022年7月3日日曜日

QGIS + PostgreSQL

ネット注文がメインになったことやコロナ禍ということで、大型書店に通う頻度が低くなりました。久しぶりに通うと、気になる図書がちらほら。この週末は2冊の 図書を読んで過ごしました。

そのうちの1冊がコチラ↓
愛知大学三遠南信地域連携研究センター「地域研究のための空間データ分析入門 -QGISとPostGISを用いて‐」

GIS は「習うより慣れろ」でしたので、この種の図書を読んだことがありません(ArcHydro 関連の図書は購入しましたが)。この図書には PostgreSQL との連携が書かれています。どのような場合に RDB を利用すべきか常々考えていましたので、思わず購入しました。

QGIS から PostgreSQL への接続は、何度か実施しています。が、実務レベルで採用したことはありません。SQL ではなく、Pythonで整形してしまってから GIS へ取り込むか、そのまま Python で表示してしまう流れの方が速かったからでしょう。紹介されていた作業内容も、データ加工の部分では C 等で作成されたライブラリで処理する Python の方が速いのかもしれません。

図書では PostgreSQL に shp を取り込み、PostGIS で空間分析まで実施されていました。このようなことができるとは知りませんでした。初めて PostGIS の役割を理解できたような気がします。Python では空間分析を GIS に任せていました。SQLで書いておけば、データ更新時の繰り返し処理の手間を省けるところが利点でしょうか。

単年で完結する作業ではデータベースを更新する必要がないため、RDB まで必要としません。一方、データを随時更新・追加したり、多くの方が データ入力に携わったりする状況では、RDB の方が良いのでしょう。
図書を読むことで、ようやく答えの1つを得ました。