AVS30（5次メッシュ）を 3次メッシュ内で平均化したいなあと思案。
全国だと EXCEL では行数が足りないので、正攻法の ArcGISで処理。
時間がかかるので Python でも追っかけ処理。
shpを読んだ後の主な処理はこれだけ。
#AVS30=0(水域)を削除
gdf=gdf[gdf['5th_CODE'] != 0]
#3次メッシュコード作成
gdf['3rd_CODE']=gdf['5th_CODE'].astype(np.int64)//100
#3次メッシュにおける平均を求める
df=gdf.groupby('3rd_CODE').mean()
空間演算が必要ないので、Python の方が早く終わりました。
便利。
任意関数の最小二乗法は EXCEL でも Pythonでも可能です。
精度、速度、記載の単純さで scipy.optimize.curve_fit の圧勝です。初期値によっては一発で求まらないこともありますが、 EXCELソルバーに比べると安定かつ高速で解を求めてくれます。変分推論は難しいですね。収束させるのにあたりのつけ方が必要になるようです。で、結局は最小二乗法に戻ってしまったり。慣れが必要なのでしょう。
Pandas の DataFrame でよく使う便利な乗算。
間違えやすいので、( ..)φメモメモ
**************************************
import numpy as np
import pandas as pd
df1=pd.DataFrame(np.ones(9).reshape(3, 3),
columns=['a','b','c'],
index=[1,2,3])
df1
| a | b | c | |
|---|---|---|---|
| 1 | 1.0 | 1.0 | 1.0 |
| 2 | 1.0 | 1.0 | 1.0 |
| 3 | 1.0 | 1.0 | 1.0 |
df2=pd.DataFrame(np.arange(9).reshape(3, 3),
columns=['a','b','c'],
index=[1,2,3])
df2
| a | b | c | |
|---|---|---|---|
| 1 | 0 | 1 | 2 |
| 2 | 3 | 4 | 5 |
| 3 | 6 | 7 | 8 |
#アダマール積
df1*df2
| a | b | c | |
|---|---|---|---|
| 1 | 0.0 | 1.0 | 2.0 |
| 2 | 3.0 | 4.0 | 5.0 |
| 3 | 6.0 | 7.0 | 8.0 |
index, column の名前、型が一致していないとダメ
df2.columns = ['d','e','f']
df1*df2
| a | b | c | d | e | f | |
|---|---|---|---|---|---|---|
| 1 | NaN | NaN | NaN | NaN | NaN | NaN |
| 2 | NaN | NaN | NaN | NaN | NaN | NaN |
| 3 | NaN | NaN | NaN | NaN | NaN | NaN |
df2.columns = ['a','b','c']
df2.index = [4,5,6]
df1*df2
| a | b | c | |
|---|---|---|---|
| 1 | NaN | NaN | NaN |
| 2 | NaN | NaN | NaN |
| 3 | NaN | NaN | NaN |
| 4 | NaN | NaN | NaN |
| 5 | NaN | NaN | NaN |
| 6 | NaN | NaN | NaN |
#indexを文字列として入力
df2.index = ['1','2','3']
df1*df2
| a | b | c | |
|---|---|---|---|
| 1 | NaN | NaN | NaN |
| 2 | NaN | NaN | NaN |
| 3 | NaN | NaN | NaN |
| 1 | NaN | NaN | NaN |
| 2 | NaN | NaN | NaN |
| 3 | NaN | NaN | NaN |
#indexを数値として入力
df2.index = [1,2,3]
df1*df2
| a | b | c | |
|---|---|---|---|
| 1 | 0.0 | 1.0 | 2.0 |
| 2 | 3.0 | 4.0 | 5.0 |
| 3 | 6.0 | 7.0 | 8.0 |
dataframe * 行、列
#df1の各行に、df2の2行目をかける
df1*df2.iloc[1]
| a | b | c | |
|---|---|---|---|
| 1 | 3.0 | 4.0 | 5.0 |
| 2 | 3.0 | 4.0 | 5.0 |
| 3 | 3.0 | 4.0 | 5.0 |
# df1の各列に、df2のc列をかける・・・よりスマートな方法がありそうですが。
(df1.T*df2['c']).T
| a | b | c | |
|---|---|---|---|
| 1 | 2.0 | 2.0 | 2.0 |
| 2 | 5.0 | 5.0 | 5.0 |
| 3 | 8.0 | 8.0 | 8.0 |
df1*df2.iloc[1,2]
| a | b | c | |
|---|---|---|---|
| 1 | 5.0 | 5.0 | 5.0 |
| 2 | 5.0 | 5.0 | 5.0 |
| 3 | 5.0 | 5.0 | 5.0 |
ようやくMCMC を概ね理解し動かせるようになったので、地下水モデル10章の不確実性に関する部分へ少しだけ立ち戻ることに。
https://phreeqc.blogspot.com/2019/07/4.html
MCMC関連の文献を探してみると、簡単な計算をされている事例がありました。第2著者以降は日本の方ですね。
Julien Boulange, Hirozumi Watanabe, Shinpei Akai (2017) A Markov Chain Monte Carlo technique for parameter estimation and inference in pesticide fate and transport modeling
尤度を計算する前に採択されたパラメータでシミュを回す必要があります。大きなモデルを扱う実務向きではないですね。
実務で扱う地下水モデルで不確実性を検討するには、PEST一択なのでしょうか？確かに、市販のパッケージソフトには組み込まれており、使い勝手も良かった印象はあります。加えて、逆解析を古くから扱っているせいか、Bayes も良く出てきます。海外では当たり前なのかな？
もう少し、文献等を読んでみましょう。
Ohki et al.(2020)
Landslide detection in mountainous forest areas using polarimetry and interferometric coherence
https://earth-planets-space.springeropen.com/articles/10.1186/s40623-020-01191-5
ALOS-2 のデータ （Lバンド、空間解像度約６m)を用いた、土砂移動箇所検出の検討事例です。北海道胆振東部地震（2018年9月）と九州北部豪雨（2017年7月）を対象としています。
提案手法は change detection の範疇ですが、画像データを利用していません。波のデータのみです。そして最後はお決りの機械学習です（といっても決定木です）。
特徴
結果
留意点
最近の国交省さんの動向を見ていますと、今後、土砂災害分野における各種衛星データの利用は拡大して行くようですね。関連マニュアルも複数出ていますし、ALOS-3へ向けてのシンポジウムでも以下のような報告がされています。
https://www.pco-prime.com/alos-3_alos-4_sympo/pdf/1505.pdf
機械学習を用いた change detection の精度は高いといえない状況ですが、いずれ向上し、災害時の人手不足問題も（change detection に関しては）回避できるでしょう。
検出アルゴリズムの構築以外は特に難しい作業はないので、備えておきましょう。
EXCELをよく使っていた頃、データに多項式をフィッティングする場面で何次式を使うかは感覚的でした。
低次だと合いませんし、高次だとノイズまで合わせてしまいます（いわゆる過学習）。ちょうどよいところを探す方法は、、、古くからありましたね。
LASSO
https://scikit-learn.org/stable/modules/linear_model.html#lasso
\begin{align*}
\min_{w}{ \frac{1}{2n_{\text{samples}}} ||y - X w||_2 ^ 2 + \alpha ||w||_1}\end{align*}
この損失関数はモデルの適合性を表す誤差項と、複雑さを表すペナルティー項で構成されています。αは求めたいパラメータwを制御する上位のパラメータ（ハイパーパラメータ）。αを大きくするほど大きなL1ノルムに対してペナルティーがかかる仕組みです。するとL1ノルムは小さくならざるを得なくなり、スパースな解が導かれる（が適合性は低くなりやすい）という仕組み。逆にαを小さくするほどL1ノルムを大きく取れるようになり、適合性が向上しやすい（が複雑で、過学習に陥りやすい）。αをどのように調整・採用するかがキモです。
その判断基準は以下の通り。
・情報量基準（AIC、BICなど）
・cross‐validation
スパースを仮定した場合、非ゼロのwを増やすと組み込まれるXも高次まで多く含まれて適合性は良くなる（が複雑になる）といった表現が理想です。が、L0ノルムを扱うと組み合わせ（計算コスト）が膨大になるとのこと。そこでα＊L1ノルムで代用しているそうです。
ちなみに、L1、L2 は機械学習でも多用されていましたね。回帰だと以下の通りです。
・L0正則化・・・非ゼロノルムの個数
・L1正則化・・・L1ノルム利用、LASSO回帰
・L2正則化・・・L2ノルム利用、Ridge回帰
・L1＋L2正則化・・・Elastic Net
scikit-learn を利用すれば、これらのモデリングは容易。ありがたいですね。
==========================================
ラグランジュの未定乗数法を適用した場合の図や説明も web 上には多くあります。そちらの方が直感的で理解し易いでしょう。
関連内容はコチラ↓
https://phreeqc.blogspot.com/2019/06/blog-post_19.html
https://phreeqc.blogspot.com/2020/02/blog-post.html
SAR データを利用した Change Detection を整理。
以前利用していたDocker は －vオプションを使用しなければ動きました。機能が追加されているようです。
https://phreeqc.blogspot.com/2019/04/sar-change-detection-on-google-earth.html
GEE コミュニティのチュートリアルにもありました。Part１は動きましたが、Part２の後半はダメでした。
https://developers.google.com/earth-engine/tutorials/community/detecting-changes-in-sentinel-1-imagery-pt-2
いずれも、GEE＋ Sentinel-1。解析エンジンもデータも海外製です。データは政府事業の Tellus 提供より揃っています。
https://phreeqc.blogspot.com/2019/02/google-earth-engine.html
国内の複数個所で Change Detection を試行しましたが、動作も問題ありません。きちんと結果は出てきます。
扱っているモノがモノなので、まだ精度面では物足りません。が、手軽に加工できるところはありがたいですね。将来的にはアルゴリズムが向上し、精度面でも耐え得るモノになるでしょう。
Tellus でも ALOS2 のデータが公開され始めました。が、未だ直近1年分くらいは公開されていません。過去分も、２～3年前くらいまで。直近のデータがないのは辛い。
しかも、開発環境を割り当ててもらえないと、自由に加工できない状況は継続中。https://phreeqc.blogspot.com/2020/09/tellusar.html
まだまだ海外製には追いつけない感が漂っています。頑張って欲しいですね。