请选择 进入手机版 | 继续访问电脑版
爱气象,爱气象家园! 

气象家园

 找回密码
 立即注册

QQ登录

只需一步,快速开始

新浪微博登陆

只需一步, 快速开始

搜索
查看: 10660|回复: 7

[混合编程] Python 温度平流 曲线不平滑

[复制链接]

新浪微博达人勋

发表于 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))#360181
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)
##对于离散数据,无法做微分会偏微分运算,采取计算方法中
##提到的差分方法,这里利用前差
##对于leveltime的选取,因为本code仅仅做示例
##我们选取level850hPa,time1948-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=[]
###r90N0的小圆半径
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)##dx90N0每一纬圈经过经度的距离差
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()

研究区域

研究区域

700hpa

700hpa
密码修改失败请联系微信:mofangbao

新浪微博达人勋

发表于 2020-7-1 12:34:15 | 显示全部楼层
赶脚好复杂,为啥子不用GRADS绘制呢?'set csmooth on'就行
密码修改失败请联系微信:mofangbao
回复 支持 反对

使用道具 举报

新浪微博达人勋

发表于 2020-7-2 22:45:08 | 显示全部楼层
{:5_235:}{:5_235:}
密码修改失败请联系微信:mofangbao
回复 支持 反对

使用道具 举报

新浪微博达人勋

 楼主| 发表于 2020-7-7 22:34:10 | 显示全部楼层
oydl1985 发表于 2020-7-1 12:34
赶脚好复杂,为啥子不用GRADS绘制呢?'set csmooth on'就行

我刚开始学 有些地方很啰嗦 后面实在做不出来就去学grads
密码修改失败请联系微信:mofangbao
回复 支持 反对

使用道具 举报

新浪微博达人勋

发表于 2020-9-15 10:44:05 | 显示全部楼层
楼主找到曲线不平滑的原因了吗,我也遇到了这类问题,想请教一下下
密码修改失败请联系微信:mofangbao
回复 支持 反对

使用道具 举报

新浪微博达人勋

发表于 2020-9-16 11:14:42 | 显示全部楼层
先把原始数据插值插的更密点,再画图就平滑一些了
密码修改失败请联系微信:mofangbao
回复 支持 反对

使用道具 举报

新浪微博达人勋

 楼主| 发表于 2020-9-24 22:08:39 | 显示全部楼层
oydl1985 发表于 2020-7-1 12:34
赶脚好复杂,为啥子不用GRADS绘制呢?'set csmooth on'就行

最后用MI了
密码修改失败请联系微信:mofangbao
回复 支持 反对

使用道具 举报

新浪微博达人勋

 楼主| 发表于 2020-9-24 22:58:22 | 显示全部楼层
小阳光 发表于 2020-9-15 10:44
楼主找到曲线不平滑的原因了吗,我也遇到了这类问题,想请教一下下

最后用MI了 王老师现成的脚本改了下
密码修改失败请联系微信:mofangbao
回复 支持 反对

使用道具 举报

您需要登录后才可以回帖 登录 | 立即注册 新浪微博登陆

本版积分规则

Copyright ©2011-2014 bbs.06climate.com All Rights Reserved.  Powered by Discuz! (京ICP-10201084)

本站信息均由会员发表,不代表气象家园立场,禁止在本站发表与国家法律相抵触言论

快速回复 返回顶部 返回列表