Rossby于1930s年代提出位涡的想法,是希望得到一个与垂直涡度相关的、类似于位温一样的守恒量。这里的“垂直”是指垂直于分层界面。
这一想法是具有重要而深远的意义的。因为这一概念不仅适用于地球上的大气和海洋,也适用于太阳这一类恒星的辐射内部。它对于理解平衡流以及一系列基本的动力过程都很重要,例如Rossby波的传播与破碎及其对大气的影响、全球尺度的遥相关、急流self-sharpening等的反摩擦现象、气旋反气旋和风暴路径的生成、风的形成等等。
感觉位涡可以理解为等熵面间单位质量的涡度
位涡是一个标量,是一个综合表征大气运动状态和热力状态的物理量,它的重要性在于绝热无摩擦运动中微团的位涡守恒,位涡守恒定理揭示了涡度变化是受到大气热力结构制约的。
位涡是由绝对涡度和静力稳定度两部分组成的,平流层和对流层顶及高空锋区附近的高位涡主要是因为温度垂直梯度很大,和涡度的大小无明显联系。
等熵面就是等位温面,等位温面向低纬倾斜
一般会以2PVU为对流层顶
平流层的位涡比对流层的位涡大很多(因为垂直温度梯度大)。相同高度下,高纬的位涡大于低纬地区(因为f)
位温随高度升高而增大,因此当等位温线上突时,表明此处位温小,对应温度小,偏冷
因为PV和位温都是随着高度的升高而增大。因此等PV面上的低位温表明此处等PV面下凹,有PV大值,对应气旋;高位温表明有反气旋
位涡的特点:
空气微团的位涡具有守恒性,只有在非绝热加热和摩擦的作用下能使位涡发生变化. 位涡守恒状况下,当空气团被压缩时,其绝对涡度减小;当空气团被拉伸时,其绝对涡度增大.可用这个解释山地背风波现象.
同样,当辐散导致空气微团水平膨胀时,绝对位涡减弱.
Rossby的位涡公式(1936年):对于无耗散流,位涡守恒。因此当两个等熵面远离时,绝对涡度的垂直分量通过涡旋拉伸而增加;当等熵面只是倾角发生变化时,绝对涡度的垂直分量守恒。
Ertel的位涡公式(1942年),用了三维的涡度矢量,位温的三维梯度。在静力平衡、垂直运动可以忽略的情况下,Ertel公式与Rossby的位涡公式类似。在准地转位涡等相比较时,这两个公式的计算结果一般都被当成是精确的位涡
Ertel的位涡公式的变形,用于验证位涡在无耗散流中是守恒的
稳定大气的PV逆算子
准地转位涡公式( 最简单但最不精确。首先只保留了水平风速,尽管此时垂直运动很明显;其次,遗弃了有垂直输送和水平输送的精确的Rossby-Ertel PV,引入了一个只有水平输送的准地转PV)
ttest检验原理公式
pytho的显著线检验一般都使用SciPy程序包,当检验两组样本均值时,有如下几种方法:
from scipy import stats
import numpy as np
tvalue,pvalue = stats.ttest_ind_from_stats(mean1, std1, numb1, \
mean2, std2, numb2, equal_var=True, alternative='two-sided')
# 输入两组数据的均值、方差和样本量以进行显著性检验
# equal_var=True表示两组数据具有相近的总体方差(默认)
# 若False,表示两组数据不具有相近的总体方差
# 一样的输入,若equal_var=True,则得到的pvalue偏小,更容易通过显著性检验
tvalue,pvalue = stats.ttest_ind(data1,data2, axis=0, equal_var=True,\
nan_policy='propagate', permutations=None, random_state=None, \
alternative='two-sided', trim=0)
# 也是用于检验两组样本均值是否显著差异,但该函数需要将两个样本全部输入,
# 若样本量较大,则比较占用内存,因此个人还是更喜欢上面那个函数
# 关于各参数的具体介绍可以看
# https://docs.scipy.org/doc/scipy/reference/generated/scipy.stats.ttest_ind.html
# 此外 SciPy 也有提供检验两组样本总体方差是否相等的函数levene与bartlett
# 输出统计值和pvalue,当pvalue大时,说明两组样本的总体方差差异大
# 相同的两组数据,bartlett的pvalue大于levene的
print(stats.bartlett(data1, data2))
print(stats.levene(data1, data2))
siglvl = 0.05 # 95% 显著性水平
var.values=np.ma.array(var.values,mask=(pvalue>siglvl))
mask = np.array([pu>siglvl,pv>siglvl]).all(axis=0)
uwnd.values=np.ma.array(uwnd.values,mask=mask)
vwnd.values=np.ma.array(vwnd.values,mask=mask)
每张画布只有一个Figure,而一个Figure可以包含多个Axes对象。一个axes对象即一个子图,这个对象包含大量的子元素,包括线条、文本、坐标轴复合对象等等。之所以和Figure进行分离的原因是,Figure用来控制尺寸等一些和元素无关的属性。Axes则用来管理内部的元素。通过Figure.add_axes()可以创建一个axes。或者,直接调用fig,ax = plt.subplots(n)一次创建这两个对象。
第三个重要的层是x(y)axis,这个层也是一个复合对象(实际上是mpl.axis.Axis和mpl.axis.Tick两个类),包含X或者Y轴的元素,包括刻度(Ticks)、刻度文字(Tick Labels ; Offset text),轴标签(Axis Label)等。这一层用来控制图形的尺度,比如x/ylim用来控制坐标轴范围。比如xlabel用来设置label的一些属性,比如字体大小、旋转、颜色、背景和透明度等。
当调用plt.plot()方法的时候,如果系统没有检测到存在Figure对象和Axes对象,则会创建一个,如果检测到,则直接使用。
对于想要重新画一幅图,而不是将图形绘制在一个Figure中并且区分成几个subplot,直接使用plt.figure()创建一个新的Figure对象,使用plt.gcf()也可以捕获这个对象。
每个图形可以有很多子图形构成,每个子图形只有一个Axes,它们可以共享坐标轴,每个子图形都有其各自的AxesSubplot子轴,包含在总的Axes中。
使用plt进行绘图固然很方便,但如果需要细微调整,则不太方便。调用axes单个元素可以获得更精细的控制,并且更面向对象一些。
import matplotlib.pyplot as plt
import numpy as np
x = np.arange(0.1, 4, 0.5)
y = np.exp(-x)
#几种创建图形 Fig 和 Axes 的方法
fig = plt.figure(figsize=(12,9),dpi=100)
axe = plt.axes(projection=cart_proj)
fig = plt.figure(figsize=(12,9),dpi=300)
axe = fig.add_axes([0.05, 0.05, 0.9, 0.65]) # 括号中列表指定了该Axes对象的左下角的x、y值、长dx,高dy
fig = plt.figure(figsize=(12,9),dpi=100)
axe = plt.subplot(1,2,i+1,projection=cart_proj) # 创建一行两列子图
fig, ax = plt.subplots(3, 3, sharex=True, sharey=True) # 返回的ax是一个图形列表
fig, ((ax1,ax2,ax3),(ax4,ax5,ax6),(ax7,ax8,ax9)) = plt.subplots(3,3)
print(ax)
# plot the linear_data on the 5th,8th subplot axes
ax[2][2].plot(x,y, '-')
ax[1][1].plot(x,y**2, '-')
# 对于一个包含多个轴的fig,可以使用 gcf().get_axis,或者直接对矩阵进行索引和操纵
for ax in plt.gcf().get_axes():
for label in ax.get_xticklabels() + ax.get_yticklabels():
label.set_visible(True)
# necessary on some systems to update the plot
plt.gcf().canvas.draw()
参考资料:
关于时间的处理,和ncl很不一样的是,Python 提供了一个 time 和 calendar 模块可以用于格式化日期和时间,想处理成任何格式的日期都可以。
python中的日期有两种类型:
import time
# 格式化成 2016-03-20 11:45:39 形式
print (time.strftime("%Y-%m-%d %H:%M:%S", time.localtime()))
# time.localtime([secs])可获得当前时间的本地时间元组,也接受时间戳,并返回当地时间下的时间元组
# time.strftime(fmt[,tupletime])只接收时间元组,并返回以字符串表示的时间
# 将格式字符串转换为时间戳
a = "Sat Mar 28 22:24:24 2016"
print (time.mktime(time.strptime(a,"%a %b %d %H:%M:%S %Y")))
# time.strptime(str,fmt)根据fmt的格式把一个时间字符串解析为时间元组
# time.mktime(tupletime) 接受时间元组并返回时间戳
time.sleep(secs) # 推迟调用线程的运行,secs指秒数。
time.process_time() # 返回当前进程执行 CPU 的时间总和,不包含睡眠时间。由于返回值的基准点是未定义的,所以,只有连续调用的结果之间的差才是有效的。
calendar.isleap(year) # 是闰年返回 True,否则为 false。
calendar.leapdays(y1,y2) # 返回在Y1,Y2两年之间的闰年总数。
pandas的 pd.date_range(start*, end*, periods*, freq)
函数可以生成任意间隔的时间序列。带*的参数表示可选参数
举例如下:
import pandas as pd
ftime = pd.date_range(start='2021-09-25 00',end='2021-09-25 12', freq='1H' \
,closed=None).strftime("%Y-%m-%d_%H:00").to_list()
'''
若没有后面转换格式的命令.strftime("%Y-%m-%d_%H:00").to_list(),得到的结果如下:
DatetimeIndex(['2021-09-25 00:00:00', '2021-09-25 01:00:00',
'2021-09-25 02:00:00', '2021-09-25 03:00:00',
'2021-09-25 04:00:00', '2021-09-25 05:00:00',
'2021-09-25 06:00:00', '2021-09-25 07:00:00',
'2021-09-25 08:00:00', '2021-09-25 09:00:00',
'2021-09-25 10:00:00', '2021-09-25 11:00:00',
'2021-09-25 12:00:00'],
dtype='datetime64[ns]', freq='H')
'''
datetime是Python处理日期和时间的另一个标准库,pd.date_range(start, end, periods*, freq)生成的时间就是这个格式的。 字符串与datetime之间的相互转换和time模块类似
from datetime import datetime, timedelta, timezone
now = datetime.now() # 获取当前datetime
print(now.strftime('%a, %b %d %H:%M'))
# 将datetime以自定义格式转换为字符串打印
dt = datetime.strptime('2015/6/1 18:19:59', '%Y/%m/%d %H:%M:%S')
# 把str转换为datetime,只能一个时间一个时间地转换
dt = datetime(2015, 4, 19, 12, 20) # 用指定日期时间创建datetime
# print得到的结果会是 2015-04-19 12:20:00
print(dt.year)
# 也可以单独查看 month, day, hour, minute, second, microsecond, tzinfo
# 但无法通过dt.hour=13修改datetime的时间,只能通过timedelta对日期时间进行加减,得到新的时间
now = datetime.now() # 获取当前datetime
a1 = now + timedelta(hours=10) #时间向后推移10个小时
a2 = now - timedelta(days=1) #日期向前推移1天
# 也可以通过两个时间点相减,获得总的时间间隔秒数
dt1 = datetime(2015, 4, 19, 12, 20)
dt2 = datetime(2015, 4, 21, 12, 20)
(dt2-dt1).total_seconds()/60.0/60.0 #间隔的小时数
dt.timestamp() # 把datetime转换为timestamp
dt = datetime.fromtimestamp(time.time()) # 把时间戳转为本地时区的datetime
dt = datetime.utcfromtimestamp(time.time()) # 把时间戳转为UTC的datetime
以ERA5为例,其时间坐标是以自1900年1月1日午夜零时以来的小时数表示,和python的时间戳不一样。
int time(time) ;
time:units = "hours since 1900-01-01 00:00:00.0" ;
time:long_name = "time" ;
time:calendar = "gregorian" ;
当用 xarray.open_dataset 打开nc文件读取时间维,好像是会自动将该时间坐标转为datetime.
但若是用netCDF库的Dataset打开nc文件,则需要使用 num2date 将时间坐标转为python时间戳(根据ECMWF给出的python例子),如下所示:
from netCDF4 import Dataset, date2num, num2date
with Dataset(f_in) as ds:
var_time = ds.variables['time']
time_avail = num2date(var_time[:], var_time.units,
calendar = var_time.calendar)
'''
print(var_time) 的结果:
<class 'netCDF4._netCDF4.Variable'>
float64 time(time)
standard_name: time
units: hours since 1994-12-01 00:00:00
calendar: proleptic_gregorian
unlimited dimensions: time
current shape = (10248,)
filling on, default _FillValue of 9.969209968386869e+36 used
print(time_avail)的结果:
[cftime.DatetimeProlepticGregorian(1994, 12, 1, 0, 0, 0, 0, has_year_zero=True)
cftime.DatetimeProlepticGregorian(1994, 12, 1, 1, 0, 0, 0, has_year_zero=True)
...
cftime.DatetimeProlepticGregorian(1996, 1, 31, 22, 0, 0, 0, has_year_zero=True)
cftime.DatetimeProlepticGregorian(1996, 1, 31, 23, 0, 0, 0, has_year_zero=True)]
'''
import xarray as xr
ds = xr.open_dataset("ERA5_NH_t_1989.nc")
t = ds['time']
'''
print(t)的结果:
<xarray.DataArray 'time' (time: 1460)>
array(['1989-01-01T00:00:00.000000000', '1989-01-01T06:00:00.000000000',
'1989-01-01T12:00:00.000000000', ..., '1989-12-31T06:00:00.000000000',
'1989-12-31T12:00:00.000000000', '1989-12-31T18:00:00.000000000'],
dtype='datetime64[ns]')
Coordinates:
* time (time) datetime64[ns] 1989-01-01 ... 1989-12-31T18:00:00
Attributes:
long_name: time
'''
# 挑选特定时间段的数据
# 方法一:使用datetime
from datetime import datetime
strt=["1989120206","1989120312"]
tim=[]
for str1 in strt:
tim.append(datetime.strptime(str1, '%Y%m%d%H'))
var = ds['t'].sel(level=250,time=tim)
# 方法二:使用isin,生成逻辑数组
stime = np.array([ds.time.dt.year.isin([1989,1981]),
ds.time.dt.dayofyear.isin([185,186,187])]).all(axis=0)
var= ds['t'].sel(time=stime)
由于用 xarray 打开 netcdf 文件得到的时间维是一个元素类型为 datetime64[ns] 的 xarray.DataArray 数据,因此无法直接使用 strftime(“%Y-%m-%d_%H:00”) 将其转换其他格式的时间字符串。
此时可以使用 np.datetime_as_string(ds.time) 将其转换为字符串元素的 numpy.ndarray,然后用 pd.to_datetime() 将其转换为 pandas 的数据类型,具体如下:
from datetime import datetime
time = pd.to_datetime(np.datetime_as_string(ds.time,unit='h'),format='%Y-%m-%dT%H')
# 如果只转换一个时刻
str_time = datetime.strptime(np.datetime_as_string(var.time[0]),
'%Y-%m-%dT%H:00:00.000000000').strftime('%b %d')
当时间坐标轴用 datetime.datetime 列表、pd.date_range 生成的数组表示时,plt中均有函数可以设置横纵坐标的时间格式,示例如下:
import matplotlib.dates as mdates
from datetime import datetime
import pandas as pd
ftime = pd.date_range(start='2021-09-16 00',end='2021-09-16 23',freq='1H',closed=None)
ftime = []
for nt in
ftime.append(datetime.strptime(nt,'%Y-%m-%d %H:00:00'))
# datetime一次只能处理一个变量,不能用于列表和数组
axe.plot(ftime, var, label=case[nc], linewidth=2)
axe.tick_params(axis='both', which='major',labelsize=label_font)
axe.xaxis.set_major_formatter(mdates.DateFormatter("%m-%d_%H:00"))
plt.setp(axe.get_xticklabels(), rotation=30, ha="right")
在用shell脚本批处理文本文件时,可能需要读取或修改ascii文件中的一些信息。在shell中有许多可以处理文本文件的工具,配合正则表达式,可以实现很多功能。现罗列如下:
上述命令的一些常见用法:
# 统计文本文件中某单词出现次数
cat tr_trs_pos.txt | grep -o "TRACK_ID" | wc -l
# 竖杠表示管道命令,将前一命令的结果作为下一个命令的输入
# grep -o 表示搜索文本并只显示匹配的部分
# wc [-clw][文件…]可以统计文件(如果没有指定文件,则会从标准输入读取数据)
# 的字符数(-c),字数(-w),行数(-l),这里就是只统计行数
# 统计文件中每一个单词出现的次数,并按照出现频率从大到小排序
sed 's/ /\n/g' file.txt | sort | uniq -c | sort -nr
# 使用sed读取file.txt文件中的所有内容,并将空格替换为换行符,然后将修改后的内容向sort输入
# sed 在不使用 -i 参数时,是不会直接将文件里的内容修改,因此此时file.txt里的内容不会被修改
# sort可以将文本内容加以排序,-n表示按照数值大小排序,-r表示以相反的顺序排序
# uniq 检查并删除文本中重复出现的行列,-c或--count 在每列旁边显示该行重复出现的次数。
# 重复的行并不相邻时,uniq 命令不起作用,因此需要先 sort
sed -i "34s/.*/${fras}/" STATS.latlng_1hr.in
# 将文件中第34行的内容替换为变量fras存储的内容
awk '{if(NF==4 && ($2 > 400 || $3 > 100 )) print FILENAME, $0}' $filname | awk 'END{print NR}'
awk 'NR==4 {print FILENAME, $2}' file1.txt | tee -a file2.txt
# 查看文本文件第一列的去重信息
awk '!seen[$1]++' hg38.all.cpgSite.bed # 输出第一次遇到的第一列,而忽略后续的相同项
awk '{ print $1 }' hg38.all.cpgSite.bed | uniq # 去重,若后面加-c,则会显示每一个元素的重复次数, 但重复行不相邻时,该命令不起作用
awk '{ print $1 }' hg38.all.cpgSite.bed | sort -u # 对第一列去重并排序
# 将文本文件按第一列特定字符串及第二列排序
sort -k 1,1V -k 2,2n ${file} | awk '{ if($1=="chr1"||$1=="chr2"||$1=="chr3"||$1=="chr4"||$1=="chr5"||$1=="chr6"||$1=="chr7"||$1=="chr8"||$1=="chr9"||$1=="chr10"||$1=="chr11"||$1=="chr12"||$1=="chr13"||$1=="chr14"||$1=="chr15"||$1=="chr16"||$1=="chr17"||$1=="chr18"||$1=="chr19"||$1=="chr20"||$1=="chr21"||$1=="chr22"||$1=="chrX"||$1=="chrY") { print } }' > ${file}.sort
# -k 1,1V: 指定按照第一列进行排序,-k1,1 表示只考虑第一列,V 表示按照版本号进行排序(版本排序,例如,10会排在2的前面而不是后面)。
# -k 2,2n: 如果第一列相同,则按照第二列进行数值(numeric)排序,-k2,2 表示只考虑第二列,n 表示按照数值进行排序。
# 详细解释可以参照该网站 https://www.cnblogs.com/51linux/archive/2012/05/23/2515299.html
ERA5提供逐小时 0.25分辨率,垂直37层的再分析资料。也提供月平均资料。但是没有日平均资料,因此如果需要日平均资料,则需要自己将24小时数据下载后再处理。关于分辨率,除了0.25°的,在下载时,好像也可以通过更改 python下载脚本中的 'grid': [1.0, 1.0],
参数进行下载。
The ERA5 dataset contains one (hourly, 31 km) high resolution realisation (referred to as “reanalysis” or “HRES”) and a reduced resolution ten member ensemble (referred to as “ensemble” or “EDA”).
官网文档介绍:https://confluence.ecmwf.int/display/CKB/ERA5%3A+data+documentation
ERA5数据中Analysis and forecast的区别:https://confluence.ecmwf.int/pages/viewpage.action?pageId=85402030
数据下载的准备工作(注册,获得API KEY,在本地创建 .cdsapirc 文件存储API KEY):https://cds.climate.copernicus.eu/api-how-to
ERA5的数据下载使用的python软件包是 cdsapi,和之前 下载ERA-Interim 所用的 ecmwfapi 不同。
数据下载链接 :
逐小时:https://cds.climate.copernicus.eu/cdsapp#!/dataset/reanalysis-era5-pressure-levels?tab=overview
逐日数据下载接口(有多个分辨率可选,但好像只能逐月的下载,且没有提供python下载脚本) https://cds.climate.copernicus.eu/cdsapp#!/software/app-c3s-daily-era5-statistics?tab=app
从上述网址获得python下载代码中没有指明 expver(experiment version),因此下载得到的数据中会多一个expver的维度。其中,expver=0001是ERA5,expver=0005是ERA5T,因此可以再加一行选项'expver': '1'
来进行筛选。
其中,ERA5T是初始发布数据,即滞后于实时时间不超过三个月的数据。如果在ERA5T中检测到严重的缺陷,该数据可能与最终的ERA5数据不同。但在实践中,这种情况不太可能发生。根据迄今为止ERA5的生产经验(以及过去的ERA-Interim),ECMWF预计这样的事件不会每几年发生一次以上,如果真的发生的话。在不太可能的情况下,需要这样的纠正,将尽快通知用户。
关于netcdf文件中ERA5T和ERA5的具体说明可以参考https://confluence.ecmwf.int/display/CUSF/ERA5+CDS+requests+which+return+a+mixture+of+ERA5+and+ERA5T+data
#!/usr/bin/env python
import cdsapi
varname = ['u_component_of_wind','v_component_of_wind']
filname = ['u','v']
c = cdsapi.Client()
for nv in range(0,len(varname),1):
for year in range(1979,2020,1):
c.retrieve(
'reanalysis-era5-pressure-levels', # 也可以是 'reanalysis-era5-single-levels-monthly-means',
{
'product_type': 'reanalysis', # 也可以是 'monthly_averaged_reanalysis'
'expver': '1',
'variable': varname[nv],
'pressure_level': [
'200', '225', '250',
'450', '500', '825',
'850',
],
'year': str(year),
'month': [
'01', '02', '03',
'04', '05', '06',
'07', '08', '09',
'10', '11', '12',
],
'day': [
'01', '02', '03',
'04', '05', '06',
'07', '08', '09',
'10', '11', '12',
'13', '14', '15',
'16', '17', '18',
'19', '20', '21',
'22', '23', '24',
'25', '26', '27',
'28', '29', '30',
'31',
],
'time': [
'00:00', '06:00', '12:00',
'18:00',
],
'area': [
90, -180, 0,
180,
],
'grid': [1.0, 1.0
],
'format': 'netcdf',
},
'/home/users/data/ERA5_subdaily_NH_%s_%d.nc'%(varname[nv],year))
下面这个网址可以将ERA5的逐月数据交互式可视化(如下图),并且提供了相应的python代码: https://cds.climate.copernicus.eu/apps/c3s/app-era5-explorer
ERA5的降水数据是逐小时的,如果要计算日累积降水,就需要将24小时的降水数据加起来,这里提供了代码: https://confluence.ecmwf.int/display/CKB/ERA5%3A+How+to+calculate+daily+total+precipitation