```
import os
import matplotlib.ticker as mticker
import netCDF4 as nc
import matplotlib.path as mpath
import cmaps
import matplotlib.pyplot as plt
import numpy as np
import matplotlib as mpl
import cartopy.crs as ccrs
import cartopy.feature as cfeature
from netCDF4 import Dataset
from cartopy.mpl.gridliner import LONGITUDE_FORMATTER, LATITUDE_FORMATTER
import regionmask
import cartopy.io.shapereader as shpreader
#-------------------------------- 1. 读取NC文件 --------------------------------#
nc_file = r"D:\biyelunwen\data\sat_swe_flow.nc"
ds = Dataset(nc_file)
# 假设数据变量名为'snow',经纬度为'lon'和'lat'
# 请根据实际NC文件结构调整变量名称和维度顺序
lon = ds.variables['longitude'][:] # 读取经度数组
lat = ds.variables['latitude'][:] # 读取纬度数组
data1 = ds.variables['flow_sat_to_swe'][:] # 读取数据,假设为三维数组(time, lat, lon)
data2 = ds.variables['flow_swe_to_sat'][:]
# 示例:取第一个时间步和创建示例数据
dss = np.array([data1, data2]) # 生成两个数据层用于绘图演示
#-------------------------------- 2. 图形参数设置 --------------------------------#
mpl.rcParams["font.family"] = 'Arial'
mpl.rcParams["mathtext.fontset"] = 'cm'
mpl.rcParams["font.size"] = 12
mpl.rcParams["axes.linewidth"] = 1
proj = ccrs.NorthPolarStereo(central_longitude=0) # 北半球极地投影
leftlon, rightlon, lowerlat, upperlat = (-180, 180, 20, 90)
img_extent = [leftlon, rightlon, lowerlat, upperlat]
#-------------------------------- 3. 创建画布和子图 --------------------------------#
fig = plt.figure(figsize=(14,6), dpi=600)
title = np.array(['(a) inf from SAT to SWE', '(b) inf from SWE to SAT']).reshape(2, 1)
for j in range(2):
# 设置子图位置
left = 0.1 + j * 0.25
ax = fig.add_axes([left, 0.5, 0.35, 0.4], projection=proj)
# 添加标题
ax.set_title(f'{title[j, 0]}', loc='center', fontsize=14)
#-------------------------------- 4. 设置地图要素 --------------------------------#
# 添加网格线
gl = ax.gridlines(crs=ccrs.PlateCarree(), draw_labels=True,
linewidth=1, color='grey', linestyle='--',
xlocs=[-180, -90, 0, 90, 180])
gl.xformatter = LONGITUDE_FORMATTER # 经度格式
gl.yformatter = LATITUDE_FORMATTER # 纬度格式
gl.rotate_labels = False # 禁止旋转标签
ax.set_extent(img_extent, ccrs.PlateCarree()) # 设置地图范围
ax.add_feature(cfeature.COASTLINE.with_scale('110m')) # 添加海岸线
#-------------------------------- 5. 裁剪圆形边界 --------------------------------#
theta = np.linspace(0, 2*np.pi, 100)
center, radius = [0.5, 0.5], 0.5
verts = np.vstack([np.sin(theta), np.cos(theta)]).T
circle = mpath.Path(verts * radius + center)
ax.set_boundary(circle, transform=ax.transAxes)
#-------------------------------- 6. 添加青藏高原轮廓 --------------------------------#
tp_shape = shpreader.Reader(r"D:\biyelunwen\TPBoundary_new(2021)\TPBoundary_new(2021).shp").geometries()
ax.add_geometries(tp_shape, crs=ccrs.PlateCarree(),
edgecolor='red', facecolor='none',
linewidth=1.5, linestyle='--', alpha=0.8)
#-------------------------------- 7. 数据处理与可视化 --------------------------------#
# 创建二维经纬度网格
lon2d, lat2d = np.meshgrid(lon, lat)
# 绘制填色图
pcm = ax.contourf(lon, lat, transform=ccrs.PlateCarree(),
cmap='RdBu_r', levels=np.arange(-4, 4, 0.4),
extend='neither')
#-------------------------------- 8. 添加颜色条 --------------------------------#
cax = fig.add_axes([0.2, 0.4, 0.4, 0.03]) # 调整颜色条位置
cbar = fig.colorbar(pcm, cax=cax, orientation='horizontal')
cbar.set_label('Snow Depth (m)') # 设置颜色条标签
plt.show()```如何画出data1,data2的北极极地投影图
最新发布