# %%
# 分析月份干旱对物候的影响强度 - 修复版
def analyze_monthly_drought_impact():
    """分析不同月份发生的干旱事件对植被物候的影响强度"""
    
    # 设置字体
    plt.rcParams['font.family'] = 'SimSun'  # 宋体
    plt.rcParams['axes.unicode_minus'] = False  # 负号正常显示
    
    # 定义中文名称映射
    cn_names = {
        'Extreme_Atmospheric_Drought': '大气干旱',
        'Extreme_Soil_Drought': '土壤干旱',
        'Extreme_Hydrological_Drought': '水文干旱',
        'VPD_TWSA_Compound': 'VPD-TWSA复合干旱',
        'TWSA_SSMI_Compound': 'TWSA-SSMI复合干旱',
        'VPD_TWSA_SSMI_Compound': '三重复合干旱',
        'Valid_Dormancy_doy_mean': '休眠期',
        'Valid_Greenup_doy_mean': '返青期',
        'Valid_LOS_doy_mean': '生长季长度',
        'Valid_Peak_doy_mean': '峰值期'
    }
    
    # 1. 准备月度干旱数据
    monthly_drought = pd.DataFrame()
    
    for year in df['Year'].unique():
        for month in range(1, 13):
            month_data = df[(df['Year'] == year) & (df['month'] == month)]
            if not month_data.empty:
                row_dict = {'Year': year, 'month': month}
                for drought_type in case2:
                    if drought_type in month_data.columns:
                        row_dict[drought_type] = month_data[drought_type].iloc[0]
                monthly_drought = pd.concat([monthly_drought, pd.DataFrame([row_dict])])
    
    # 2. 准备物候数据
    pheno_data_copy = pheno_data.copy()
    pheno_data_copy['Year'] = pheno_data_copy['Year'].astype(int)
    
    # 3. 分析每个月份干旱对物候的影响
    monthly_impact = pd.DataFrame()
    
    # 考虑前一年的月份对当年物候的影响
    for lag in [0, 1]:  # 0表示当年，1表示前一年
        for month in range(1, 13):
            for drought_type in case2:
                for phase in phases_doy:
                    # 创建月度干旱与物候的配对数据
                    paired_data = pd.DataFrame()
                    
                    for year in pheno_data_copy['Year'].unique():
                        target_year = year - lag  # 对应的干旱年份
                        
                        # 获取该年该月的干旱状态
                        drought_status = monthly_drought[
                            (monthly_drought['Year'] == target_year) & 
                            (monthly_drought['month'] == month)
                        ]
                        
                        # 获取物候数据
                        pheno_value = pheno_data_copy[pheno_data_copy['Year'] == year]
                        
                        if not drought_status.empty and not pheno_value.empty:
                            paired_data = pd.concat([paired_data, pd.DataFrame({
                                'drought': [drought_status[drought_type].iloc[0]],
                                'phenology': [pheno_value[phase].iloc[0]]
                            })])
                    
                    # 计算相关性
                    if len(paired_data) > 3:  # 确保有足够的数据点
                        from scipy import stats
                        corr, p_value = stats.pearsonr(paired_data['drought'], paired_data['phenology'])
                        
                        monthly_impact = pd.concat([monthly_impact, pd.DataFrame({
                            'lag': [lag],
                            'month': [month],
                            'drought_type': [drought_type],
                            'phenology_phase': [phase],
                            'correlation': [corr],
                            'p_value': [p_value]
                        })])
    
    # 4. 可视化分析结果
    
    # 设置月份标签
    month_labels = ['一月', '二月', '三月', '四月', '五月', '六月', 
                   '七月', '八月', '九月', '十月', '十一月', '十二月']
    
    # 4.1 热图展示月份干旱对物候的影响
    for phase in phases_doy:
        # 创建一个新的图形，包含多个子图
        fig, axes = plt.subplots(2, 3, figsize=(15, 10), dpi=300)
        axes = axes.flatten()
        
        for i, drought_type in enumerate(case2):
            ax = axes[i]
            
            # 筛选数据
            drought_data = monthly_impact[
                (monthly_impact['phenology_phase'] == phase) & 
                (monthly_impact['drought_type'] == drought_type)
            ]
            
            # 创建一个完整的月份-滞后年数网格
            full_grid = pd.DataFrame(
                index=[0, 1],  # 滞后年数
                columns=range(1, 13)  # 月份
            )
            
            # 填充相关系数
            for _, row in drought_data.iterrows():
                full_grid.loc[row['lag'], row['month']] = row['correlation']
            
            # 绘制热图
            im = sns.heatmap(full_grid, ax=ax, cmap='RdBu_r', vmin=-1, vmax=1, 
                           annot=True, fmt='.2f', cbar=True)
            
            # 设置标题和标签
            ax.set_title(f'{cn_names.get(drought_type, drought_type)}', fontsize=12, fontweight='bold')
            ax.set_xlabel('月份')
            ax.set_ylabel('滞后年数')
            
            # 设置x轴标签 - 修复刻度问题
            ax.set_xticks(np.arange(0.5, 12.5))
            ax.set_xticklabels(month_labels, rotation=45)
            
            # 设置y轴标签
            ax.set_yticks([0.5, 1.5])
            ax.set_yticklabels(['当年', '前一年'])
            
            # 标记显著性
            for _, row in drought_data.iterrows():
                if row['p_value'] < 0.05:
                    ax.text(row['month'] - 0.5, row['lag'] + 0.5, '*', 
                           ha='center', va='center', color='black', fontsize=15)
        
        # 调整布局
        plt.tight_layout()
        plt.subplots_adjust(top=0.9)
        plt.suptitle(f'{cn_names.get(phase, phase)}对月度干旱的敏感性', fontsize=16, fontweight='bold')
        
        # 保存图像
        plt.savefig(f'./干旱对植被的影响/{cn_names.get(phase, phase)}月度干旱敏感性.jpeg', dpi=800, bbox_inches='tight')
        plt.show()
    
    # 4.2 找出影响最大的月份
    # 按物候期分组，找出每种干旱类型影响最大的月份
    max_impact = pd.DataFrame()
    
    for phase in phases_doy:
        phase_data = monthly_impact[monthly_impact['phenology_phase'] == phase]
        
        for drought_type in case2:
            drought_data = phase_data[phase_data['drought_type'] == drought_type]
            
            # 找出相关性绝对值最大的月份
            if not drought_data.empty:
                max_idx = drought_data['correlation'].abs().idxmax()
                max_row = drought_data.loc[max_idx]
                
                max_impact = pd.concat([max_impact, pd.DataFrame({
                    'phenology_phase': [phase],
                    'drought_type': [drought_type],
                    'max_impact_month': [max_row['month']],
                    'lag': [max_row['lag']],
                    'correlation': [max_row['correlation']],
                    'p_value': [max_row['p_value']]
                })])
    
    # 4.3 可视化最大影响月份
    plt.figure(figsize=(12, 8), dpi=300)
    
    # 创建表格数据
    table_data = []
    for phase in phases_doy:
        row = [cn_names.get(phase, phase)]
        for drought_type in case2:
            subset = max_impact[
                (max_impact['phenology_phase'] == phase) & 
                (max_impact['drought_type'] == drought_type)
            ]
            if not subset.empty:
                month = int(subset['max_impact_month'].iloc[0])
                lag = int(subset['lag'].iloc[0])
                corr = subset['correlation'].iloc[0]
                p_val = subset['p_value'].iloc[0]
                
                month_label = month_labels[month-1]
                lag_label = "当年" if lag == 0 else "前一年"
                sig_mark = "*" if p_val < 0.05 else ""
                
                cell_text = f"{lag_label}{month_label}\n(r={corr:.2f}){sig_mark}"
                row.append(cell_text)
            else:
                row.append("N/A")
        table_data.append(row)
    
    # 创建表格
    col_labels = ['物候期'] + [cn_names.get(d, d) for d in case2]
    
    ax = plt.subplot(111, frame_on=False)
    ax.xaxis.set_visible(False)
    ax.yaxis.set_visible(False)
    
    table = plt.table(cellText=table_data, colLabels=col_labels, loc='center', cellLoc='center')
    table.auto_set_font_size(False)
    table.set_fontsize(12)
    table.scale(1, 1.5)
    
    # 设置单元格颜色
    for i in range(len(table_data)):
        for j in range(1, len(col_labels)):
            cell = table[(i+1, j)]
            if "N/A" not in cell.get_text().get_text():
                # 提取相关系数
                import re
                corr_match = re.search(r'r=([-\d.]+)', cell.get_text().get_text())
                if corr_match:
                    corr = float(corr_match.group(1))
                    if corr > 0.5:
                        cell.set_facecolor('#ffcccc')  # 浅红色
                    elif corr < -0.5:
                        cell.set_facecolor('#ccccff')  # 浅蓝色
                    
                    # 添加显著性标记的颜色
                    if "*" in cell.get_text().get_text():
                        cell.set_text_props(weight='bold')
    
    plt.title('./干旱对植被的影响/各物候期对月度干旱的最大敏感性', fontsize=16, fontweight='bold')
    plt.tight_layout()
    plt.savefig('./干旱对植被的影响/物候期月度干旱最大敏感性.jpeg', dpi=800, bbox_inches='tight')
    plt.show()
    
    # 5. 创建月度干旱影响强度柱状图
    plt.figure(figsize=(15, 10), dpi=300)
    
    for i, phase in enumerate(phases_doy):
        plt.subplot(2, 2, i+1)
        
        # 筛选数据 - 只看当年的影响
        phase_data = monthly_impact[
            (monthly_impact['phenology_phase'] == phase) & 
            (monthly_impact['lag'] == 0)
        ]
        
        # 按干旱类型分组
        grouped_data = []
        for drought_type in case2:
            drought_data = phase_data[phase_data['drought_type'] == drought_type]
            if not drought_data.empty:
                # 按月份排序
                drought_data = drought_data.sort_values('month')
                grouped_data.append({
                    'drought_type': drought_type,
                    'months': drought_data['month'].values,
                    'correlations': drought_data['correlation'].values,
                    'p_values': drought_data['p_value'].values
                })
        
        # 设置柱状图位置
        bar_width = 0.15
        r = np.arange(12)
        
        # 绘制柱状图
        for j, data in enumerate(grouped_data):
            bars = plt.bar(r + j*bar_width, data['correlations'], width=bar_width, 
                          label=cn_names.get(data['drought_type'], data['drought_type']))
            
            # 标记显著性
            for k, p_val in enumerate(data['p_values']):
                if p_val < 0.05:
                    plt.text(r[k] + j*bar_width, data['correlations'][k] + 0.02, '*', 
                            ha='center', va='bottom', fontsize=12)
        
        # 设置图表属性
        plt.title(f'{cn_names.get(phase, phase)}', fontsize=14, fontweight='bold')
        plt.xlabel('月份')
        plt.ylabel('相关系数')
        plt.xticks(r + bar_width*2.5, month_labels)
        plt.axhline(y=0, color='black', linestyle='-', alpha=0.3)
        plt.legend(loc='best')
        plt.grid(axis='y', linestyle='--', alpha=0.3)
    
    plt.tight_layout()
    plt.savefig('./干旱对植被的影响/月度干旱影响强度柱状图.jpeg', dpi=800, bbox_inches='tight')
    plt.show()
    
    return monthly_impact, max_impact

# 执行分析
monthly_impact, max_impact = analyze_monthly_drought_impact()  {
	"name": "TypeError",
	"message": "Image data of dtype object cannot be converted to float",
	"stack": "---------------------------------------------------------------------------
TypeError                                 Traceback (most recent call last)
Cell In[13], line 292
    289     return monthly_impact, max_impact
    291 # 执行分析
--> 292 monthly_impact, max_impact = analyze_monthly_drought_impact()

Cell In[13], line 116, in analyze_monthly_drought_impact()
    113     full_grid.loc[row['lag'], row['month']] = row['correlation']
    115 # 绘制热图
--> 116 im = sns.heatmap(full_grid, ax=ax, cmap='RdBu_r', vmin=-1, vmax=1, 
    117                annot=True, fmt='.2f', cbar=True)
    119 # 设置标题和标签
    120 ax.set_title(f'{cn_names.get(drought_type, drought_type)}', fontsize=12, fontweight='bold')

File d:\\Python310\\lib\\site-packages\\seaborn\\matrix.py:459, in heatmap(data, vmin, vmax, cmap, center, robust, annot, fmt, annot_kws, linewidths, linecolor, cbar, cbar_kws, cbar_ax, square, xticklabels, yticklabels, mask, ax, **kwargs)
    457 if square:
    458     ax.set_aspect(\"equal\")
--> 459 plotter.plot(ax, cbar_ax, kwargs)
    460 return ax

File d:\\Python310\\lib\\site-packages\\seaborn\\matrix.py:306, in _HeatMapper.plot(self, ax, cax, kws)
    303     kws.setdefault(\"vmax\", self.vmax)
    305 # Draw the heatmap
--> 306 mesh = ax.pcolormesh(self.plot_data, cmap=self.cmap, **kws)
    308 # Set the axis limits
    309 ax.set(xlim=(0, self.data.shape[1]), ylim=(0, self.data.shape[0]))

File d:\\Python310\\lib\\site-packages\\matplotlib\\__init__.py:1473, in _preprocess_data.<locals>.inner(ax, data, *args, **kwargs)
   1470 @functools.wraps(func)
   1471 def inner(ax, *args, data=None, **kwargs):
   1472     if data is None:
-> 1473         return func(
   1474             ax,
   1475             *map(sanitize_sequence, args),
   1476             **{k: sanitize_sequence(v) for k, v in kwargs.items()})
   1478     bound = new_sig.bind(ax, *args, **kwargs)
   1479     auto_label = (bound.arguments.get(label_namer)
   1480                   or bound.kwargs.get(label_namer))

File d:\\Python310\\lib\\site-packages\\matplotlib\\axes\\_axes.py:6434, in Axes.pcolormesh(self, alpha, norm, cmap, vmin, vmax, shading, antialiased, *args, **kwargs)
   6430 coords = np.stack([X, Y], axis=-1)
   6432 kwargs.setdefault('snap', mpl.rcParams['pcolormesh.snap'])
-> 6434 collection = mcoll.QuadMesh(
   6435     coords, antialiased=antialiased, shading=shading,
   6436     array=C, cmap=cmap, norm=norm, alpha=alpha, **kwargs)
   6437 collection._scale_norm(norm, vmin, vmax)
   6439 coords = coords.reshape(-1, 2)  # flatten the grid structure; keep x, y

File d:\\Python310\\lib\\site-packages\\matplotlib\\collections.py:2146, in QuadMesh.__init__(self, coordinates, antialiased, shading, **kwargs)
   2144 kwargs.setdefault(\"pickradius\", 0)
   2145 super().__init__(coordinates=coordinates, shading=shading)
-> 2146 Collection.__init__(self, **kwargs)
   2148 self._antialiased = antialiased
   2149 self._bbox = transforms.Bbox.unit()

File d:\\Python310\\lib\\site-packages\\matplotlib\\collections.py:206, in Collection.__init__(self, edgecolors, facecolors, linewidths, linestyles, capstyle, joinstyle, antialiaseds, offsets, offset_transform, norm, cmap, pickradius, hatch, urls, zorder, **kwargs)
    203 self._offset_transform = offset_transform
    205 self._path_effects = None
--> 206 self._internal_update(kwargs)
    207 self._paths = None

File d:\\Python310\\lib\\site-packages\\matplotlib\\artist.py:1216, in Artist._internal_update(self, kwargs)
   1209 def _internal_update(self, kwargs):
   1210     \"\"\"
   1211     Update artist properties without prenormalizing them, but generating
   1212     errors as if calling `set`.
   1213 
   1214     The lack of prenormalization is to maintain backcompatibility.
   1215     \"\"\"
-> 1216     return self._update_props(
   1217         kwargs, \"{cls.__name__}.set() got an unexpected keyword argument \"
   1218         \"{prop_name!r}\")

File d:\\Python310\\lib\\site-packages\\matplotlib\\artist.py:1192, in Artist._update_props(self, props, errfmt)
   1189             if not callable(func):
   1190                 raise AttributeError(
   1191                     errfmt.format(cls=type(self), prop_name=k))
-> 1192             ret.append(func(v))
   1193 if ret:
   1194     self.pchanged()

File d:\\Python310\\lib\\site-packages\\matplotlib\\collections.py:2020, in _MeshData.set_array(self, A)
   2015     if shape not in ok_shapes:
   2016         raise ValueError(
   2017             f\"For X ({width}) and Y ({height}) with {self._shading} \"
   2018             f\"shading, A should have shape \"
   2019             f\"{' or '.join(map(str, ok_shapes))}, not {A.shape}\")
-> 2020 return super().set_array(A)

File d:\\Python310\\lib\\site-packages\\matplotlib\\cm.py:452, in ScalarMappable.set_array(self, A)
    450 A = cbook.safe_masked_invalid(A, copy=True)
    451 if not np.can_cast(A.dtype, float, \"same_kind\"):
--> 452     raise TypeError(f\"Image data of dtype {A.dtype} cannot be \"
    453                     \"converted to float\")
    455 self._A = A
    456 if not self.norm.scaled():

TypeError: Image data of dtype object cannot be converted to float"
}