TimeseriesPlots_AverageForChinaRegion_All.py

from paths import data_path
import pandas as pd
import matplotlib.pyplot as plt
from tqdm import tqdm
import os
import warnings
warnings.filterwarnings('ignore')

# Get the absolute path of the current script
script_path = os.path.abspath(__file__)
script_dir = os.path.dirname(script_path)
os.chdir(script_dir)

# ====================== 全局样式配置（核心美化） ======================
# 设置字体（优先Times New Roman，兼容中文显示）
plt.rcParams["font.family"] = ["Times New Roman", "SimHei", "serif"]
plt.rcParams['axes.unicode_minus'] = False  # 解决负号显示问题
plt.rcParams['figure.dpi'] = 100  # 默认分辨率
plt.rcParams['savefig.dpi'] = 300  # 保存图片分辨率
plt.rcParams['axes.linewidth'] = 1.2  # 坐标轴边框宽度
plt.rcParams['grid.alpha'] = 0.3  # 网格透明度
plt.rcParams['grid.linewidth'] = 0.8  # 网格线宽度

# --- 配置项 ---
REGION_MASK_FILE = data_path("DataFusion_China_CleanAir", "ChinaGrid_Region_Mask.csv")
DATA_DIR = data_path("DataFusion_China_CleanAir", "Data", "AtFOutput_WithoutCV_All")
OUTPUT_DIR = data_path("DatGATaFusion_China_CleanAir", "Picture", "ChinaRegion_Timeseries_All")

# 要绘制的指标（需和CSV列名完全匹配）
metrics = ['model', 'vna_pm25', 'evna_pm25', 'avna_pm25', 'CHAP', 'TAP']

# 美化配置
method_display_names = {
    'model': 'Model',
    'vna_pm25': 'VNA',
    'evna_pm25': 'eVNA',
    'avna_pm25': 'aVNA',
    'CHAP': 'CHAP',
    'TAP': 'TAP',
}

# 专业配色（符合学术图表规范）
method_colors = {
    'model': '#66c2a5',    # Teal (青绿)
    'vna_pm25': '#fc8d62', # Orange (橙色)
    'evna_pm25': '#8da0cb',# Lilac (淡紫)
    'avna_pm25': '#e78ac3',# Pink (粉色)
    'CHAP': '#377eb8',     # Blue (蓝色)
    'TAP': '#4daf4a',      # Green (绿色)
}

# 不同标记点
method_markers = {
    'model': '*',
    'vna_pm25': 'o',
    'evna_pm25': 's',
    'avna_pm25': '^',
    'CHAP': '>',
    'TAP': 'D',
}

# 地区显示名
region_display_names = {
    'BTHS': 'Beijing-Tianjin-Hebei-Shandong (BTHS)',
    'YRD': 'Yangtze River Delta (YRD)',
    'PRD': 'Pearl River Delta (PRD)',
    'Others': 'Other Regions',
    'China': 'Mainland China (All)'
}

# ====================== 数据加载 ======================
print(f"Reading region mask from {REGION_MASK_FILE}...")
try:
    region_map = pd.read_csv(REGION_MASK_FILE)
    if 'Region' not in region_map.columns:
        raise ValueError("Region mask file must contain 'Region' column")
    region_map = region_map[['ROW', 'COL', 'Region']].dropna()
    region_map['ROW'] = region_map['ROW'].astype(int)
    region_map['COL'] = region_map['COL'].astype(int)
    unique_regions = sorted(region_map['Region'].unique())
    print(f"Found regions: {unique_regions}")
except Exception as e:
    print(f"Error loading region mask: {e}")
    exit()

# 加载年度数据
years = list(range(2014, 2021))
all_data = pd.DataFrame()

print("Reading data for each year...")
for year in tqdm(years):
    file_name = f'{year}_AtF(Annual)_All_OnlyChina.csv'
    file_path = os.path.join(DATA_DIR, file_name)
    
    if os.path.exists(file_path):
        try:
            data = pd.read_csv(file_path)
            if 'Year' not in data.columns:
                data['Year'] = year
            
            # 合并地区信息
            data = data.merge(region_map, on=['ROW', 'COL'], how='left')
            all_data = pd.concat([all_data, data], ignore_index=True)
            
        except Exception as e:
            print(f"Error reading {file_path}: {e}")
    else:
        print(f"File {file_path} not found.")

if all_data.empty:
    print("No data loaded. Exiting.")
    exit()

# 数据校验：检查列是否存在
missing_cols = [m for m in metrics if m not in all_data.columns]
if missing_cols:
    print(f"Warning: Missing columns in data: {missing_cols}")
    metrics = [m for m in metrics if m not in missing_cols]
    if not metrics:
        print("No valid metrics to plot. Exiting.")
        exit()

# ====================== 绘图逻辑（核心修复+美化） ======================
os.makedirs(OUTPUT_DIR, exist_ok=True)
print(f"Generating time series plots to {OUTPUT_DIR}...")

regions_to_plot = unique_regions + ['China']
periods = all_data['Period'].unique() if 'Period' in all_data.columns else ['Annual']

for period in tqdm(periods, desc="Processing periods"):
    # 按周期筛选数据
    if 'Period' in all_data.columns:
        period_data = all_data[all_data['Period'] == period].copy()
    else:
        period_data = all_data.copy()

    if period_data.empty:
        continue

    # ====================== 新增：计算统一的Y轴范围 ======================
    print(f"Calculating unified Y-axis limits for period: {period}...")
    period_min_val = float('inf')
    period_max_val = float('-inf')
    
    # 预遍历所有地区和指标，找出最大最小值
    for r_name in regions_to_plot:
        if r_name == 'China':
            r_data = period_data
        else:
            r_data = period_data[period_data['Region'] == r_name]
            
        if r_data.empty: continue
        
        for method in metrics:
            if method not in r_data.columns: continue
            
            # 计算该地区该指标的年度均值（与绘图逻辑一致）
            means = r_data.groupby('Year')[method].mean()
            if not means.empty:
                period_min_val = min(period_min_val, means.min())
                period_max_val = max(period_max_val, means.max())
    
    # 计算统一的Y轴范围（带缓冲）
    if period_min_val != float('inf') and period_max_val != float('-inf'):
        padding = (period_max_val - period_min_val) * 0.15
        padding = max(padding, 2.0)
        unified_ylim_min = max(0, period_min_val - padding)
        unified_ylim_max = period_max_val + padding
    else:
        unified_ylim_min, unified_ylim_max = 0, 100 # 默认值
        
    print(f"Unified Y-axis range: [{unified_ylim_min:.2f}, {unified_ylim_max:.2f}]")

    for region_name in tqdm(regions_to_plot, desc="Processing regions", leave=False):
        # 筛选地区数据
        if region_name == 'China':
            region_period_data = period_data.copy()
        else:
            region_period_data = period_data[period_data['Region'] == region_name].copy()
        
        if region_period_data.empty:
            print(f"Warning: No data for region {region_name}, skipping.")
            continue
        
        # 创建图表（调整尺寸为更专业的比例）
        fig, ax = plt.subplots(figsize=(10, 6))
        
        min_val = float('inf')
        max_val = float('-inf')
        
        # ========== 修复缩进：这个循环必须在region内部 ==========
        # 绘制每个指标
        for method in metrics:
            if method not in region_period_data.columns:
                continue
            
            # 按年份计算均值
            yearly_data = region_period_data.groupby('Year')[method].mean().reset_index()
            
            if yearly_data.empty:
                print(f"Warning: No data for method {method} in region {region_name}")
                continue
            
            # 获取美化配置
            display_name = method_display_names.get(method, method)
            marker = method_markers.get(method, 'o')
            color = method_colors.get(method, '#888888')
            
            # 绘制折线图
            ax.plot(yearly_data['Year'].values, yearly_data[method].values,
                    marker=marker, markersize=9, markeredgecolor='white', markeredgewidth=1,
                    linestyle='-', linewidth=2.5, color=color, label=display_name,
                    alpha=0.9)
            
            # 更新y轴范围
            min_val = min(min_val, yearly_data[method].min())
            max_val = max(max_val, yearly_data[method].max())
        
        # ========== 图表美化 ==========
        # Y轴范围调整（使用统一范围）
        ax.set_ylim(unified_ylim_min, unified_ylim_max)
        
        # 设置坐标轴
        display_region = region_display_names.get(region_name, region_name)
        ax.set_title(f'{period}: {display_region} PM2.5 Time Series (2014-2020)', 
                     fontsize=16, fontweight='bold', pad=20)
        ax.set_xlabel('Year', fontsize=14, labelpad=10)
        ax.set_ylabel('PM2.5 (μg/m³)', fontsize=14, labelpad=10)
        
        # X轴设置
        ax.set_xticks(years)
        ax.set_xticklabels(years, fontsize=12)
        ax.tick_params(axis='x', rotation=0)  # 年份不需要旋转
        
        # Y轴设置
        ax.tick_params(axis='y', labelsize=12)
        
        # 网格线（只显示水平网格，更简洁）
        ax.grid(True, axis='y', linestyle='--', alpha=0.7)
        ax.grid(False, axis='x')
        
        # 边框优化（只保留上/右边框的刻度，隐藏边框线）
        ax.spines['top'].set_visible(False)
        ax.spines['right'].set_visible(False)
        ax.spines['left'].set_color('#CCCCCC')
        ax.spines['bottom'].set_color('#CCCCCC')
        
        # 图例（放置在右上角，带背景美化）
        legend = ax.legend(loc='upper right', fontsize=12, frameon=True, 
                           bbox_to_anchor=(0.98, 0.98), borderaxespad=0,
                           framealpha=1, shadow=True)
        legend.get_frame().set_facecolor('white')
        legend.get_frame().set_edgecolor('#EEEEEE')
        
        # 调整布局
        plt.tight_layout()
        
        # 保存图片
        filename = f"{OUTPUT_DIR}/{period}_{region_name}_PM25_Timeseries.png"
        plt.savefig(filename, dpi=300, bbox_inches='tight', facecolor='white')
        plt.close()

print("All time series plots generated successfully!")