极地投影及全球相关

题目

课程作业
题目：选取base points, P(45N,165W),A（65N, 20W）,利用1950-2020年月平均 DJF 500mb 高度场计算Figure 8.4 c,d 中的one-point correlation maps. 注意第一年DJF为1950年12月，1951年1， 2月，这三个月是连续的， 记为1950/51 DJF。计算相关系数前去除每个月的70年的气候平均， 得到每个月的异常，再把DJF连接起来， 每个格点得到70*3=210个月的异常Z500时间序列。

导入相关库

import xarray as xr
import pandas as pd
import numpy as np

# 绘图
import matplotlib.pyplot as plt
import matplotlib.path as mpath
import matplotlib.ticker as mticker
import cartopy.crs as ccrs      # 创建一个投影
import cartopy.feature
from cartopy.util import add_cyclic_point
from cartopy.mpl.ticker import (LongitudeFormatter, LatitudeFormatter)
from cartopy.mpl.ticker import LongitudeFormatter, LatitudeFormatter
import cartopy.mpl.ticker as cticker
from matplotlib.colors import ListedColormap,BoundaryNorm
from cartopy.mpl.gridliner import LONGITUDE_FORMATTER, LATITUDE_FORMATTER
import cmaps

# 导入以下包从而使得可以在jupyter中的一个cell输出多个结果
from IPython.core.interactiveshell import InteractiveShell
InteractiveShell.ast_node_interactivity = "all"

数据处理

import xarray as xr
import pandas as pd
import numpy as np

# 绘图
hgt = xr.open_dataset(r'E:/hgt.mon.mean.nc')
Hgt_12 = hgt.sel(level = 500).isel(time=np.arange(35,864,12))['hgt']
Hgt_1 = hgt.sel(level = 500).isel(time=np.arange(36,865,12))['hgt']
Hgt_2 = hgt.sel(level = 500).isel(time=np.arange(37,866,12))['hgt']
# 减去气候平均态
Hgt_12_anomaly = Hgt_12- Hgt_12.mean(axis=0)
Hgt_1_anomaly = Hgt_1 - Hgt_1.mean(axis=0)
Hgt_2_anomaly = Hgt_2- Hgt_2.mean(axis=0)
Hgt_anomaly_DJF = xr.merge([Hgt_12_anomaly,Hgt_1_anomaly,Hgt_2_anomaly])
Hgt_P_anomaly = xr.merge([Hgt_12_anomaly,Hgt_1_anomaly,Hgt_2_anomaly]).sel(lat = 45,lon = 195)
Hgt_A_anomaly = xr.merge([Hgt_12_anomaly,Hgt_1_anomaly,Hgt_2_anomaly]).sel(lat = 65,lon = 340)
Hgt_P_corr = np.zeros([hgt['hgt'].shape[2],hgt['hgt'].shape[3]])
Hgt_A_corr = np.zeros([hgt['hgt'].shape[2],hgt['hgt'].shape[3]])

for lat in range(0,hgt['hgt'].shape[2]):
    for lon in range(0,hgt['hgt'].shape[3]):
        Hgt_P_corr[lat,lon] = np.corrcoef(Hgt_P_anomaly['hgt'].values,Hgt_anomaly_DJF.isel(lat=lat,lon=lon)['hgt'].values)[0,1]


for lat in range(0,hgt['hgt'].shape[2]):
    for lon in range(0,hgt['hgt'].shape[3]):
        Hgt_A_corr[lat,lon] = np.corrcoef(Hgt_A_anomaly['hgt'].values,Hgt_anomaly_DJF.isel(lat=lat,lon=lon)['hgt'].values)[0,1]

绘图

import xarray as xr
import pandas as pd
import numpy as np

# 绘图
# 设置画布
fig = plt.figure(figsize=(16, 10))
# 选定经纬度
lat = np.array(hgt['lat'].values)
lon = np.array(hgt['lon'].values)
# 创建一个投影
proj = ccrs.NorthPolarStereo(central_longitude=90)              # 极地投影
leftlon, rightlon, lowerlat, upperlat = (-180,180.1,15,90)      #仅画15°E-90°E部分
img_extent = [leftlon, rightlon, lowerlat, upperlat]
ax = fig.add_axes([0.1, 0.1, 0.5, 0.5],projection = ccrs.NorthPolarStereo())
ax.gridlines()
ax.set_extent(img_extent, ccrs.PlateCarree())#通过圆柱投影的范围限制地图范围，这样设置地图参数较为方便
ax.add_feature(cfeature.COASTLINE.with_scale('50m'))
# 画圆所需圆框
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)#设置axes边界，为圆形边界，否则为正方形的极地投影
# 经纬网格
gl = ax.gridlines(draw_labels=True, color='grey', alpha=0.5, linestyle='--',      # draw_labels=True表示显示经度标签，color表示内部虚线的颜色
                      xlocs=np.arange(-180, 181, 30), ylocs=[],
                      x_inline=False)
gl.xformatter = LongitudeFormatter()  # 使横坐标转化为经纬度形式
gl.yformatter = LatitudeFormatter()
#################以下步骤添加数据循环，防止白条##################
# 消除环形白带
cycleHgt500, cycle_lon = add_cyclic_point(np.array(Hgt_P_corr), coord=lon)
cycle_LON, cycle_LAT = np.meshgrid(cycle_lon, lat)
c = ax.contour(cycle_LON,cycle_LAT,cycleHgt500, zorder=0, levels =np.arange(-1,1,0.1) , extend = 'both',transform=ccrs.PlateCarree(),cmap =cmap1)
#a = plt.colorbar(c)
b = plt.text(340,65,'A',transform=ccrs.PlateCarree())
#plt.savefig('A.jpg',dpi=500)     
# 注释：对点P同样此操作

one-point correlation maps，A点相关，红色表示正相关，蓝色表示负相关