Mekong River Basin Precipitation and Elevation Analysis

import xarray as xr
import geopandas as gpd
import matplotlib.pyplot as plt
from matplotlib.ticker import FormatStrFormatter

# 加载bcc数据
dataset = xr.open_dataset(r'I:\CMIP6HR\highresSST-present\bcc.nc')
pr_mean = dataset['pr_mean']
pr_frequency = dataset['pr_frequency']
pr_intensity = dataset['pr_intensity']

# 加载ERA5数据
era5_dataset = xr.open_dataset(r'I:\ERA5\era5.nc')
tp_average = era5_dataset['tp_average']
tp_frequency = era5_dataset['tp_frequency']
tp_intensity = era5_dataset['tp_intensity']

# 加载mekongDEM数据
dem_dataset = xr.open_dataset(r'I:\mekongDEM\demnc.nc')
mekong_clip = dem_dataset['mekong_clip']

# 加载shp文件
shapefile = gpd.read_file(r'I:\MekongRiverBasin_whole\MekongRiverBasin.shp')
lat_min, lat_max = shapefile.bounds['miny'].min(), shapefile.bounds['maxy'].max()
lon_min, lon_max = shapefile.bounds['minx'].min(), shapefile.bounds['maxx'].max()

# 筛选出shpfile内部的经纬度范围
pr_mean = pr_mean.sel(lat=slice(lat_min, lat_max), lon=slice(lon_min, lon_max))
pr_frequency = pr_frequency.sel(lat=slice(lat_min, lat_max), lon=slice(lon_min, lon_max))
pr_intensity = pr_intensity.sel(lat=slice(lat_min, lat_max), lon=slice(lon_min, lon_max))
mekong_clip = mekong_clip.sel(lat=slice(lat_min, lat_max), lon=slice(lon_min, lon_max))

# 计算每一行格点的平均值
pr_mean_avg = pr_mean.mean(dim='lon')
pr_frequency_avg = pr_frequency.mean(dim='lon')
pr_intensity_avg = pr_intensity.mean(dim='lon')
mekong_clip_avg = mekong_clip.mean(dim='lon')

# 创建Figure和Axes对象
fig, axes = plt.subplots(1, 3, figsize=(12, 5), sharey=False)

# 绘制pr_mean变量的纬度变化曲线
axes[0].plot(pr_mean_avg.coords['lat'], pr_mean_avg)
axes[0].set_xlabel('', fontsize=12)
axes[0].set_ylabel('precipitation average/(mm/d)', fontsize=12)
axes[0].set_title('', fontsize=12)
axes[0].xaxis.set_major_formatter(FormatStrFormatter('%d°N'))
axes[0].set_ylim([0, 20])  # 设置y轴刻度范围

# 绘制pr_frequency变量的纬度变化曲线
axes[1].plot(pr_frequency_avg.coords['lat'], pr_frequency_avg)
axes[1].set_xlabel('', fontsize=12)
axes[1].set_ylabel('precipitation frequency/%', fontsize=12)
axes[1].set_title('', fontsize=12)
axes[1].xaxis.set_major_formatter(FormatStrFormatter('%d°N'))
axes[1].set_ylim([0, 40])  # 设置y轴刻度范围

# 绘制pr_intensity变量的纬度变化曲线
axes[2].plot(pr_intensity_avg.coords['lat'], pr_intensity_avg)
axes[2].set_xlabel('', fontsize=12)
axes[2].set_ylabel('precipitation intensity/(mm/h)', fontsize=12)
axes[2].set_title('', fontsize=12)
axes[2].xaxis.set_major_formatter(FormatStrFormatter('%d°N'))
axes[2].set_ylim([0, 4.0])  # 设置y轴刻度范围

# 绘制mekong_clip变量的纬度变化曲线
axes2 = axes[2].twinx()
axes2.plot(mekong_clip_avg.coords['lat'], mekong_clip_avg, color='black')
axes2.set_ylabel('Elevation/m', fontsize=12)
axes2.set_ylim([0, 500])  # 设置y轴刻度范围
axes2.spines['right'].set_position(('outward', 50))  # 调整坐标轴位置

plt.tight_layout()
plt.show()

This code loads precipitation data from BCC and ERA5 datasets, as well as elevation data from a DEM file for the Mekong River Basin. It then calculates the average precipitation and frequency along with elevation profiles. The results are visualized in three subplots, with the elevation profile displayed on a secondary y-axis in the third subplot. The code uses the xarray and geopandas libraries for data manipulation and matplotlib for visualization. You can adjust the y-axis limits and the position of the secondary y-axis as needed.