短期气候Python绘图——平均高度图、距平图、纬偏图

序言:

今天开个新坑，是针对短期气候方面的一些绘图总结，均来自作者本人平时课程实践作业。

一、要求

(1)计算1948-2007年（60年）1月的平均高度场，绘制环流平均图；

(2)计算2008年1月的高度距平，绘制高度距平场图（相对于1948-2007年共60年的平均）；

(3)计算2008年1月的高度场纬偏值，绘制环流纬偏图

二、资料说明

NCEP/NCAR 1948－2012年1～12月的500hPa月平均高度场资料，范围（90°S-90°N,0-360°E)，网格距2.5°×2.5°,纬向格点数144，经向格点数73，资料为GRD格式，资料从南到北、自西向东排列，每月为一个记录，按年逐月排放。（注：数据为二进制）

三、方法说明

1、计算（1948-2007）1月平均高度场，计算平均：

 (i,j)为经、纬格点m为样本长度，k为年份。

2、计算2008年1月的高度距平，计算距平：

(i,j)为经、纬格点m为样本长度，k为年份 。

3、计算2008年1月的高度场纬偏值，计算纬偏：

(i,j)为经、纬格点m为纬圈格点数，k为年份。

四、代码

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

from pylab import *
mpl.rcParams['font.sans-serif'] = ['SimHei']
plt.rcParams['axes.unicode_minus'] = False

# 数据读取与处理
filename = 'hgt500.grd'
data = np.fromfile(filename,dtype=np.float32) # 使用np.fromfile读取二进制文件，数据类型为32位浮点数
data = data.reshape((780, 73, 144)) # 使用reshape函数将读取的数据重构为三维，注意第一维是时间，第二维是纬度，第三维是经度
# 纬度顺序如上安排是由于数据排列顺序与reshape函数特性决定的，在此不再赘述。
data = np.array(data) # 将数据转为数组
janmhgt = np.zeros((73,144)) # 1948-2007年1月平均高度场
anohgt = np.zeros((73,144)) # 2008年高度距平场
anohgt[:,:] = data[732,:,:] # 2008年高度场
h08min = int(np.min(data[732,:,:]))
h08max = int(np.max(data[732,:,:]))
# 平均场
for i in range(0,720,12):
    janmhgt = janmhgt + data[i,:,:]
janmhgt = janmhgt/60.
mmin = int(np.min(janmhgt))
mmax = int(np.max(janmhgt))
# 距平场
anohgt = anohgt - janmhgt
amin = int(np.min(anohgt))
amax = int(np.max(anohgt))
#纬偏
mlat = np.zeros(73) # 2008年纬度平均
devlat = np.zeros((73, 144)) # 2008年纬偏
for i in range(73):
    sumlat = 0
    for j in range(144):
        sumlat = sumlat+data[732,i,j]
    mlat[i] = sumlat/144.
for i in range(73):
    for j in range(144):
        devlat[i,j] = data[732,i,j] - mlat[i]
lmin = int(np.min(devlat))
lmax = int(np.max(devlat))
#print(mlat[57])
lon = np.arange(0, 360, 2.5)
lat = np.arange(-90, 92.5, 2.5)

# 设置绘图函数，参数分别为高度场数据，绘图名称，等值线下限，等值线上限，通过设置等值线上下限来限制绘图的数值范围以及设置等值线间隔
def draw(hgt, name, m, n):
    # 绘制地图底图
    fig = plt.figure(figsize=(10, 6))
    proj = ccrs.PlateCarree(central_longitude=180)
    leftlon, rightlon, lowerlat, upperlat = (0, 360, -90, 90)
    img_extent = [leftlon, rightlon, lowerlat, upperlat]

    ax = fig.add_axes([0.06, 0., 0.9, 1], projection=proj)
    c = ax.contour(lon, lat, hgt,levels=np.arange(m, n, 40),transform=ccrs.PlateCarree(), cmap=jet())
    ax.add_feature(cfeature.COASTLINE.with_scale('50m'))
    ax.add_feature(cfeature.LAKES)
    ax.set_extent(img_extent, crs=ccrs.PlateCarree())
    ax.set_xticks(np.arange(leftlon, rightlon + 20, 20), crs=ccrs.PlateCarree())
    ax.set_yticks(np.arange(lowerlat, upperlat + 20, 20), crs=ccrs.PlateCarree())
    lon_formatter = cticker.LongitudeFormatter()
    lat_formatter = cticker.LatitudeFormatter()
    ax.xaxis.set_major_formatter(lon_formatter)
    ax.yaxis.set_major_formatter(lat_formatter)
    ax.set_title(name, loc='center', fontsize=16)
    xlocs = np.arange(leftlon, rightlon, 20)
    ylocs = np.arange(lowerlat, upperlat, 20)
    #plt.colorbar(c, shrink = 0.4, pad = 0.06)
    plt.clabel(c, fontsize=8)
    #plt.savefig(r'C:\Users\HUAWEI\Desktop\700hPa露点温度差20时.png')
    plt.show()

# 引用函数绘图
draw(janmhgt, 'Jan average 500hPa geopotential height from 1948 to 2007(unit:gpm)', mmin, mmax)
draw(anohgt, 'The 500 hPa geopotential height anomalies in Jan 2008(unit:gpm)', amin, amax)
draw(devlat, 'The latitude deviations of 500hPa height in Jan 2008(unit:gpm)', lmin, lmax)

五、绘图结果

六、结语

从绘制平均高度场图、距平图以及纬偏图，我们可以分析得到许多有用的信息，比如平均图冬季北半球500hPa高度场明显的三槽三脊结构；距平图则是能分析出2008年500hPa一些槽脊系统的强度相较于气候态的变化，偏强抑或是偏弱，或者没有明显变化；纬偏图则是能很方便的分析槽脊系统，极值区域对应着槽脊系统，极大值区域对应有脊，极小值区域对应有槽。通过以上分析我们能够得到500hPa高度场的气候态信息以及目前形势相对于气候态的变化。同时我们也可以用来绘制其他高度或者其他要素，不一定局限在500hPa或是高度场。

由于本人水平限制，多有错误，欢迎批评指正！