首页
学习
活动
专区
圈层
工具
发布
首页
学习
活动
专区
圈层
工具
MCP广场
社区首页 >专栏 >推荐|Python绘制3D动态对流层顶(Tropopause)

推荐|Python绘制3D动态对流层顶(Tropopause)

作者头像
气象学家
发布于 2020-02-17 09:30:27
发布于 2020-02-17 09:30:27
1.6K00
代码可运行
举报
文章被收录于专栏:气象学家气象学家
运行总次数:0
代码可运行

作者

Mathew Barlow: Professor of Climate Science University of Massachusetts Lowell

工具

GFS, the nomads server, python, and the python packages numpy, matplotlib, cartopy, scipy, and netcdf4

potential-vorticity: Python code (gfs_pv_1.2.py) for dynamic tropopause (DT) calculations: DT pressure, DT potential temperature (theta), PV on the 330K isentropic surface, and a PV and theta cross-section at the latitude where the tropopause is lowest in the domain. The date and time need to be set in the beginning of the code; the domain can be changed there as well. gfs_pv_1.2_3D.py is the same code but with additional (poor) 3D plots This code has a DOI and is citable: DOI The data source is the online GFS analysis, and the date and time need to be set within the period of available data. The program can take a few minutes to run because it is accessing data over the internet.

代码

https://github.com/mathewbarlow/potential-vorticity

具体参考以上链接

代码语言:javascript
代码运行次数:0
运行
AI代码解释
复制
#
# run on python 3.7
#
# python code for some calculations related to the dynamic tropopause (DT)  -
# DT pressure, DT potential temperature, 330K PV, 
# and a cross-section of PV at the latitude where the tropopause is lowest -
# all based on the GFS analysis available online.  As the data is accessed
# online, the program can take a while to run.
#
# the date and lat-lon range can be set below
#
# (poorly) coded by Mathew Barlow
# initial release: 14 Nov 2017
# last updated: 10 Oct 2019
#
# this code has *not* been extensively tested and has been 
# awkwardly translated from other coding languages, so if you find
# any errors or have any suggestions or improvements, including for
# the plotting, please let me know at Mathew_Barlow@uml.edu . Thanks!
#
# Support from NSF AGS-1623912 is gratefully acknowledged
#

import numpy as np
import netCDF4
import matplotlib.pyplot as plt
import matplotlib.ticker as tick
from mpl_toolkits.mplot3d import axes3d
import cartopy.crs as ccrs
from scipy.ndimage import gaussian_filter
from cartopy.feature import NaturalEarthFeature
from cartopy.mpl.gridliner import LONGITUDE_FORMATTER, LATITUDE_FORMATTER
from cartopy.mpl.ticker import LatitudeFormatter, LongitudeFormatter

from datetime import datetime


# VALUES TO SET *************************************************
# set date, lat-lon range, and PV-value definition of tropopause
mydate='20191016'
myhour='06'
(lat1,lat2)=(20,70)
(lon1,lon2)=(120,240.1)
tpdef=2   # definition of tropopause in PVU
#****************************************************************


#constants
re=6.37e6
g=9.81
cp=1004.5
r=2*cp/7
kap=r/cp
omega=7.292e-5
pi=3.14159265

# open dataset, retreive variables, close dataset

url='https://nomads.ncep.noaa.gov/dods/gfs_0p25/gfs'+\
mydate+'/gfs_0p25_'+myhour+'z_anl'

file = netCDF4.Dataset(url)

lat_in  = file.variables['lat'][:]
lon_in  = file.variables['lon'][:]
lev = file.variables['lev'][:]

pres2pv_in = file.variables['pres2pv'][0,:,:]
pressfc_in = file.variables['pressfc'][0,:,:]

nlev = lev.size
nx = lon_in.size
ny = lat_in.size

u_in = np.full((nlev, ny, nx), None)
v_in = np.full((nlev, ny, nx), None)
t_in = np.full((nlev, ny, nx), None)
hgt_in = np.full((nlev, ny, nx), None)

ilev = 0
while ilev < nlev:
    print(ilev)
    u_in[ilev, :, :] = file.variables['ugrdprs'][0, ilev, :, :]
    ilev = ilev + 1

ilev = 0
while ilev < nlev:
    v_in[ilev, :, :] = file.variables['vgrdprs'][0, ilev, :, :]
    ilev = ilev + 1

ilev = 0
while ilev < nlev:
    t_in[ilev, :, :] = file.variables['tmpprs'][0, ilev, :, :]
    ilev = ilev + 1

ilev = 0
while ilev < nlev:
    hgt_in[ilev, :, :] = file.variables['hgtprs'][0, ilev, :, :]
    ilev = ilev + 1

#t_in = file.variables['tmpprs'][0,:,:,:]
#u_in = file.variables['ugrdprs'][0,:,:,:]
#v_in = file.variables['vgrdprs'][0,:,:,:]
#hgt_in = file.variables['hgtprs'][0,:,:,:]

file.close()

# get array indices for lat-lon range
# specified above
iy1 = np.argmin( np.abs( lat_in - lat1 ) )
iy2 = np.argmin( np.abs( lat_in - lat2 ) ) 
ix1 = np.argmin( np.abs( lon_in - lon1 ) )
ix2 = np.argmin( np.abs( lon_in - lon2 ) )  

# select specified lat-lon range
t=t_in[:,iy1:iy2,ix1:ix2]
lon=lon_in[ix1:ix2]
lat=lat_in[iy1:iy2]
u=u_in[:,iy1:iy2,ix1:ix2]
v=v_in[:,iy1:iy2,ix1:ix2]
hgt=hgt_in[:,iy1:iy2,ix1:ix2]
pres2pv=pres2pv_in[iy1:iy2,ix1:ix2]
pressfc=pressfc_in[iy1:iy2,ix1:ix2]

# some prep work for derivatives
xlon,ylat=np.meshgrid(lon,lat)

# define potential temperature and Coriolis parameter
theta=t*(1.E5/(lev[:,np.newaxis,np.newaxis]*100))**kap
f=2*omega*np.sin(ylat*pi/180)

lon = np.array(lon, dtype='float')
lat = np.array(lat, dtype='float')
lev = np.array(lev, dtype='float')
u = np.array(u, dtype='float')
v = np.array(v, dtype='float')
hgt = np.array(hgt, dtype='float')
pres2pv = np.array(pres2pv, dtype='float')
pressfc = np.array(pressfc, dtype='float')
theta = np.array(theta, dtype='float')
f = np.array(f, dtype='float')

# calculate derivatives

def ddp(f):
# handle unevenly-spaced levels with 2nd order
# Lagrange interpolation
# except for top and bottom, where use forward diff
    lev3=lev.reshape(lev.size,1,1)*100
    dpp=lev3-np.roll(lev3,-1,axis=0)
    dpm=lev3-np.roll(lev3,1,axis=0)
    fp=np.roll(f,-1,axis=0)
    fm=np.roll(f,1,axis=0)
    ddp_f=(
        fm*dpp/( (dpp-dpm)*(-dpm) ) +
        f*(dpp+dpm)/( dpm*dpp ) +
        fp*dpm/( (dpm-dpp)*(-dpp) )
        )
    ddp_f[0,:,:]=(f[1,:,:]-f[0,:,:])/(lev3[1,:,:]-lev3[0,:,:])
    ddp_f[-1,:,:]=(f[-1,:,:]-f[-2,:,:])/(lev3[-2,:,:]-lev3[-1,:,:])
    return(ddp_f)

def ddx(f):
# use center-difference, assuming evenly spaced lon
# except for side-boundaries, where use forward diff
    x=(re*np.cos(ylat*np.pi/180)*np.pi/180)*lon
    x3=x.reshape(1,x.shape[0],x.shape[1])
    dx3=np.roll(x3,-1,axis=2)-np.roll(x3,1,axis=2)
    ddx_f=(np.roll(f,-1,axis=2)-np.roll(f,1,axis=2))/dx3
    ddx_f[:,:,0]=(f[:,:,1]-f[:,:,0])/(x3[:,:,1]-x3[:,:,0])
    ddx_f[:,:,-1]=(f[:,:,-2]-f[:,:,-1])/(x3[:,:,-2]-x3[:,:,-1])
    return(ddx_f)

def ddy(f):
# use center-difference, assuming evenly spaced lon
# except for N/S boundaries, where use forward diff
    y=(re*np.pi/180)*lat
    y3=y.reshape(1,y.shape[0],1)
    dy3=np.roll(y3,-1,axis=1)-np.roll(y3,1,axis=1)
    ddy_f=(np.roll(f,-1,axis=1)-np.roll(f,1,axis=1))/dy3
    ddy_f[:,0,:]=(f[:,1,:]-f[:,0,:])/(y3[:,1,:]-y3[:,0,:])
    ddy_f[:,-1,:]=(f[:,-2,:]-f[:,-1,:])/(y3[:,-2,:]-y3[:,-1,:])
    return(ddy_f)


#lev3=lev.reshape(lev.size,1,1)
#ddp_theta=np.gradient(theta,lev3*100,axis=0)
#ddx_theta=np.gradient(theta,axis=2)/dx
#ddy_theta=np.gradient(theta,axis=1)/dy

gf=1

ddp_theta=ddp(theta)
ddp_u=ddp(gaussian_filter(u,sigma=gf))
ddp_v=ddp(gaussian_filter(v,sigma=gf))

ddx_theta=ddx(theta)
ddy_theta=ddy(theta)
ddx_v=ddx(gaussian_filter(v,sigma=gf))
ddy_ucos=ddy(gaussian_filter(u,sigma=gf)*np.cos(ylat*pi/180))

# calculate contributions to PV and PV
absvort=ddx_v-(1/np.cos(ylat*pi/180))*ddy_ucos+f
pv_one=g*absvort*(-ddp_theta)
pv_two=g*(ddp_v*ddx_theta-ddp_u*ddy_theta)
pv=pv_one+pv_two

# calculate pressure of tropopause, Fortran-style (alas!)
# as well as potential temperature (theta) and height
#
# starting from 10hPa and working down, to avoid
# more complicated vertical structure higher up
#
nx=ix2-ix1+1
ny=iy2-iy1+1
nz=lev.size
nzs=np.argwhere(lev==50.0)[0,0]
tp=np.empty((ny-1,nx-1))*np.nan   # initialize as undef
tp_theta=np.empty((ny-1,nx-1))*np.nan   # initialize as undef
tp_hgt=np.empty((ny-1,nx-1))*np.nan   # initialize as undef

for ix in range(0,nx-1):
    for iy in range(0,ny-1):
        for iz in range(nzs,0,-1):
            if pv[iz,iy,ix]/1e-6<=tpdef:
                if np.isnan(tp[iy,ix]):
                    tp[iy,ix]=(
                    (lev[iz]*(pv[iz+1,iy,ix]-tpdef*1e-6)
                    -lev[iz+1]*(pv[iz,iy,ix]-tpdef*1e-6))/
                    (pv[iz+1,iy,ix]-pv[iz,iy,ix])
                    )

                    tp_theta[iy,ix]=(
                    ((lev[iz]-tp[iy,ix])*theta[iz+1,iy,ix]+
                    (tp[iy,ix]-lev[iz+1])*theta[iz,iy,ix])/
                    (lev[iz]-lev[iz+1])
                    )

                    tp_hgt[iy,ix]=(
                    ((lev[iz]-tp[iy,ix])*hgt[iz+1,iy,ix]+
                    (tp[iy,ix]-lev[iz+1])*hgt[iz,iy,ix])/
                    (lev[iz]-lev[iz+1])
                    )

# calculate PV on the 330K isentropic surface
# (also not in a pythonic way)
nx=ix2-ix1+1
ny=iy2-iy1+1
nz=lev.size
pv330=np.empty((ny-1,nx-1))*np.nan   # initialize as undef
for ix in range(0,nx-1):
    for iy in range(0,ny-1):
        for iz in range(nz-2,0,-1):
            if theta[iz,iy,ix]>=330:
                if theta[iz-1,iy,ix]<=330:
                    if np.isnan(pv330[iy,ix]):
                        pv330[iy,ix]=(
                        ((330-theta[iz-1,iy,ix])*pv[iz,iy,ix]+
                        (theta[iz,iy,ix]-330)*pv[iz-1,iy,ix])/
                        (theta[iz,iy,ix]-theta[iz-1,iy,ix])
                        )


# slight smoothing of result
# (appears to work better than smoothing u,v,t first)
tp=gaussian_filter(tp,sigma=1)
tp_theta=gaussian_filter(tp_theta,sigma=1)
pv330=gaussian_filter(pv330,sigma=1)

# define spatial correlation function for testing results
def scorr(a,b):
    abar=np.mean(a)
    bbar=np.mean(b)
    covar=sum((a-abar)*(b-bbar))
    avar=sum((a-abar)**2)
    bvar=sum((b-bbar)**2)
    r=covar/np.sqrt(avar*bvar)
    return(r)

# identify latitude of lowest tropopause
maxloc=np.argwhere(tp==np.amax(tp))
latmax=lat[maxloc[0,0]]


# now make some plots - these badly need to be improved

states = NaturalEarthFeature(category='cultural', 
    scale='50m', facecolor='none', 
    name='admin_1_states_provinces_shp')

# get date for plotting
fdate=datetime.strptime(mydate, '%Y%m%d').strftime('%d %b %Y')

plt.close(fig='all')

print('got here')

nframe=30
iframe=0
while iframe<=nframe:
    plt.figure(iframe,figsize=plt.figaspect(0.5))

    pressfc_smooth=gaussian_filter(pressfc,sigma=1)
    ax=plt.gca(projection='3d')

    surf=ax.plot_surface(xlon,ylat,tp,cmap="coolwarm",alpha=1,
                       rstride=1,cstride=1,
                       linewidth=0, antialiased=False)

    ax.plot_surface(xlon,ylat,pressfc_smooth/100,color="lightgray",
                       rstride=1,cstride=1,
                       linewidth=0, antialiased=False)

    ax.set_zlim(1000,100)
    ax.set_xlim(lon1,lon2)
    ax.set_ylim(lat1,lat2)
    ax.view_init(elev=90 - iframe*90/nframe,azim=-90)

    plt.title('2PVU Dynamic Tropopause over topography\n'+myhour+'Z '+fdate)
    plt.colorbar(surf)

    plt.savefig('goo'+ '{:04d}'.format(iframe)+'.png', bbox_inches='tight',
                    dpi=300)

    iframe=iframe+1

彩蛋

1.thorpe_tropopause

https://github.com/mathewbarlow/thorpe_tropopause

2.stratospheric-polar-vortex

https://github.com/mathewbarlow/stratospheric-polar-vortex

本文参与 腾讯云自媒体同步曝光计划,分享自微信公众号。
原始发表:2020-02-10,如有侵权请联系 cloudcommunity@tencent.com 删除

本文分享自 气象学家 微信公众号,前往查看

如有侵权,请联系 cloudcommunity@tencent.com 删除。

本文参与 腾讯云自媒体同步曝光计划  ,欢迎热爱写作的你一起参与!

评论
登录后参与评论
暂无评论
推荐阅读
编辑精选文章
换一批
Python绘制气象实用地图(附代码和测试数据)
前面的推文对于常用的Python绘图工具都有了一些介绍,在这里就不赘述了。本文主要就以下几个方面:“中国区域绘图”、“包含南海”、“兰伯特投影带经纬度标签”、“基于salem的mask方法”、“进阶中国区域mask方法”、“进阶省份mask方法”。对日常的实用需求能够在一定程度上满足。后续就Python在气象常用的统计方法(显著性检验)、合成分析、多变量叠加绘图再进行推送,敬请期待!
MeteoAI
2019/08/12
16.3K3
Python绘制气象实用地图(附代码和测试数据)
NumPyML 源码解析(四)
The wrappers.py module implements wrappers for the layers in layers.py. It includes
ApacheCN_飞龙
2024/02/17
4090
PbootCMS 后台登录界面“3D 云”背景
1.找到后台登录模板文件:/apps/admin/view/default/index.html
Savalone
2020/04/21
9.8K3
3dreshaper_3d曲面屏幕是什么意思
版权声明:本文内容由互联网用户自发贡献,该文观点仅代表作者本人。本站仅提供信息存储空间服务,不拥有所有权,不承担相关法律责任。如发现本站有涉嫌侵权/违法违规的内容, 请发送邮件至 举报,一经查实,本站将立刻删除。
全栈程序员站长
2022/11/04
4930
3dreshaper_3d曲面屏幕是什么意思
文末彩蛋 | 这个 Request URL 长得好不一样
很明显这不是加密的数据,这是一张图片 base64 后的结果,第一次写爬虫朋友遇到这样的请求,可能需要琢磨一下这是什么东西。
全栈程序员站长
2022/11/04
3510
文末彩蛋 | 这个 Request URL 长得好不一样
微信小程序中-[渲染层网络层错误] pages/card/card.wxss 中的本地资源图片无法通过 WXSS 获取-解决办法
2. base64 将图片转换为base64,转换图片网址,转换后将得到的字符放在background-image属性值原位置处
全栈程序员站长
2022/09/12
1.3K0
微信小程序中-[渲染层网络层错误] pages/card/card.wxss 中的本地资源图片无法通过 WXSS 获取-解决办法
君子不玩物丧志,亦常以借物调心,网站集成二次元网页小组件(widget)石蒜模拟器,聊以赏玩
    传世经典《菜根谭》中有言曰:“徜徉于山林泉石之间,而尘心渐息;夷犹于诗书图画之内,而俗气潜消。故君子虽不玩物丧志,亦常借物调心。”意思是,徜徉在林泉山石之间,能够摒弃杂念,留意诗词歌画之中,可以尽弃俗见。所以说君子虽然不会玩物丧志,也常常要借一些优雅的小物件来调理情绪,二次元网页小组件(widget)就是这样的小物件,功能上无甚大观,却可以博君一晒。
用户9127725
2022/09/29
1.4K0
君子不玩物丧志,亦常以借物调心,网站集成二次元网页小组件(widget)石蒜模拟器,聊以赏玩
java图片转二进制流_java将文件转化成二进制流
二进制流的主要编码格式是base64码。可以在网上找一些在线转base64编码的网站进行尝试转换。
全栈程序员站长
2022/11/18
2.4K0
将图片转换为Base64编码字符串、解析Base64编码字符串后生成图片「建议收藏」
当图片转换为base64编码字符串后,其中包含大量的+号,如果我们将上述base64编码字符串通过网络传输给其他接口,那么服务器在解析数据时会把+号当成连接符,然后自动将+号转换为空格,所以为保证数据的准确性,我们需要将空格转换成+号,转换方法如下:
全栈程序员站长
2022/09/14
1.2K0
【数据分析报告】携程客户分析与流失预测
携程作为中国领先的综合性旅行服务公司,每天向超过2.5亿会员提供全方位的旅行服务,因此每天都会产生海量的用户行为数据,这些数据蕴含着丰富的信息资源。另外,客户是企业的重要资源,也是企业的无形资产,客户的流失,也就意味着资产的流失,因此客户流失率是考量业务成绩的一个非常关键的指标。
全栈程序员站长
2022/11/02
7.4K0
【数据分析报告】携程客户分析与流失预测
移动端页面适配方案(viewport)[通俗易懂]
优点是可以使用px布局,不用额外进行rem或者vw等等单位的换算了 缺点是如果是无滚动条的页面在某些设备上(例如平板这种宽高3比4的,折叠屏8比7的)由于宽高比不同有些区域会被挤到视口之外从而导致一些体验上的问题,不过demo2也给出了解决方案;
全栈程序员站长
2022/09/13
6640
Data URI scheme「建议收藏」
data URI scheme 允许我们使用内联(inline-code)的方式在网页中包含数据,目的是将一些小的数据,直接嵌入到网页中,从而不用再从外部文件载入。常用于将图片嵌入网页。
全栈程序员站长
2022/11/02
6110
base64编码图片 生成图片,返回地址[通俗易懂]
发布者:全栈程序员栈长,转载请注明出处:https://javaforall.cn/158546.html原文链接:https://javaforall.cn
全栈程序员站长
2022/09/14
2.1K0
vim的配置文件_vim编辑文件命令
版权声明:本文内容由互联网用户自发贡献,该文观点仅代表作者本人。本站仅提供信息存储空间服务,不拥有所有权,不承担相关法律责任。如发现本站有涉嫌侵权/违法违规的内容, 请发送邮件至 举报,一经查实,本站将立刻删除。
全栈程序员站长
2022/11/10
1.2K0
AJAX学习(一)AJAX基础
ajax在网络应用开发上运用很广泛,它能够达到局部刷新的效果,也就是页面的某一个组件或功能上进行客户端和服务端的数据交互来实现数据的刷新,而不需要整个页面重载,这样可以提升用户的使用感,缩短等待的时间。 ajax的可以用的地方很多,因此是一个很重要的知识点。所以在此写下有关于我对ajax的学习的感悟和应用的一些实例和大家分享,也希望自己对它能够更加了解
全栈程序员站长
2022/11/02
1.7K0
AJAX学习(一)AJAX基础
web添加图片的代码_html保存图片到本地
转载于:https://www.cnblogs.com/larryzeal/p/5991182.html
全栈程序员站长
2022/11/02
3.1K0
web添加图片的代码_html保存图片到本地
使用javascript实现对于chineseocr的API调用「建议收藏」
发布者:全栈程序员栈长,转载请注明出处:https://javaforall.cn/125721.html原文链接:https://javaforall.cn
全栈程序员站长
2022/07/22
9130
使用javascript实现对于chineseocr的API调用「建议收藏」
Hexo添加樱花动态效果背景
前言 最开始玩Hexo博客的时候光主题就选择了半天,当时很中意Sakura,一部分原因就是自带的樱花🌸动态效果,戳到我了简直。但是苦于搞了好久一直出问题,就放弃了。随后转战Butterfly,同样也时非常中意。 也是偶然看到一篇文章才让我明白,有了各路大神的搬运,🌸动效也不会再是Sakura独占了,哈哈哈哈。 效果可以查看这个小姐姐的博客,我当时就是一眼看到,瞬间爱上!!!! 链接在这里:https://cungudafa.gitee.io/ 操作 很简单,只需在_config.butterfly.yml文
花猪
2022/02/16
4.7K0
base编码器_base100编码
Base64编码 是一种基于 64 个可打印字符来表示二进制数据的方法。目前 Base64 已经成为网络上常见的传输 8 位二进制字节代码的编码方式之一。
全栈程序员站长
2022/11/10
4880
base编码器_base100编码
2021第四届浙江省大学生网络与信息安全竞赛预赛部分Writeup
纯签到题,题目给了一个网址,直接burpsuite抓包,在响应头上拿到flag
全栈程序员站长
2022/09/14
8240
2021第四届浙江省大学生网络与信息安全竞赛预赛部分Writeup
推荐阅读
相关推荐
Python绘制气象实用地图(附代码和测试数据)
更多 >
LV.0
南京大学博士生
交个朋友
加入腾讯云官网粉丝站
蹲全网底价单品 享第一手活动信息
领券
问题归档专栏文章快讯文章归档关键词归档开发者手册归档开发者手册 Section 归档
本文部分代码块支持一键运行,欢迎体验
本文部分代码块支持一键运行,欢迎体验