点データの等高線マッピング作成時のマスク処理に関して

#軸がついてるやつ








#使い方
#CSV読み込みのところに適切なものを入れる
#df→データが入っているファイル名
#lat→そのデータを取得した緯度
#lon→そのデータを取得した経度
#moist→水分データ
#すべてダブルクォーテーション(")の間に先頭行の値を入れる（先頭行にデータを入れると処理されない）
#緯度経度のデータはRTKで取得したものを使用すると正確性が保たれる

#すべて入力したら実行する
#実行するとpngファイルが生成される。（スケール付き）
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
from matplotlib.ticker import ScalarFormatter, FormatStrFormatter
from scipy.interpolate import griddata
import alphashape
from shapely.geometry import Point
from shapely.ops import unary_union

from scipy.spatial import ConvexHull
from matplotlib.path import Path


title_size = 16
label_size = 14
tick_size = 12
cbar_size = 14
# ========================
# 1. CSV読み込み
# ========================
df = pd.read_csv("サンプルファイル.csv")

lat = df["緯度"].values
lon = df["経度"].values
moist = df["VWC%"].values

# ========================
# 2. ConvertHull で境界作成
# ========================
from scipy.spatial import ConvexHull

points = np.vstack((lon, lat)).T

hull = ConvexHull(points)
hull_points = points[hull.vertices]




lon_valid = lon
lat_valid = lat
moist_valid = moist

# ========================
# 3. グリッド作成（CSVの範囲そのまま）
# ========================
grid_res = 200
lon_min, lon_max = lon.min(), lon.max()
lat_min, lat_max = lat.min(), lat.max()

lat_grid = np.linspace(lat_min, lat_max, grid_res)
lon_grid = np.linspace(lon_min, lon_max, grid_res)
lon_mesh, lat_mesh = np.meshgrid(lon_grid, lat_grid)

# ========================
# 4. 補間（α-shape内の点のみ）
# ========================
moist_grid = griddata(
    (lon_valid, lat_valid),
    moist_valid,
    (lon_mesh, lat_mesh),
    method='cubic'
)

# ========================
# 5.ConvertHullマスク
# ========================
from matplotlib.path import Path

path = Path(hull_points)

mask = path.contains_points(
    np.vstack((lon_mesh.flatten(), lat_mesh.flatten())).T
).reshape(lon_mesh.shape)

moist_grid_masked = np.where(mask, moist_grid, np.nan)

# ========================
# 6. 描画（オフセット表示を無効化して生の緯度経度を表示）
# ========================
fig, ax = plt.subplots(figsize=(10,8))

contour = ax.contourf(
    lon_mesh, lat_mesh, moist_grid_masked,
    levels=15,
    cmap='viridis'
)
ax.contour(
    lon_mesh, lat_mesh, moist_grid_masked,
    levels=15,
    colors='black',
    linewidths=0.8   # 少し太く
)


cbar = plt.colorbar(contour, ax=ax)
cbar.set_label("VWC(%)", fontsize=cbar_size)
cbar.ax.tick_params(labelsize=tick_size)

ax.set_title("Moisture Contour Map", fontsize=title_size)
ax.set_xlabel("Longitude", fontsize=label_size)
ax.set_ylabel("Latitude", fontsize=label_size)

ax.tick_params(labelsize=tick_size)


# ▼ ここが重要：CSVの緯度経度を「そのまま」軸に使うための設定
ax.set_xlim(lon_min, lon_max)
ax.set_ylim(lat_min, lat_max)

# オフセット表示・指数表記を無効化
ax.ticklabel_format(style='plain', useOffset=False, axis='both')   # もっとも簡単
# あるいは、より厳密に：
ax.xaxis.set_major_formatter(ScalarFormatter(useOffset=False))
ax.yaxis.set_major_formatter(ScalarFormatter(useOffset=False))
# 桁数固定したい場合は次のように（例：小数点以下6桁）
ax.xaxis.set_major_formatter(FormatStrFormatter('%.3f'))
ax.yaxis.set_major_formatter(FormatStrFormatter('%.3f'))

# アスペクトを地理座標に近づけたい場合（任意）
# ax.set_aspect('equal', adjustable='box')

ax.scatter(lon, lat, c='black', s=20, edgecolors='white', label='Sampling Points',linewidths=0.5)
ax.legend(fontsize=12)

plt.tight_layout()
plt.savefig(
    "orp_map_slide.png",
    dpi=300,
    bbox_inches="tight"
)
plt.show()