任意の座標(i,j)で、隣接する座標(i+1,j)(i-1,j)間での流速差を求めたい

Question

前提・実現したいこと

二次元のxy座標に分布する流速VxとVyをCSVファイルから読み込み、最終的には流速から圧力を求める式につなげたいと考えております
まずは各座標(i,j)での流速データVx(i,j)、Vy(i,j)の中心差分を求めるため、任意の座標ごとのx方向に隣り合う座標の流速データVx(i+1,j)、Vx(i-1,j)を求めようとしましたが、うまく取得できません

解決した際は各xy方向での流速の差
Vx(i+1,j)-Vx(i-1,j)
Vy(i+1,j)-Vy(i-1,j)
Vx(i,j+1)-Vx(i,j-1)
Vy(i,j+1)-Vy(i,j-1)
を求める予定です

取得データ（計算格子）の間隔が画像座標で16pxであるため、実際にx方向で隣り合う座標は(i+16,j)と(i-16,j)、y方向では(i,j+16)と(i,j-16)となります
ご教授いただければ幸いです

発生している問題・エラーメッセージ

データの形は以下の通りです
csvファイルの流速データをpd.read_csvで読み込みました

import pandas as pd
df1 = pd.read_csv("./data/30deg_10x10plate_10f.csv", skiprows=2, usecols=[0, 1, 2, 3, 6, 7, 8],
                      names=['Pair', 'Num', 'X_px', 'Y_px', 'Vx_mm/s', 'Vy_mm/s', 'V_mm/s'])

下記が出力されます

          Pair    Num  X_px  Y_px  Vx_mm/s  Vy_mm/s  V_mm/s
0       Pair01      2    16     0      0.0      0.0     0.0
1       Pair01      3    32     0      0.0      0.0     0.0
2       Pair01      4    48     0      0.0      0.0     0.0
3       Pair01      5    64     0      0.0      0.0     0.0
4       Pair01      6    80     0      0.0      0.0     0.0
...        ...    ...   ...   ...      ...      ...     ...
110394  Pair10  11036  1840  1456      0.0      0.0     0.0
110395  Pair10  11037  1856  1456      0.0      0.0     0.0
110396  Pair10  11038  1872  1456      0.0      0.0     0.0
110397  Pair10  11039  1888  1456      0.0      0.0     0.0
110398  Pair10  11040  1904  1456      0.0      0.0     0.0

[110399 rows x 7 columns]

Pairはシチュエーションの数、Numは計算格子点の番号、X_px,Y_pxはピクセル単位での座標、Vx_mm/s Vy_mm/s V_mm/sはx方向の流速、y方向の流速、xy方向の流速(流速ベクトルの長さ)になります

該当のソースコード

# %% environment settings
import pandas as pd
df1 = pd.read_csv("./data/30deg_10x10plate_10f.csv", skiprows=2, usecols=[0, 1, 2, 3, 6, 7, 8],
                      names=['Pair', 'Num', 'X_px', 'Y_px', 'Vx_mm/s', 'Vy_mm/s', 'V_mm/s'])
print(df1)
# 各座標での(x±16,y)と(x,y±16)の流速データを並べればいい？
df1['X_px+16'] = df1['X_px'] + 16
df1['Y_px+16'] = df1['Y_px'] + 16
df1['X_px-16'] = df1['X_px'] - 16
df1['Y_px-16'] = df1['Y_px'] - 16
df1['X_px_2'] = df1['X_px']
df1['Y_px_2'] = df1['Y_px']
df1['X_px_2'] = df1['X_px']
df1['Y_px_2'] = df1['Y_px']
df1['Vx_xpx+16'] = df1['Vx_mm/s'].loc[(df1['X_px'] == df1['X_px+16']) & (df1['Y_px'] == df1['Y_px_2'])]
df1['Vx_xpx-16'] = df1['Vx_mm/s'].loc[(df1['X_px'] == df1['X_px-16']) & (df1['Y_px'] == df1['Y_px_2'])]
df1['Vx_ypx+16'] = df1['Vx_mm/s'].loc[(df1['Y_px'] == df1['Y_px+16']) & (df1['X_px'] == df1['X_px_2'])]
df1['Vx_ypx-16'] = df1['Vx_mm/s'].loc[(df1['Y_px'] == df1['Y_px-16']) & (df1['X_px'] == df1['X_px_2'])]
df1['Vy_xpx+16'] = df1['Vy_mm/s'].loc[(df1['X_px'] == df1['X_px+16']) & (df1['Y_px'] == df1['Y_px_2'])]
df1['Vy_xpx-16'] = df1['Vy_mm/s'].loc[(df1['X_px'] == df1['X_px-16']) & (df1['Y_px'] == df1['Y_px_2'])]
df1['Vy_ypx+16'] = df1['Vy_mm/s'].loc[(df1['Y_px'] == df1['Y_px+16']) & (df1['X_px'] == df1['X_px_2'])]
df1['Vy_ypx-16'] = df1['Vy_mm/s'].loc[(df1['Y_px'] == df1['Y_px-16']) & (df1['X_px'] == df1['X_px_2'])]

print(df1)
df1.to_csv(f"./output/vvv.csv")

試したこと

任意のxy座標に計算格子の間隔(±16px)と一致する座標の流速データの行を追加しようとしましたが、すべて空欄で出力されました
また複数の状況(Pair1-10)が混在しているせいかと考え

df1 = df1.groupby('Pair')

を追加してみましたが、

Traceback (most recent call last):
  File "C:/Users/ユーザー名/Pressure/Predict_pressure_from_velocity.py", line 17, in <module>
    df1['Vx_xpx+16'] = df1['Vx_mm/s'].loc[(df1['X_px'] == df1['X_px+16']) & (df1['Y_px'] == df1['Y_px_2'])]
  File "C:\Users\ユーザー名\anaconda3\envs\PIVcalc\lib\site-packages\pandas\core\groupby\groupby.py", line 752, in __getattr__
    raise AttributeError(
AttributeError: 'SeriesGroupBy' object has no attribute 'loc'

となりました

そもそもpandasでやろうとしていること自体が間違いで、numpyで処理すべきなのでしょうか？
まだ独学でPythonに触れ始めてから日が浅く、何か根本的な部分で致命的に間違えているように思えてなりません

ご教授いただければ幸いです
よろしくお願いいたします

補足情報（FW/ツールのバージョンなど）

Python3.8
Anaconda3
PyCharm 2021.1.2 (Community Edition)

Comments

[bsdfan]

Pairはどのように扱いたいのでしょうか？
Pair毎に分けて計算するのか、混ぜて計算するのか。

[okari]

ありがとうございます
Pair毎で分けて計算できればと考えております

Pair1が時系列で一つ目の流速分布
Pair2が二つ目・・・
と続いており、まずはそれぞれの流速分布から座標ごとでxy方向の差分をとりたく…

最初に1つのCSVファイルをからPairごとに分けて出力してから処理しようとも考えましたが、
差分を求め圧力への変換処理が終わった後、今度は共通する座標ごとに時系列間隔（Pair1-5）の中央値や平均値を求めるため、座標に割り振られた計算格子点の番号(Num)を用いて
df1 = df1.groupby('Num')
df1 = df1.rolling(window=5, min_periods=5).median()
等で処理する予定です

そのため、Pairが混在するファイルからでもPair毎に分けて各座標の差分を計算できればと思います

Answer 1

groupbyの使い方と、pandasでのデータの選択の仕方がちゃんと理解できていないと思われます。

.loc[(df1['X_px'] == df1['X_px+16']) & ...] では、同じ行のX_pxと、X_px+16が同じになっている行を選ぶという意味なので、そんな行は存在しないので空になります。
やりたいことは、X_pxに対して、X_px+16が同じになる行をとってきたいということだと思うのですが、それには merge が必要になります。
(pd.merge(df, df, left_on='X_px', right_on='X_px+16')みたいな)

そもそもpandasでやろうとしていること自体が間違いで、numpyで処理すべきなのでしょうか？

pandasだけでやろうとするとかなり面倒になりそうなので、私ならばnumpyも組み合わせて使います。
Pairでgroupbyした上で、
Vx等は、X_pxとY_pxを 1/16 したものをインデックスにもつ2次元のndarrayにして、処理しています。

import numpy as np

def func(df):
    # X_px, Y_px を数値インデックス化
    i = df['X_px'] // 16
    j = df['Y_px'] // 16

    # Vx を2次元配列に
    vx = np.full((i.max() + 1, j.max() + 1), np.nan)
    vx[i, j] = df['Vx_mm/s']

    # 差分
    vxx = np.full_like(vx, np.nan)
    vxy = np.full_like(vx, np.nan)
    vxx[1:-1] = vx[2:] - vx[:-2]
    vxy[:, 1:-1] = vx[:, 2:] - vx[:, :-2]

    # データフレームに戻す
    df['Vxx'] = vxx[i, j]
    df['Vxy'] = vxy[i, j]

    return df

df = df.groupby('Pair').apply(func)

print(df)

Comments

[okari]

ありがとうございます！！
いただいた関数を組み込んだところ、空欄にならず値が出力されました

なるほど、計算格子の間隔の分だけ割ったうえでnp.fullで配列を初期化し、
vx[i, j] = df['Vx_mm/s']で配列ごとの値を割り振ればよかったのですね…！

また、pandasを用いた場合の解決策も同時にご提示いただきありがとうございます
ご指摘の通り、すべてpandasでやろうとするときっと次の段階（流速の差分を圧力を求めるポアソン方程式に組みこんで値が収束するまで反復計算）で面倒になる予感がしておりました…

同様の関数をVyx、Vyyで作成し、二次元中心差分を各座標でとったあと、渦度を求める式に入れて求め、値が正しいか確認してみます
ご教授いただき誠にありがとうございました！