Python计算气候场、距平场、均方差场（后附气象水文科研猫交流群）

import pandas as pd
import numpy as np
import math
import matplotlib.pyplot as plt
import cartopy
import cartopy.crs as ccrs
import cartopy.feature as cfeature
import cartopy.mpl.ticker as cticker
from pylab import *
from matplotlib import rcParams
config = {"font.family":'Times New Roman',"font.size": 12,"mathtext.fontset":'stix'}
rcParams.update(config)
##数组建立
h500=np.zeros((17,37,48))
ave=np.zeros((17,37,12))
juping=np.zeros((17,37,48))
std=np.zeros((17,37,12))
##数据读取
data=pd.read_table('F:/Rpython/lp37/data37/H500.DAT',header=None,skiprows=1,sep='\s+')
data=np.array(data)
for k in range(0,48):
    for i in range(0,17):
        for j in range(0,37):
            h500[i,j,k]=data[i+k*18,j]
##气候平均场
for k in range(0,12):
    for i in range(0,17):
        for j in range(0,37):
            ave[i,j,k]=0.25*(h500[i,j,k]+h500[i,j,k+12]+h500[i,j,k+24]+h500[i,j,k+36])
#距平计算
for k in range(0,48):
    juping[:,:,k]=h500[:,:,k]-ave[:,:,k%12]
##均方差
for k in range(0,12):
    for i in range(0,17):
        for j in range(0,37):
            std[i,j,k]=math.sqrt((1/4)*(juping[i,j,k]**2+juping[i,j,k+12]**2+juping[i,j,k+24]**2+juping[i,j,k+36]**2))
t=ave#画图中转变量
##经纬度
lon=np.arange(60,151,2.5)
lat=np.arange(0,41,2.5)
##气候场绘图（共12张）
for k in range(0,12):
    ave=t
    #建立画布
    proj = ccrs.PlateCarree() # 设置投影
    fig, f2_ax1 = plt.subplots(figsize=(15,15), subplot_kw=dict(projection=proj))
    leftlon, rightlon, lowerlat, upperlat = (60,150,0,40)
    #绘制
    ave=f2_ax1.contour(lon,lat,ave[:,:,k],cmap='gist_rainbow',extend='both')
    #在画布的绝对坐标建立子图
    f2_ax1.set_extent([leftlon, rightlon, lowerlat, upperlat], crs=ccrs.PlateCarree())
    #海岸线，50m精度
    f2_ax1.add_feature(cfeature.COASTLINE.with_scale('50m'))
    #湖泊数据
    f2_ax1.add_feature(cfeature.LAKES, alpha=0.5)
    #地理坐标标签格式
    f2_ax1.set_xticks(np.arange(leftlon,rightlon+10,10), crs=ccrs.PlateCarree())
    f2_ax1.set_yticks(np.arange(lowerlat,upperlat+10,10), crs=ccrs.PlateCarree())
    lon_formatter = cticker.LongitudeFormatter()
    lat_formatter = cticker.LatitudeFormatter()
    f2_ax1.xaxis.set_major_formatter(lon_formatter)
    f2_ax1.yaxis.set_major_formatter(lat_formatter)
    f2_ax1.set_title('jupipng %d mounth'%(k+1),loc='center',fontsize =15)
    f2_ax1.stock_img()
    # shrink 控制 colorbar 长度，pad 控制colorbar和图的距离
    cb = fig.colorbar(ave,shrink=0.42,aspect=30,fraction=.03,pad=0.01)
    fig.savefig('F:/Rpython/lp37/plot13/f1/%d.png'%(k+1),bbox_inches='tight')

##距平场绘图（共48张
for k in range(0,48):
    #建立画布
    proj = ccrs.PlateCarree() # 设置投影
    fig, f2_ax1 = plt.subplots(figsize=(15,15), subplot_kw=dict(projection=proj))
    leftlon, rightlon, lowerlat, upperlat = (60,150,0,40)
    #绘制
    ave=f2_ax1.contourf(lon,lat,juping[:,:,k],cmap='gist_rainbow',extend='both')
    #在画布的绝对坐标建立子图
    f2_ax1.set_extent([leftlon, rightlon, lowerlat, upperlat], crs=ccrs.PlateCarree())
    #海岸线，50m精度
    f2_ax1.add_feature(cfeature.COASTLINE.with_scale('50m'))
    #湖泊数据
    f2_ax1.add_feature(cfeature.LAKES, alpha=0.5)
    #以下6条语句是定义地理坐标标签格式
    f2_ax1.set_xticks(np.arange(leftlon,rightlon+10,10), crs=ccrs.PlateCarree())
    f2_ax1.set_yticks(np.arange(lowerlat,upperlat+10,10), crs=ccrs.PlateCarree())
    lon_formatter = cticker.LongitudeFormatter()
    lat_formatter = cticker.LatitudeFormatter()
    f2_ax1.xaxis.set_major_formatter(lon_formatter)
    f2_ax1.yaxis.set_major_formatter(lat_formatter)
    f2_ax1.set_title('jupipng %.1f'%(1982+int(k/12)+(k+1)*0.1),loc='center',fontsize =15)
    # shrink 控制 colorbar 长度，pad 控制colorbar和图的距离
    cb = fig.colorbar(ave,shrink=0.42,aspect=30,fraction=.03,pad=0.01)
    plt.savefig('F:/Rpython/lp37/plot13/f2/%d.png'%(k+1),bbox_inches='tight')

##均方差绘图（共12张）
for k in range(0,12):
    #建立画布
    proj = ccrs.PlateCarree() # 设置投影
    fig, f2_ax1 = plt.subplots(figsize=(15,15), subplot_kw=dict(projection=proj))
    leftlon, rightlon, lowerlat, upperlat = (60,150,0,40)
    #绘制
    ave=f2_ax1.contourf(lon,lat,std[:,:,k],cmap='gist_rainbow',extend='both')
    #在画布的绝对坐标建立子图
    f2_ax1.set_extent([leftlon, rightlon, lowerlat, upperlat], crs=ccrs.PlateCarree())
    #海岸线，50m精度
    f2_ax1.add_feature(cfeature.COASTLINE.with_scale('50m'))
    #湖泊数据
    f2_ax1.add_feature(cfeature.LAKES, alpha=0.5)
    #以下6条语句是定义地理坐标标签格式
    f2_ax1.set_xticks(np.arange(leftlon,rightlon+10,10), crs=ccrs.PlateCarree())
    f2_ax1.set_yticks(np.arange(lowerlat,upperlat+10,10), crs=ccrs.PlateCarree())
    lon_formatter = cticker.LongitudeFormatter()
    lat_formatter = cticker.LatitudeFormatter()
    f2_ax1.xaxis.set_major_formatter(lon_formatter)
    f2_ax1.yaxis.set_major_formatter(lat_formatter)
    f2_ax1.set_title('jupipng %d mounth'%(k+1),loc='center',fontsize =15)
    # shrink 控制 colorbar 长度，pad 控制colorbar和图的距离
    cb = fig.colorbar(ave,shrink=0.42,aspect=30,fraction=.03,pad=0.01)
    fig.savefig('F:/Rpython/lp37/plot13/f3/%d.png'%(k+1),bbox_inches='tight')