Python 温度平流曲线不平滑

dousha12321 · 发表于 2020-7-1 09:40:38

登录后查看更多精彩内容~

您需要登录才可以下载或查看，没有帐号？立即注册

x

用python计算fnl 1*1度的grib2数据温度平流，但是不平滑，求解。图片和脚本附上，谢谢！
#!/usr/bin/env python3
# -*- coding: utf-8 -*-
"""
Created on Sat May  9 14:59:45 2020

@author: yiny
"""

import pygrib
import matplotlib.pyplot as plt
import numpy as np
from mpl_toolkits.basemap import Basemap
import math
data = pygrib.open('fnl_20190510_06_00.grib2')
#data.messages ## 表示文件中总共有多少条数据
data.seek(0)  ##ygrib所提供的处理方式类似二进制数据处理，其提供了一些处理二进制数据的方法，比如.seek .tell .rewind .read等方法。
#.seek  跳到第几条记录；.tell 当前所在记录位置；.rewind  回到第一条记录； .read 读取指定个数记录
#for grb in data:
#print(grb)
print('所有此再分析资料数据列表信息输出完毕！')

print('1.所有')
temp = data.select(name='Temperature')
print(temp)
temp= data.select(name='Temperature',typeOfLevel='surface', level=0)

print('2.某高度层temp')## 指定获取层的类型，可以选择获取指定类型层的数据，这对应了不同压力层的的数据
temp = data.select(name='Temperature')[22]#22-700HPA 25-850
#print(temp)，返回的 temp 变量包含了一些方法和属性，比如 temp.data() 包含了温度数据以及经纬度信息，temp.latl
print(temp.latlons())
print(temp.projparams)
#print(temp.data())
UGRD = data.select(name='U component of wind')[23]#23-700HPA 26-850
VGRD= data.select(name='V component of wind')[23]
lons = temp.data()[2][0,:]
lats = temp.data()[1][:,0]
print('lons个数，lats个数：',len(lons),len(lats))#360，181
tc = np.array(temp.data()[0])
u=np.array(UGRD.data()[0])
v=np.array(VGRD.data()[0])
#print(temp.data()[0])
#print(UGRD.data()[0])
#print(VGRD.data()[0])

print('3.计算')##水平温度平流计算公式为 -(udt/dx+vdt/dy)
##对于离散数据,无法做微分会偏微分运算,采取计算方法中
##提到的差分方法，这里利用前差
##对于level和time的选取,因为本code仅仅做示例
##我们选取level为850hPa,time为1948-01-01
##即level[2];time[0]
##还有一个重要问题:dx为某一纬度,通过2.5经度的距离
##即不同纬度dx不等(考虑球坐标);但是对于dy从赤道到极点是固定的
#--------------------------------------#
#首先计算dy
R=6371*1000#地球半径6371km
theta=math.radians(1)##角度转弧度,数据经纬度间隔为2.5°
dy=R*theta
#print('3.计算',dy)
#计算dx
dx=[]
###r为90N到0的小圆半径
for i in range(len(lats[0:181])):
#print(i,lats)
arg=math.radians(lats)
r=R*math.cos(arg)
#print(r)
dx.append(r*math.radians(1))
dx=np.array(dx)##dx为90N到0每一纬圈经过1°经度的距离差
dx[0]=0;dx[180]=0
#print('3.计算',dx)
##先求dt/dx
###因为是圆形区域,所以前差并不会导致横向数据量减少
dtdx=np.empty((181,360))
for i in range(len(lats)):
for j in range(len(lons)):
      #print('4.绘图', i,j,tc[i,j])
      if (i==0) or (i==180):
         dtdx[i,j]=0
      else:
         if j==359:
            dtdx[i,0]=(tc[i,0]-tc[i,j])/dx
         else:
            dtdx[i,j+1]=(tc[i,j+1]-tc[i,j])/dx
###再求dt/dy
dtdy=np.zeros((181,360))
for j in range(len(lons)):
for i in range(len(lats)):
      if i==180:
         dtdy[i,j]=0
      else:
         dtdy[i+1,j]=(tc[i+1,j]-tc[i,j])/dy
###得出温度平流 -(udt/dx+vdt/dy)
######!!!!!!!唯独87.5往上无数据
######!!!!!!!唯独87.5往上无数据
######!!!!!!!唯独87.5往上无数据
dTdt=np.zeros((181,360))
for i in range(len(lats)):#i:[0,180]
#print(i,lats, dTdt[0,180],dTdt[i-1,180],dTdt[180,180])#dTdt[0,180]和[180,180]=0
for j in range(len(lons)):
      dTdt[i,j]=-(u[i,j]*dtdx[i,j]+v[i,j]*dtdy[i,j])
   # print('1111111111111111111',i,j,dTdt[i,j])
#print('4.绘图',dTdt[0:180,:])
#print('.............................................',lats[0:181])
#print("23222222222222222",dTdt[0:181,:])

print('4.绘图')#以上述读取的温度数据为例，绘制中国及周边地区的温度分布：绘图：dTdt[1:72,:]

lons, lats = np.meshgrid(lons, lats)
fig, ax = plt.subplots(figsize=(12, 9))
m = Basemap(projection=temp.projparams.get('proj'), lon_0=180, llcrnrlat=0, llcrnrlon=70, urcrnrlat=60, urcrnrlon=140,ax=ax)
x, y = m(lons[:],lats[0:181])
#print(x,y)
m.drawcoastlines()
#画线
dTdt[0:181,:]=dTdt[0:181,:]*100000
cset = m.contour(x, y, dTdt[0:181,:],100,cmap=plt.cm.gist_ncar)#gnuplot2
_ = m.ax.set_xlim([70, 140])
_ = m.ax.set_ylim(10, 60)
_ = m.readshapefile('/home/qxt/pylinux/shape/china', 'china', linewidth=1.5)
_ = m.readshapefile('/home/qxt/pylinux/shape/cn_province', 'cn_province', linewidth=1.0)
print("vvvvvvvvvvvvvvv")
_ = m.drawparallels(np.arange(10, 61, 10), labels=[1,0,0,0], linewidth=0.1, fontsize=16, fmt="%.2f", dashes=[2,2])
_ = m.drawmeridians(np.arange(70, 141, 15), labels=[0,0,0,1], linewidth=0.1, fontsize=16, fmt="%.2f", dashes=[2,2])
print("vvvvvvvvvvvvvvv2")
cb=plt.colorbar(cset,pad=0.01, shrink=0.75)
cb.ax.set_ylabel('20190510_06Temperature Advection($\circ$C/s)e-05', fontdict=dict(family='Times New Roman', fontsize=20))
#contour = plt.contour(x, y, dTdt[0:181,:],100,colors='b')
plt.clabel(cset,fontsize=5,colors='k')
plt.show()
print("vvvvvvvvvvvvvvv3")

con = m.contourf(x, y, dTdt[0:181,:],np.arange(-1.1e-04, 9e-05), cmap=plt.cm.RdBu_r)
print('6')
#绘制经纬度网格：

#_ = m.readshapefile('中国行政区_包含南海九段线', 'china', linewidth=1.5)

#parallels = np.arange(0.,90,10.)
#m.drawparallels(parallels,labels=[1,0,0,0],fontsize=10)
#以10度为间隔画出0度到北纬90度纬线，并且在图像左侧设置纬线标签
#meridians = np.arange(180.,360.,10.)
#以10度为间隔画出西经180度到本初子午线经线，并且在图像下侧设置经线标签
#m.drawmeridians(meridians,labels=[0,0,0,1],fontsize=10)

cb = fig.colorbar(con, pad=0.02, shrink=0.95)
cb.set_ticks(np.arange(-1.1e-04, 9e-05, 2e-06))
cb.set_ticklabels(np.arange(-1.1e-04, 9e-05, 2e-06))
cb.ax.tick_params(labelsize=16)  ## 控制colorbar ticklables字体大小
  # set colorbar yaxis label font
plt.show()