代码收藏家技术教程 2024-10-28

Python气象绘图——绘制海温填色图

一、作业要求

读入月平均的海温数据，用填色图绘制带地图底图2023年冬季平均的SST与冬季平均气候态SST的差值图，注意设置填色的范围和间隔，用不同色系区分冷暖异常，并确保0值附近为白色，给出colorbar，地图的中心经度设置为180度，冬季平均气候态可选1991-2020年冬季（DJF）的平均。

二、数据来源

数据是ERA5的全球月平均海温数据（1940-01至2023-07），ERA5数据是由欧洲中期天气预报中心（ECMWF）开发的一种全球再分析数据集，粗浅的解释，再分析数据的意思是该数据集不是实际测量的数据，而是通过数值模式运行出来，用来模拟真实大气机制的数据，所以其所描述的大气与真实大气有一定的差别。

下面是下载网站，有很多类型的数据，可以按照自己的需求下载，下载攻略就不再赘述。
Search results (copernicus.eu)https://cds.climate.copernicus.eu/cdsapp#!/search?type=dataset

三、代码

1、引入库

import numpy as np
import matplotlib.pyplot as plt
import xarray as xr
import cartopy.crs as ccrs
import cartopy.feature as cfeature
import cartopy.mpl.ticker as cticker

2、读取数据以及数据处理

filename = r'ERA5_sst_1940_202307.nc'
f1 = xr.open_dataset(filename)
print(f1) # 打印出来查看数据信息
sst = f1.sst.loc[f1.time.dt.month.isin([1,2,12])].loc['1991-01-01':'2020-12-01'] # 读取从1991年到2020年的冬季（12，1，2月）海温
sst_23 = f1.sst.loc[f1.time.dt.month.isin([1,2,12])].loc['2022-12-01':'2023-07-01'] # 读取23年冬季海温，即22年12月，23年1、2月。
lon = sst.longitude # 读取经度
lat = sst.latitude # 读取纬度

打印结果：可以看到数据的时间范围，经纬度格点划分，变量名等。由此来读取我们需要的数据。

aversst = np.zeros((181,360))  #冬季气候态
aversst_23 = np.zeros((181,360))  #23年冬季平均
anosst = np.zeros((181,360)) # 差值
for i in range(90):
    aversst = aversst + sst.values[i,:,:]
aversst = aversst/90.
for i in range(3):
    aversst_23 = aversst_23 +sst_23.values[i,:,:]
aversst_23 = aversst_23/3.
anosst = aversst_23 - aversst

3、绘图

# 检验数据读取绘图
fig = plt.figure(figsize=(12, 5))
leftlon, rightlon, lowerlat, upperlat = (-180, 180, -90, 90)
proj = ccrs.PlateCarree(central_longitude=180)  # 指定投影为经纬度投影。
ax1 = fig.add_axes([0.1, 0.1, 0.9, 0.8], projection=proj)  # 绘制第一个子图
ax1.set_extent([leftlon, rightlon, lowerlat, upperlat], crs=proj)
# 设置经纬度标签
ax1.set_xticks(np.arange(leftlon, rightlon + 30, 30), crs=proj)
ax1.set_yticks(np.arange(lowerlat, upperlat + 10, 10), crs=proj)
lon_formatter = cticker.LongitudeFormatter()
lat_formatter = cticker.LatitudeFormatter()
ax1.xaxis.set_major_formatter(lon_formatter)
ax1.yaxis.set_major_formatter(lat_formatter)
c = ax1.contourf(lon.values, lat.values, anosst, levels=np.arange(-5,6,0.1),cmap=plt.get_cmap('RdBu_r'),
                 transform=ccrs.PlateCarree(),extend='both')
cbar = plt.colorbar(c, shrink=0.7)
coastline_50m = cfeature.COASTLINE.with_scale('50m')
ax1.add_feature(coastline_50m)
land_facecolor = 'white'
land = cfeature.NaturalEarthFeature('physical', 'land', '50m',facecolor=land_facecolor)
ax1.add_feature(land)
cbar.set_label('Surface Sea Temperature')
plt.show()