Python气象绘图——一页多图以及不同投影坐标系

安载寿

已于 2024-12-01 15:00:42 修改

阅读量658

点赞数 10

文章标签： python 开发语言

于 2024-12-01 14:59:24 首次发布

本文链接：https://blog.csdn.net/2202_75559127/article/details/144169715

版权

一、作业要求

一页多图，利用所给数据绘制：1. Robinson投影的全球气温图；2.极地投影的50°N以北区域的位势高度图；3.兰伯特投影的东亚地区(70°E-140°E, 15°N-55°N)气温图；4.PlateCarree投影的欧亚地区(50°W-140°E, 0°-90°N)风场图。需要设置colorbar，经纬度坐标，参考矢量等参数。每个子图用a，b，c，d标注

二、数据来源

数据来源于欧洲中期天气预报中心（ECMWF）的再分析数据（ERA5）。数据包括位势高度，气温以及纬向风和经向风。

数据读取以及数据信息。

# 数据读取
filename = r'ERA5_GP_T_UV.nc'
data = xr.open_dataset(filename)
hgt, temp = data.z, data.t
u, v = data.u, data.v
lon = data.longitude.values
lat = data.latitude.values
print(data)

ECMWF再分析数据的下载网址如下：

https://cds.climate.copernicus.eu/cdsapp#!/search?type=dataset

三、代码

import numpy as np
import matplotlib.pyplot as plt
import xarray as xr
from matplotlib import gridspec
import cartopy.crs as ccrs
import cartopy.feature as cfeature
import cartopy.mpl.ticker as cticker
import matplotlib.path as mpath
from cartopy.mpl.gridliner import LONGITUDE_FORMATTER, LATITUDE_FORMATTER

plt.rcParams['font.sans-serif']= ['Microsoft YaHei'] # 设置“微软雅黑”，图上显示出中文
plt.rcParams['axes.unicode_minus'] = False # 设置中文后，解决坐标轴上负号显示问题

# 数据读取
filename = r'ERA5_GP_T_UV.nc'
data = xr.open_dataset(filename)
hgt, temp = data.z, data.t
u, v = data.u, data.v
lon = data.longitude.values
lat = data.latitude.values
print(data)

# 设置布局
fig = plt.figure(figsize=(12,8))
gs = gridspec.GridSpec(2,4)
gs.update(wspace=0.3, hspace=0.3)
ax = {}

# 1、Robinson投影
leftlon, rightlon, lowerlat, upperlat = (lon[0], lon[-1], lat[0], lat[-1])
proj = ccrs.Robinson(central_longitude=180)
ax['t_r'] = fig.add_subplot(gs[0,:2],projection=proj)
ax['t_r'].set_global()
#ax['t_r'].set_stock_img()
coastline_50m = cfeature.COASTLINE.with_scale('50m')
ax['t_r'].add_feature(coastline_50m)
# 使用gridlines添加刻度标签
gl = ax['t_r'].gridlines(draw_labels=True,
                   #xlocs=np.arange(leftlon, rightlon, 10),
                   #ylocs=np.arange(lowerlat, upperlat, 5),
                   linewidth=1, color='gray', alpha=0.8, linestyle='--',
                   xformatter=LONGITUDE_FORMATTER,
                   yformatter=LATITUDE_FORMATTER)
c = ax['t_r'].contourf(lon, lat,temp[0,:,:]-273.15,cmap=plt.get_cmap('RdBu_r'),
                 transform=ccrs.PlateCarree(),extend='both')
cbar = plt.colorbar(c, shrink=0.5)
ax['t_r'].annotate('(a)', xy=(0.01, 1.1), color='black',xycoords='axes fraction',
                   ha='left', va='top', fontsize=14)

# 2、极射赤道投影
proj =ccrs.NorthPolarStereo(central_longitude=180)
leftlon, rightlon, lowerlat, upperlat = (-180,180,50,90)
ax['z'] = fig.add_subplot(gs[0,2:],projection=proj)
ax['z'].set_extent([leftlon, rightlon, lowerlat, upperlat],crs=ccrs.PlateCarree())
gl = ax['z'].gridlines(draw_labels=True,
                   #xlocs=np.arange(leftlon, rightlon, 40),
                   #ylocs=np.arange(lowerlat, upperlat, 5),
                   linewidth=1, color='gray', alpha=0.8, linestyle='--',
                   xformatter=LONGITUDE_FORMATTER,
                   yformatter=LATITUDE_FORMATTER)
# 裁切圆形边框
theta = np.linspace(0, 2*np.pi, 100) # 利用linspace创建了一个从0到2Π的等差数数列，总共100个点，代表圆上的角度
center, radius = [0.5, 0.505], 0.52 # 圆心（center）被设置为坐标（0.5，0.5），这是相对于Axes的坐标系统，原点在左下角，半径为0.51
verts = np.vstack([np.sin(theta), np.cos(theta)]).T # 计算圆上的点的坐标，分别为y和x坐标，然后通过np.vstack和.T将他们组合为数组
circle = mpath.Path(verts * radius + center) # 将点缩放合适大小，加上圆心移动到正确位置，利用Path类表示前面计算的圆形边界
ax['z'].set_boundary(circle, transform=ax['z'].transAxes) # 设置Axes的边界，通过set_boundary来设置自定义的边界，这里传递了circle Path对象
# transform=ax['z'].transAxes说明圆形边界是相对于Axes的大小和位置
ax['z'].add_feature(coastline_50m)
c1= ax['z'].contourf(lon, lat, hgt[0,:,:]/100.,levels=np.arange(480,600,10),cmap=plt.get_cmap('coolwarm'),
                 transform=ccrs.PlateCarree(),extend='both')
#plt.clabel(c1, fontsize=8)
cbar = plt.colorbar(c1, shrink=0.5)
cbar.set_label('位势高度 unit:hPa', rotation=90)
ax['z'].annotate('(b)', xy=(0.01, 0.95), color='black',xycoords='axes fraction',
                   ha='left', va='top', fontsize=14)

# 3、兰伯特投影
standard_parallels = (30, 60)
leftlon, rightlon, lowerlat, upperlat = (70, 140, 15, 55)  # 更改经纬度范围以适应兰伯特投影
proj = ccrs.LambertConformal(central_longitude=(70 + (140)) / 2, standard_parallels=standard_parallels)  # 设定兰伯特投影
ax['t_l'] = fig.add_subplot(gs[1,:2], projection=proj)
ax['t_l'].set_extent([leftlon, rightlon, lowerlat, upperlat], crs=ccrs.PlateCarree())  # 设置地图范围，使用经纬度CRS

# 使用gridlines添加刻度标签
gl = ax['t_l'].gridlines(draw_labels=True,
                   xlocs=np.arange(leftlon, rightlon + 10, 10),
                   ylocs=np.arange(lowerlat, upperlat + 10, 10),
                   linewidth=1, color='gray', alpha=0.8, linestyle='--',
                   xformatter=LONGITUDE_FORMATTER,
                   yformatter=LATITUDE_FORMATTER)
gl.right_labels = False
ax['t_l'].add_feature(coastline_50m)
lonshape = np.arange(70,140.25,0.25)
latshape = np.arange(55,15-0.25,-0.25)
c3 = ax['t_l'].contourf(lonshape, latshape,temp[0,141:302,281:562]-273.15,cmap=plt.get_cmap('RdBu_r'),
                 transform=ccrs.PlateCarree(),extend='both')
cbar = plt.colorbar(c3, shrink=0.5)
ax['t_l'].annotate('(c)', xy=(0.04, 0.95), color='black',xycoords='axes fraction',
                   ha='left', va='top', fontsize=14)

# 4、PlateCarree投影
proj = ccrs.PlateCarree(central_longitude=35)
leftlon, rightlon, lowerlat, upperlat = (-50, 140, 0, 90)
ax['uv'] = fig.add_subplot(gs[1,2:],projection=proj)
ax['uv'].set_extent([leftlon, rightlon, lowerlat, upperlat],crs=ccrs.PlateCarree())
c4 = ax['uv'].quiver(lon[::10], lat[::10], u.values[0,::10,::10], v.values[0,::10,::10],
               color='navy',angles='xy',scale=300,width=0.002,transform=ccrs.PlateCarree())
#由于风场数据过于密集导致绘制的图像无法呈现出风场矢量，所以横纵坐标以及U、V间隔十个用以绘图。
#ax['uv'].set_title('2020年夏季850hPa风场矢量图',fontsize=15,loc='center')
plt.quiverkey(c4,X=0.98,Y=1.02,U=5,angle=90,labelpos='W',label='Wind:5m/s',color='r',labelcolor='r')

ax['uv'].add_feature(coastline_50m)
ax['uv'].set_xticks(np.arange(leftlon, rightlon+10, 20), crs=ccrs.PlateCarree())
ax['uv'].set_yticks(np.arange(lowerlat, upperlat+10, 10), crs=ccrs.PlateCarree())
lon_formatter = cticker.LongitudeFormatter()
lat_formatter = cticker.LatitudeFormatter()
ax['uv'].xaxis.set_major_formatter(lon_formatter)
ax['uv'].yaxis.set_major_formatter(lat_formatter)
gl = ax['uv'].gridlines(draw_labels=False,
                   xlocs=np.arange(leftlon, rightlon + 20, 20),
                   ylocs=np.arange(lowerlat, upperlat + 20, 20),
                   linewidth=1, color='gray', alpha=0.8, linestyle='--',
                   xformatter=LONGITUDE_FORMATTER,
                   yformatter=LATITUDE_FORMATTER)
ax['uv'].annotate('(d)', xy=(0.04, 0.95), color='black',xycoords='axes fraction',
                   ha='left', va='top', fontsize=14)

plt.show()