#!/usr/bin/python3
#!/usr/bin/env python3
# coding: utf-8
# 2025/05/29 made by Amano, and modefied by Tamagawa
# read LiDAR, hourly mean and plot
import subprocess
import os
import pandas as pd
import matplotlib.pyplot as plt
import pandas as pd
import numpy as np
import matplotlib.colors as mcolors
from matplotlib.ticker import FuncFormatter, LogLocator
directory_7z = r'/pub3/LASER/2022-2025屋上LIDAR/STA' # 解凍したいファイルが入っているディレクトリ
exdirectory = r'/pub3/LASER/2022-2025屋上LIDAR/STA/extract' # 解凍後の保存場所
os.makedirs(exdirectory, exist_ok=True)
# 1時間平均処理のために読み込むディレクトリのパス
# exdirectory
# 1時間平均値を置くディレクトリ
output_dir_1hr = r'/pub3/LASER/2022-2025屋上LIDAR/STA/1hr'
output_dir_Graph = output_dir_1hr + r'Graph'
# CSV file の保存ファイル
csvfilename = os.path.join(output_dir_1hr, 'hourly_data.csv')
# 保存先ディレクトリが存在しない場合は作成
os.makedirs(output_dir_1hr, exist_ok=True)
os.makedirs(output_dir_Graph, exist_ok=True)
def extract_7z_file(file_path, extract_to_folder):
command = [
'/bin/7z', # 7-Zipのインストールパス
'x', file_path, '-o'+exdirectory ,'-y'
]
subprocess.run(command, check=True)
print( "Before extract 7z" )
for filename in os.listdir(directory_7z):
if filename.endswith('.7z'):
print( filename )
file_path = os.path.join(directory_7z, filename)
filenamesta = os.path.splitext( file_path )[0] # sta ファイル名
#
extract_7z_file(file_path, exdirectory)
# if( os.path.isfile( filenamesta )) :
# # extract_7z_file(file_path, exdirectory)
# print(f'skip '+filename)
# else:
# extract_7z_file(file_path, exdirectory)
# print(f'解凍しました: '+filename+'|| ` '+filenamesta)
# # In[ ]:
#
# 対象ファイルのリストを作成（.sta ファイルを対象）
file_paths = [os.path.join(exdirectory, f) for f in os.listdir(exdirectory) if f.endswith('.sta')]
print("1")
# 空のリストにデータフレームを格納
dataframes = []
print("空のリストにデータフレームを格納")
# 各ファイルを読み込み、リストに追加
print( file_paths )
for file_path in file_paths:
try:
# print( file_path )
# データを読み込み、リストに追加
df = pd.read_csv(file_path, delimiter='\t', encoding="Shift-JIS", skiprows=41, parse_dates=True)
dataframes.append(df)
# print(f"{file_path} 読み込み成功")
except Exception as e:
print(f"{file_path} の読み込み中にエラーが発生しました: {e}")
# データフレームを1つに統合
if dataframes:
combined_df = pd.concat(dataframes, ignore_index=True)
# データフレームの範囲を確認
print("データフレームの先頭行:")
print(combined_df.head())
print("データ範囲:")
if 'Timestamp (end of interval)' in combined_df.columns:
print(combined_df['Timestamp (end of interval)'].min(), combined_df['Timestamp (end of interval)'].max())
else:
print("指定された列が存在しません")
else:
print("読み込まれたデータフレームがありません")
# combined_df
# 必要な列を確認
time_column = 'Timestamp (end of interval)' # 時間を表す列の名前
if time_column in combined_df.columns:
try:
#
print( combined_df.shape, combined_df.shape[0] )
# 時間列をdatetime型に変換（必要であれば）
combined_df[time_column] = pd.to_datetime(combined_df[time_column])
# 時間を1時間単位に丸める
combined_df['Hour'] = combined_df[time_column].dt.floor('h')
# 数値列を1時間ごとに集計（例: 平均）
hourly_data = combined_df.groupby('Hour').mean(numeric_only=True)
# 結果のデータフレームを表示
print("1時間ごとのデータ:")
print(hourly_data)
csvpath = os.path.join(output_dir_1hr, 'hourly_data.csv' )
# 必要なら結果をCSVに保存
hourly_data.to_csv(csvpath, index=True, encoding='utf-8')
except Exception as e:
print(f"1時間ごとのデータ作成中にエラーが発生しました: {e}")
else:
print(f"{time_column} 列がデータフレームに存在しません")
# グラフの列リスト
wind_speed_columns = [
'40m Wind Speed (m/s)', '60m Wind Speed (m/s)', '80m Wind Speed (m/s)',
'100m Wind Speed (m/s)', '120m Wind Speed (m/s)', '140m Wind Speed (m/s)',
'160m Wind Speed (m/s)', '180m Wind Speed (m/s)', '200m Wind Speed (m/s)',
'220m Wind Speed (m/s)', '240m Wind Speed (m/s)', '260m Wind Speed (m/s)',
'280m Wind Speed (m/s)', '300m Wind Speed (m/s)'
]
# y軸の固定値（高度）
heights = [40, 60, 80, 100, 120, 140, 160, 180, 200, 220, 240, 260, 280, 300]
# y軸の値に+30を加えたリスト
adjusted_heights = [h + 30 for h in heights]
# カスタムフォーマッタを作成
def integer_log_formatter(x, pos):
return f"{int(x):d}" if x >= 1 else ""
# 色相環の24色を生成（12:00が赤になるよう調整）
shift = 12 # 12:00を基準に色相をシフト
colors_hsv = [mcolors.hsv_to_rgb(((h + shift) % 24 / 24.0, 1.0, 1.0)) for h in range(24)]
# daily_groupsをhourly_data全体に基づいて作成
daily_groups = hourly_data.groupby(hourly_data.index.date) # 全ての日付を対象に処理
# グラフ描画と保存
for date, group in daily_groups:
plt.figure(figsize=(10, 6))
output_file = os.path.join(output_dir_Graph, f'Windspeed_{date.strftime("%Y%m%d")}.png')
print(output_file)
if((os.path.isfile( output_file )) and (os.path.getsize( output_file ) > 300000) ) :
print( 'exist ', os.path.getsize( output_file ) )
else:
print( "ELSE PLOT SAVE FIG"+output_file )
# 各時刻ごとのデータをプロット
for timestamp, row in group.iterrows():
x_values = row[wind_speed_columns].values # 風速データ x_values に
hour = timestamp.hour
# plt.savefig(output_file, dpi=300, bbox_inches='tight')
# 時刻に対応する色相環の色を取得
color = colors_hsv[hour]
plt.plot(x_values, adjusted_heights, marker='o', color=color, label=f'Time: {timestamp.strftime("%H:%M")}')
# グラフ設定
formatted_date = date.strftime('%Y/%m/%d')
plt.title(f'Windspeed (Gifu Univ) {formatted_date}')
plt.xlabel('Wind Speed (m/s)')
plt.ylabel('Height (m) (+30 (m))')
plt.yscale('log') # y軸を対数スケールに設定
# y軸の目盛りを調整（整数目盛りを使用）
plt.gca().yaxis.set_major_locator(LogLocator(base=10.0, subs=np.arange(1, 10) * 0.1, numticks=10))
plt.gca().yaxis.set_major_formatter(FuncFormatter(integer_log_formatter))
plt.grid(which='both', linestyle='--', linewidth=0.5) # 対数スケール用のグリッド設定
# 凡例の重複を防ぐため、ユニークなラベルだけ表示
handles, labels = plt.gca().get_legend_handles_labels()
unique_labels = dict(zip(labels, handles))
plt.legend(unique_labels.values(), unique_labels.keys(), title="Time of Day", bbox_to_anchor=(1.05, 1), loc='upper left')
plt.tight_layout()
plt.savefig(output_file, dpi=300, bbox_inches='tight')
# plt.show()
plt.close() # メモリを節約するためにグラフを閉じる
まだ良くわかってないが、動いてはいる。(2025.5.29 玉川) 学生用サンプルプログラム