Mekong River Basin Precipitation and Elevation Analysis with xarray and geopandas

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变量的纬度变化曲线
axes[0].plot(pr_mean_avg.coords["lat"], mekong_clip_avg, color='black')
axes[1].plot(pr_frequency_avg.coords["lat"], mekong_clip_avg, color='black')
axes[2].plot(pr_intensity_avg.coords["lat"], mekong_clip_avg, color='black')

# 添加高程刻度
axes2 = axes[0].twinx()
axes2.plot([], [])
axes2.set_ylabel("Elevation/m", fontsize=12)
axes2.set_ylim([0, 500])
axes2.spines['right'].set_position(('outward', 0))

plt.tight_layout()
plt.show()

This script analyzes precipitation data from BCC and ERA5 datasets and elevation data from a DEM file within the Mekong River Basin using xarray and geopandas. It visualizes the latitude-dependent trends in precipitation average, frequency, and intensity, along with elevation.

Key Features:

• Loads precipitation data from BCC and ERA5 datasets using xarray.
• Loads elevation data from a DEM file using xarray.
• Loads shapefile for the Mekong River Basin using geopandas to define the analysis area.
• Calculates the average precipitation and elevation values for each latitude within the analysis area.
• Creates three subplots to visualize precipitation average, frequency, and intensity against latitude.
• Adds a secondary y-axis to the first subplot to display elevation values.
• Uses FormatStrFormatter to format latitude axis labels.
• Adjusts the y-axis limits to optimize visualization.

This script provides a comprehensive analysis of precipitation and elevation within the Mekong River Basin, offering insights into the relationship between these variables. Feel free to modify the script to explore different aspects of the data or to apply the same techniques to other regions and datasets.