【科研绘图】记录一次论文结果复现

复现原论文中的图片是科研的基本功之一，它不仅验证了研究结果的可靠性，确保了科学工作的准确性和可重复性，还深刻地评估了方法的有效性，体现了对原始研究的尊重和对科学过程的严谨态度。这个过程不仅提高了研究的透明度，促进了科学交流和进步，更在细致的操作中提升了个人和团队的技能。通过精确复现图形，研究者能够在继承和创新中不断推动科学领域的发展，为后续的研究奠定坚实的基础，体现了科研工作中的严谨与追求卓越。

任务：

给定海水温度数据（cmems_mod_glo_phy_my_0.083deg_P1M-m_thetao_70.00W-54.00W_70.00S-61.00S_0.49-902.34m_1993-01-01-2021-06-01.nc）和研究区域矢量（Ocean_roi_03.shp），请使用这些数据绘制如下示例图（黑线和灰色线不需要）。

原论文：

Wallis B J, Hogg A E, Meredith M P, et al. Ocean warming drives rapid dynamic activation of marine-terminating glacier on the west Antarctic Peninsula[J]. Nature Communications, 2023, 14(1): 7535.

绘制思路：

目标是从 NetCDF 文件中读取海洋温度数据，并结合 Shapefile 文件定义的区域掩膜，筛选出特定深度和地理区域内的数据，计算温度异常，并绘制结果图。具体步骤如下：

1. 导入库：

   • matplotlib.pyplot 用于绘图。
   • numpy 用于数组操作。
   • xarray 用于处理 NetCDF 数据。
   • geopandas 用于处理 Shapefile 文件。
   • salem 用于地理数据操作和掩膜。

2. 读取数据：

   • 使用 xarray 打开 NetCDF 文件，读取海洋温度数据 (thetao) 和深度数据 (depth)。
   • 使用 geopandas 读取 Shapefile 文件，获取感兴趣区域的边界框。

3. 筛选数据：

   • 根据 Shapefile 边界框的经纬度范围，筛选出在该区域内的温度数据。
   • 将深度数据限制在 0 到 500 米之间。

4. 创建掩膜：

   • 使用 salem 创建 Shapefile 区域的掩膜，并将其应用到筛选后的温度数据上。

5. 处理空白数据：

   • 对掩膜后的数据进行插值处理，填补时间和深度维度上的空白数据。

6. 补全深度数据：

   • 如果深度数据的最大值小于 500 米，则使用最近邻方法补全数据，使其覆盖 0 到 500 米的范围。

7. 计算温度平均值和异常值：

   • 计算每个月的平均温度。
   • 计算每个数据点的温度异常，即实际温度减去月平均温度。

8. 定义颜色映射表：

   • 自定义颜色映射表，将温度异常值映射到相应的颜色。

9. 设置颜色条：

   • 定义颜色条的分界线，使颜色条反映温度异常值的范围（-2°C 到 2°C）。

10. 创建图形：

    • 创建一个大小为 17x5 英寸的图形。
    • 提取月份和深度数据。
    • 计算每个深度层的平均温度异常值。

11. 绘制温度异常图：

    • 使用 pcolormesh 方法绘制温度异常值的网格图。
    • 添加颜色条，并设置颜色条的标签和刻度。

12. 设置坐标轴标签和刻度：

    • 设置纵轴标签为“Depth (m)”，字体为 Times New Roman，字号为 20。
    • 反转纵轴，使深度从 500 米到 0 米。
    • 设置纵轴刻度为每 100 米一档，字体为 Times New Roman，字号为 20。
    • 设置横轴为特定年份刻度，字体为 Times New Roman，字号为 20。

复现结果：

源代码：

import matplotlib.pyplot as plt
import numpy as np
import xarray as xr
import geopandas as gpd
import salem
from matplotlib.colors import LinearSegmentedColormap, BoundaryNorm
import warnings
warnings.filterwarnings("ignore")

# 打开 NetCDF 文件
file_path = 'D:/PyProject/Plot/cmems_mod_glo_phy_my_0.083deg_P1M-m_thetao_70.00W-54.00W_70.00S-61.00S_0.49-902.34m_1993-01-01-2021-06-01.nc'
data = xr.open_dataset(file_path)

# 打开 Shapefile
shapefile_path = 'D:/PyProject/Plot/Ocean_roi_03.shp'
shape = gpd.read_file(shapefile_path)

# 获取 Shapefile 的边界框
minx, miny, maxx, maxy = shape.total_bounds

# 筛选出在 Shapefile 边界框内的温度数据
temperature_data = data['thetao'].sel(longitude=slice(minx, maxx), latitude=slice(miny, maxy))

# 假设深度数据存在于 depth 变量中
depth = data['depth']
temperature_data = temperature_data.sel(depth=slice(0, 500))

# 创建一个 Shapefile 区域的掩模
roi = salem.read_shapefile(shapefile_path)
masked_temperature_data = temperature_data.salem.roi(shape=roi)

# 插值处理空白数据
masked_temperature_data = masked_temperature_data.interpolate_na(dim='time', method='linear')
masked_temperature_data = masked_temperature_data.interpolate_na(dim='depth', method='linear')

# 补全深度到500米
max_depth = masked_temperature_data['depth'].max().values
if max_depth < 500:
    additional_depths = np.arange(max_depth + 1, 501)
    additional_data = masked_temperature_data.isel(depth=-1).expand_dims('depth').reindex({'depth': additional_depths}, method='nearest')
    masked_temperature_data = xr.concat([masked_temperature_data, additional_data], dim='depth')

# 计算每月的平均温度
monthly_mean_temperature = masked_temperature_data.groupby('time.month').mean('time')

# 计算每个点的温度异常
temperature_anomaly = masked_temperature_data.groupby('time.month') - monthly_mean_temperature

# 自定义颜色映射表
colors = [
    (0.0, (0.055, 0.224, 0.422)),  # 深蓝色，代表-2度
    (0.25, '#6699cc'),
    (0.5, 'white'),
    (0.75, '#ff3300'),
    (1.0, (0.4, 0.122, 0.118))  # 深红色，代表2度
]
custom_cmap = LinearSegmentedColormap.from_list('custom_cmap', colors)

# 设置颜色条的分界线
boundaries = np.linspace(-2, 2, 17)  # 16个色块
norm = BoundaryNorm(boundaries, custom_cmap.N, clip=True)

# 创建图形
fig, ax = plt.subplots(figsize=(17, 5))

# 提取月份和深度数据
months = masked_temperature_data['time'].values
depths = masked_temperature_data['depth'].values

# 计算每个深度层的平均温度异常
temperature_anomaly_mean = temperature_anomaly.mean(dim=['longitude', 'latitude']).transpose()

# 创建pcolormesh图，设置颜色条的范围
c = ax.pcolormesh(months, depths, temperature_anomaly_mean, cmap=custom_cmap, norm=norm, shading='auto')

# 添加颜色条
cbar = fig.colorbar(c, ax=ax, boundaries=boundaries, extend='both')
cbar.set_label('Potential Temp Anomaly (C)', fontsize=20, fontname='Times New Roman')

# 设置颜色条刻度
cbar.set_ticks([-2, -1, 0, 1, 2])
cbar.ax.set_yticklabels(['-2', '-1', '0', '1', '2'], fontsize=20, fontname='Times New Roman')

# 设置轴标签和标题
# ax.set_xlabel('Year', fontsize=12, fontname='Times New Roman')
ax.set_ylabel('Depth (m)', fontsize=20, fontname='Times New Roman')
# ax.set_title('Monthly Temperature Anomaly', fontsize=14, fontname='Times New Roman')

# 反转纵轴并设置刻度
ax.set_ylim(500, 0)
ax.set_yticks([500, 400, 300, 200, 100, 0])
ax.set_yticklabels([500, 400, 300, 200, 100, 0], fontsize=20, fontname='Times New Roman')

# 设置特定的年份刻度和标签
years = [1995, 2000, 2005, 2010, 2015, 2020]
ax.set_xticks([np.datetime64(str(year)) for year in years])
ax.set_xticklabels(years, fontsize=20, fontname='Times New Roman')

# 设置字体
plt.xticks(fontsize=20, fontname='Times New Roman')

# 保存图形，提高分辨率
fig.savefig('temperature_anomaly.png', dpi=600)  # 设置图像分辨率为600 DPI

# 显示图形
plt.show()