Python有五个标准的数据类型:Numbers(数字),String(字符串),List(列表),Tuple(元组),Dictionary(字典)
# 列表是写在方括号 [] 之间、用逗号分隔开的元素列表
list = [ 'abcd', 786 , 2.23, [1,2,3]]
tinylist = [123, 'runoob']
list2 = [i for i in range(30) if i % 3 == 0]
# 得到列表 [0, 3, 6, 9, 12, 15, 18, 21, 24, 27]
names = ['Bob','Tom','alice','Jerry','Wendy','Smith']
new_names = [name.upper()for name in names if len(name)>3]
# 得到列表 ['ALICE', 'JERRY', 'WENDY', 'SMITH']
print (len(list)) # 列表元素个数
print (list[3][0]) # 嵌套列表索引
print (list[1:3]) # 从第二个开始输出到第三个元素
print (list[2:]) # 输出从第三个元素开始的所有元素
print (list[-1::-1]) # 翻转列表
print (tinylist * 2) # 输出两次列表
print (list + tinylist) # 连接列表
元组与列表类似,但元组的元素不能修改。
元组写在小括号 () 里,元素之间用逗号隔开
元组的索引及切片方式同list
# 字典用 {} 标识,它是一个无序的 键(key) : 值(value) 的集合
# 同一个字典中,键(key)必须是唯一且不可变的,因此数字,字符串或元组都可作为键,但列表不行。
dict1 = {}
dict1['one'] = "1 - 菜鸟教程"
dict1[2] = "2 - 菜鸟工具"
tinydict = {'name': 'runoob','code':1, 'site': 'www.runoob.com'}
dict2 = {x: x**2 for x in (2, 4, 6)}
# 得到的字典为{2: 4, 4: 16, 6: 36}
print (len(dict2)) # 计算字典元素个数,即键的总数。
print (dict['one']) # 输出键为 'one' 的值
print (tinydict.keys()) # 输出所有键
print (tinydict.values()) # 输出所有值
del tinydict['Name'] # 删除键 'Name'
tinydict.clear() # 清空字典
python列表和元组顺序也是从0开始的
a=input("please input password")
print("password you input is ",password)
print(type(a)) # input variable is usually str
b=int(a)
print(type(b))
a=2
b='casename'
print('%s month is %03d'%(b,a))
# 格式化输出结果:casename month is 002
python中除号用/表示,但是和C语言不同的是/得到的值总是浮点数,例如:5 / 5结果是1.0。 python中整除用//表示是,//表示两数相除,向下取整,例如8 // 5 结果是1,,至于数据类型是整数还是浮点数则取决于用到的数据类型。
不管是条件语句还是循环,命令最后都有冒号。
if a == 0 and b != 14:
b = b+a
elif not a <= -1:
b = b-a
else:
b = b*a
list1 = ['a','b','coin']
if 'a' in list1:
print('you are right')
for i in range(-1,-10,-4):
print("Number of cycles: %d"%i)
# 输出结果:
# Number of cycles: -1
# Number of cycles: -5
# Number of cycles: -9
name = 'chengdu'
for i in name:
print(i,end='\t')
# '\t'是转义字符,表示水平制表符,即 Tab 键,一般相当于四个空格
# 输出结果:c h e n g d u
b = [np.arange(d-lagh,d+lagh+1) for d in range(2,10,3)]
# 列表生成式
coins = ["Bitcoin", "Ethereum", "Cardano"]
prices = [48000,2585,2]
for coin, price in zip(coins, prices):
print(f"${price} for 1 {coin}")
for i, coin in enumerate(coins):
price = prices[i]
print(f"${price} for 1 {coin}")
i = 0
while i < 5 :
print('cycle time now : %d, i = %d'%(i+1,i))
i += 1
else :
print('cycle end, now i = %d'%i)
break语句跳出 for 和 while 的循环体
continue 跳过当前循环,直接进行下一轮循环
pass 是空语句,一般用作占位语句,不做任何事情
除了正常定义的位置参数(也就是必选参数)外,还可以使用默认参数、可变参数(传入的参数个数可变)和关键字参数,使得函数定义出来的接口,不但能处理复杂的参数,还可以简化调用者的代码。
https://www.liaoxuefeng.com/wiki/1016959663602400/1017261630425888
import netCDF4 as nc
f=nc.Dataset("/home/users/qd201969/ERA5-1HR/stat_apr-sep1979-2019.nc") # best to use absolute path
f.variables.keys() # get the list of the variables
a=f.variables['mstr'] # get the data and vari info
a=f.variables['mstr'][:] # only get the data
lat=f.variables['latitude'][:]
lon=f.variables['longitude'][:]
indx=np.argwhere((lon<100) & (lon>10))
indy=np.argwhere((lat<50) & (lat>10))
从xarray走向netCDF处理: https://www.jianshu.com/p/86f02bc58265
import xarray as xr
f = xr.open_dataset('EC-Interim_monthly_2018.nc')
a=f['mstr'] # get the data and vari info
a=f.mstr # two method to read the variable
b=a.loc[84.86:90,352.5:360] # use lat and lon read data
c=b.mean(dim=['lon','lat'])
CDO, climate data operator的缩写。提供了600多个常见的操作,能快速处理nc、grid等常见数据.
常见的功能包括:
参考:https://blog.sciencenet.cn/blog-1081898-1275862.html
转换数据类型或者转换文件类型
cdo -b F64 copy input.nc output.nc
# 将input.nc文件中的变量转换为 double,并另存为output.nc。如果不加F,则是将所有变量转换为double。加F则只是将floating data转为double
cdo -f nc copy input.grib output.nc
# 将grib文件转换为nc文件
cdo -r -copy ff_stat_[1-9].nc ff_stat_1[0-2].nc outfilename
# 每月数据放在一个文件中,将1-12月共12个文件合并,-r表示合并后添加时间维,拼接后原始文件依然存在
cdo -cat uwind.1985-07.daily.nc uwind.1985-08.daily.nc outfilename
# 将两个文件中的变量按原有的时间维拼接,拼接后,原始文件依然存在
# 当输入文件过大且变量为short变量时,容易产生报错,解决方法是改用以下的代码,将short变量转换为float,-mergetime的作用等同于-cat
cdo -b F32 -mergetime ERA5_precip_1980-[1-9].nc ERA5_precip_1980-1[0-2].nc ERA5_precip_1hr_dec-jan1980.nc
cdo selmonth,1 ERA5_NH_z_1981.nc ERA5_NH_z_1981-01.nc
# 提取一个月的数据,注意选项与参数间的逗号
cdo sellevel,850 ERA5_NH_z_1981.nc ERA5_NH_850z_1981.nc
# 提取某个高度的数据,注意选项与参数间的逗号
cdo expr,’speed=sqrt(sqr(uwnd)+sqr(vwnd));var2=ts-273.15;’ infile outfile
# infile中有变量uwnd,vwnd,ts。由这三个变量计算新变量并存储入 outfile
除了处理数据外,cdo也能快速查看 grid、nc 等数据文件信息,常用参数有:
cdo sinfon ERA5_wind10_2016.grib #输出数据文件的简短信息
参考自:https://cloud.tencent.com/developer/article/1618318
dimension:维度
coordinate:坐标
一个维度可以被赋予多个坐标,因此一般坐标数量多于纬度数
data = np.random.rand(4, 3)
locs = ["IA", "IL", "IN"]
times = pd.date_range("2000-01-01", periods=4)
foo = xr.DataArray(data, coords=[times, locs], dims=["time", "space"])
foo = xr.DataArray(data, coords=[("time", times), ("space", locs)])
foo.attrs['units']='mm/day' # 修改或创建数据属性参数
del foo.attrs['projection'] # 删去某个属性
foo.drop('time',dim=None) # 删去time维对应的坐标,但保留维度
foo.rename({'dayofyear':'time'}) # 修改某个纬度的名字
foo.coords['time'] = times # 修改某个维度的坐标
foo.dtype # dtype('float32')
foo.dims # ('time', 'space')
var.shape # (10248, 1, 256, 512)
var.min()
temp = 15 + 8 * np.random.randn(2, 2, 3)
precip = 10 * np.random.rand(2, 2, 3)
lon = [[-99.83, -99.32], [-99.79, -99.23]]
lat = [[42.25, 42.21], [42.63, 42.59]]
ds = xr.Dataset(
{
"temperature": (["x", "y", "time"], temp),
"precipitation": (["x", "y", "time"], precip),
},
coords={
"lon": (["x", "y"], lon),
"lat": (["x", "y"], lat),
"time": pd.date_range("2014-09-06", periods=3),
"reference_time": pd.Timestamp("2014-09-05"),
},
)
# Dictionary like methods to update a dataset
ds = xr.Dataset()
ds["temperature"] = (("x", "y", "time"), temp)
ds["temperature_double"] = (("x", "y", "time"), temp * 2)
ds["precipitation"] = (("x", "y", "time"), precip)
ds.coords["lat"] = (("x", "y"), lat)
ds.coords["lon"] = (("x", "y"), lon)
ds.coords["time"] = pd.date_range("2014-09-06", periods=3)
ds.coords["reference_time"] = pd.Timestamp("2014-09-05")
ds = da.to_dataset(name="temp") #将da这个DataArray加入到ds这个数据组中,并命名为 temp
ds.to_netcdf(fileout+varname+".nc","w") #存储ds文件
# the use of pipe
plt.plot((2 * ds.temperature.sel(x=0)).mean("y"))
(ds.temperature.sel(x=0).pipe(lambda x: 2 * x).mean("y").pipe(plt.plot))
To remove a dimension or a variable, you can use below method, which will return a new dataset. Any variables using that dimension are dropped:
ds_new = ds.drop_vars("temperature")
ds_new = ds.drop_dims("time")
ds_new.coords['month']=range(1,13,1) # add a new dimension
import xarray as xr
import numpy as np
lonl=0
lonr=150
lats=15
latn=70
ds = xr.open_dataset("ERA5_precip_1hr_dec-jan1989.nc")
lat = ds.latitude
lon = ds.longitude
ilon = lon[(lon>=lonl) & (lon<=lonr)]
ilat = lat[(lat>=lats) & (lat<=latn)]
term = ds['tp'][0,:,:]
term = ds['tp'].sel(time=ds.time.dt.year.isin(1989),longitude=ilon,latitude=ilat)
# <xarray.DataArray 'tp' (time: 8760, latitude: 221, longitude: 601)>
term.mean(('latitude','latitude'))
term.sel(time=slice("2000-01-01", "2000-01-02"))
var = term.sel(time=term.time.dt.month.isin(1)).mean("time")
term.loc[ct,:,:] = term.sel(time=ct).mean('time')
# loc是定位,sel是截取部分变量
# calculate monthly data
var = term.groupby(term.time.dt.month)
# DataArrayGroupBy, grouped over 'month'
# 12 groups with labels 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12.
print(var.sum("time"))
# <xarray.DataArray 'tp' (month: 12, latitude: 221, longitude: 601)>
ds1=ds.groupby(ds.time.dt.month).mean('time')
ds1=ds.groupby('time.month').mean('time')
ds1=da.groupby('time.dayofyear').mean('time') # 365天
ds_weighted = (ds * weights).groupby('time.season').sum(dim='time') # 'DJF' 'JJA' 'MAM' 'SON'
fs = xr.open_mfdataset(file_string, combine='nested', concat_dim='time', parallel=True)
# file_string 可以是一个字符串“path/to/my/files/*.nc”,也可以是文件路径列表
# combine 有两种选项 "by_coords" 和 "nested"
# concat_dim 表示将文件沿某一维度连接,只有 combine='nested' 时需要提供这一选项。
# 若concat_dim是已存在的维度,则文件沿该维度连接,若concat_dim是不存在的维度,则文件沿新的维度连接
var = fs['temp'].data
term=[var[i:i+3,:,:].sum(axis=0) for i in range(0,25587,3)] # 得到一个list
var1 = np.array(term).reshape(8529*nlat*nlon) #转成一维数组
主要依赖的库是 https://pypi.org/project/cfgrib/ ,该库可以把grib文件投影到netcdf格式。
可以搭配xarray使用,使用方法如下。可以一次读取一个文件,也可以用xr.open_mfdataset读取多个文件,其他用法同netcdf
ds = xr.open_dataset('example.grib', engine='cfgrib')
print(ds)
# 目前cfgrib库无法同时读取多个typeOfLevel,因此需要根据提示筛选我们需要的数据。
# 这个typeOfLevel可以是isobaricInhPa,surface,depth_below_land
ds = xr.open_dataset('example.grib2', engine="cfgrib", backend_kwargs={'filter_by_keys': {'typeOfLevel': 'isobaricInhPa'}})
print(ds)
默认情景下,cfgrib会输出一个index文档(后缀为.idx),该索引文档可随时删除它们,并在出现问题时重试。 当grib文件所处目录不支持写入文件时,即cfgrib无法在该目录下保存索引文档,此时便会出现如下报错信息,但该报错信息不影响grib数据的读取。 也可在backend_kwargs中加入’indexpath’:’‘选项来选择不输出该索引文档,此时上述报错信息便会消失。
Can't create file '/home/metctm1/array_hq133/data/drv_fld/era5/20220803-sl.grib.923a8.idx'
Traceback (most recent call last):
File "/home/metctm1/array/soft/anaconda3/lib/python3.8/site-packages/cfgrib/messages.py", line 343, in from_indexpath_or_filestream
with compat_create_exclusive(indexpath) as new_index_file:
File "/home/metctm1/array/soft/anaconda3/lib/python3.8/contextlib.py", line 113, in __enter__
return next(self.gen)
File "/home/metctm1/array/soft/anaconda3/lib/python3.8/site-packages/cfgrib/messages.py", line 264, in compat_create_exclusive
fd = os.open(path, os.O_WRONLY | os.O_CREAT | os.O_EXCL)
PermissionError: [Errno 13] Permission denied: '/home/metctm1/array_hq133/data/drv_fld/era5/20220803-sl.grib.923a8.idx'
Can't read index file '/home/metctm1/array_hq133/data/drv_fld/era5/20220803-sl.grib.923a8.idx'
Traceback (most recent call last):
File "/home/metctm1/array/soft/anaconda3/lib/python3.8/site-packages/cfgrib/messages.py", line 353, in from_indexpath_or_filestream
index_mtime = os.path.getmtime(indexpath)
File "/home/metctm1/array/soft/anaconda3/lib/python3.8/genericpath.py", line 55, in getmtime
return os.stat(filename).st_mtime
FileNotFoundError: [Errno 2] No such file or directory: '/home/metctm1/array_hq133/data/drv_fld/era5/20220803-sl.grib.923a8.idx'
若要循环处理多个文件,可以先利用list来存储数据,然后再转成numpy进行数学运算
var = []
fn_stream = subprocess.check_output(
'ls %s/%d07*-%s.grib'%(
indir[0],year,suffix), shell=True).decode('utf-8')
fn_list = fn_stream.split()
for itm in fn_list:
ds = xr.open_dataset(itm,engine='cfgrib',backend_kwargs={
'filter_by_keys': {'typeOfLevel': 'surface'}})
var.append(ds[varname].data.tolist())
var = np.array(var).mean(axis=0)
da = xr.DataArray(np.arange(16).reshape(4, 4), dims=["x", "y"])
da.where(da.x + da.y < 4)
# 保留满足条件的变量,未满足条件的变量设为nan
# da.where(da.x + da.y < 4,0),会将未满足条件的变量设为0
# output is
# <xarray.DataArray (x: 4, y: 4)>
# array([[ 0., 1., 2., 3.],
# [ 4., 5., 6., nan],
# [ 8., 9., nan, nan],
# [12., nan, nan, nan]])
# Dimensions without coordinates: x, y
da = xr.DataArray([1, 2, 3, 4, 5], dims=["x"])
lookup = xr.DataArray([-1, -2, -3, -4, -5], dims=["x"])
da.where(lookup.isin([-2, -4]), drop=True) # output array([2., 4.])
# 加 drop=True后,未满足条件的变量会被删去,数据结构发生变化
# when done repeatedly, this type of indexing is significantly slower than using sel().
from scipy import stats
def new_linregress(x, y):
# Wrapper around scipy linregress to use in apply_ufunc
slope, intercept, r_value, p_value, std_err = stats.linregress(x, y)
return np.array([slope, p_value])
a1 = xr.apply_ufunc(new_linregress, t, t,
input_core_dims=[['time'], ['time']],
output_core_dims=[["parameter"]],
vectorize=True,
dask="parallelized",
output_dtypes=['float64'],
)
# dimension of t is (time: 10248, lat: 3, lon: 4)
# a1 is (lat: 3, lon: 4, parameter: 2)
ts=np.zeros([2,3,10248],dtype=float)
locs = ['EA','NA']
month = ['DJF','MAM','JJA']
ts1 = xr.DataArray(ts, coords=[locs,month,t.time], dims=["space","month","time"])
a1 = xr.apply_ufunc(new_linregress, t, ts1,
input_core_dims=[['time'], ['time']],
output_core_dims=[["parameter"]],
vectorize=True,
dask="parallelized",
output_dtypes=['float64'],
)
# a1 is (lat: 3, lon: 4, space: 2, month: 3, parameter: 2)
绝热与非绝热经常搞混的几个单词: diabatic 非绝热的,传热的 adiabatic 绝热的,不传热的 adiabatically
The structure and dynamics of decadal anomalies in the wintertime midlatitude North Pacifc ocean–atmosphere system are examined in this study, using the NCEP/NCAR atmospheric reanalysis, HadISST SST and Simple Ocean Data Assimilation data for 1960–2010.
瞬变涡旋对气候态的影响在大气环流的理论和观测研究中都是一个基本问题。TE可以通过重新分配热量和涡度影响大气状态。因此TE对气候态的平衡可以用TE通量的辐合辐散表示。
但是,静力平衡和准地转约束要求TE通量的产生必须伴随次级环流,而这又会影响气候态。因此若要全面理解TE对平均流的净效应,除了要考虑TE通量的辐合辐散外,也要考虑由TE引起的环流的影响。
以往的研究通过假设TE通量的辐合辐散为热量和涡度的源汇,通过求解方程,计算出定常响应,来解决上述问题。这些研究表明TE会减弱由地形强迫出来的定常波的振幅。但这些耗散效应的定量估计对模式假设非常敏感。
除了寻找时间平均方程的有限振幅解外,也可以通过考察当加入TE强迫时的瞬时初始大气响应。且最好用准地转系统来研究这个问题,这可以保证计算结果满足静力平衡及地转平衡。如果用原始方程来研究这个问题,会获得一个有很多不真实小扰动的初始倾向。
因此Lau于1984年利用时间平均的准地转位涡方程,通过设定合适的边界条件(其中热力强迫项的边界条件根据静力平衡和时间平均温度公式推导得到),求解出由TE驱动出来的三维分布的准地转位势倾向,然后根据地转风关系和静力平衡求解地转风倾向及温度倾向。这些初始倾向满足静力平衡及地转平衡,因此可以将由eddy引起的环流和TE通量散度同时考虑进去。
由于Lau1984年的文章只考虑了准地转大气对TE强迫的响应,因此可以将平流项和摩擦项等忽略。
反算泊松算子求解位势倾向时,Lau1984用的是球谐波函数差分迭代(我也不太懂),而现在多用超松弛迭代法。
缺点:
\begin{split}
\bar \sigma_0 = \frac{R}{C_p} \frac{\bar T}{p} - \frac{\partial \bar T}{\partial p} \\
\bar \sigma_1 = \frac{R \bar \sigma_0}{p}
\end{split}
迭代计算位势倾向时,由于热力项的垂直边界条件和瞬变涡旋强迫的不同,导致其计算结果也不同。理想情况下的解析解情况:
一般在风暴轴处伴随有向极的eddy热量通量和eddy涡度通量。而这向极的eddy热量通量,会增强上层东风,增强下层西风。向极的eddy涡度通量则会增强整层的西风。因此eddy热量通量和eddy涡度通量对纬向风的作用在底层一样(增强西风),高空则相反。
翻译自:https://www.ecmwf.int/en/newsletter/159/meteorology/global-reanalysis-goodbye-era-interim-hello-era5
2020年第一季度时,ERA5被发布,以替代自2006年开始发布的ERA-Interim。ERA5基于4D-Var data assimilation using Cycle 41r2 of the Integrated Forecasting System (IFS),该系统于2016年开始在ECMWF运行。
相比于ERA-Interim, ERA5的优势在于模式物理、动力核、数据同化的改进,且分辨率从ERA-Interim的79km提高到31km。此外,ERA5可以提供逐小时输出及一个不确定性估计。The uncertainty information is obtained from a 10-member ensemble of data assimilations with 3-hourly output at half the horizontal resolution (63 km grid spacing).且ERA5输出的气象变量更丰富,例如多了100m风。
现在已经到了逐步淘汰ERA-Interim的时候了,因为ERA-Interim使用几种新的重要观测数据的能力有限,越来越难以维持,也不会迁移到欧洲中期天气预报中心在Bologna的未来数据中心。
在ERA5中使用的观测数据从1979年0.75MB/day增长到2018年24MB/day。其中,卫星辐射是最主要的增长数据。2009年后的降水数据中还同化了雷达观测。
相比于ERA-Interim,ERA5利用all-sky方法同化对湿度敏感的卫星频道。这可以在多云和降水地区提供新的信息,也纠正了早期在多雨条件下辐射同化技术导致20世纪90年代全球海洋ERA-Interim异常降水的一个问题。
传统数据和卫星数据的特征描述、相互校准和处理的改进逐步改进了历史观测的质量,并扩大其地理和时间覆盖范围。ERA5利用了几个再处理的卫星数据集,这些数据来自欧洲、美国和日本的空间机构和研究所。其中包括大气运动风矢量;臭氧、无线电掩星和测高数据;来自散射计的土壤湿度和风数据;以及对海洋湿度敏感的卫星数据的SSMI记录。
ERA-Interim无法吸收来自最新卫星仪器的观测数据(例如IASI和CrIS的高光谱数据)和地表雷达数据。此外,由于某些仪器和通道的失灵,ERA-Interim使用的观测数据逐渐减少。ERA5还可以处理自2013年以来越来越多地被地面压力、高空风和温度数据所采用的BUFR格式。而ERA-Interim无法处理这种形式,因此,它所吸收的观测数量急剧减少。在2018年底,ERA5每天使用约2400万次观测,大约是ERA-Interim的5倍。总得来说,ERA5使用的观测数据比ERA-Interim多很多。
ERA5使用的模型输入强迫包括太阳总辐照度强迫、臭氧、温室气体和为世界气候研究计划(WCRP) CMIP5倡议开发的一些气溶胶,包括平流层硫酸盐气溶胶。这是对ERA-Interim的一个重大改进,例如,ERA-Interim没有考虑到主要火山爆发造成的平流层硫酸盐气溶胶。
海表温度(SST)和海冰覆盖的演变基于不同时期的若干产品:英国气象局哈德利中心的海温和海冰HadISST2产品,EUMETSAT OSI-SAF再分析海冰数据,以及英国气象局的海温和海冰OSTIA产品,该产品也用于ECMWF的业务预报系统。细节可以在Hirahara et al.(2016)中找到。使用这些数据的目的是产生一个合并的数据集:
再分析是将模式信息与观测信息混合在一起的数据同化系统的结果,可以提供随时覆盖全球的数据。再分析数据总体上比前卫星时代更准确,因为当时的观测相对稀疏。ERA5使用的10个集合的4D-Var系统提供了对数据不确定性的估计。这些估计取决于数据同化中使用的短期预测中与流体相关的不确定性。它们还很大程度上取决于观测覆盖范围。
从ERA5再分析开始的再预测显示,相比于ERA-Interim,在技术范围内增加了大约一天的提前预报天数。此时,可以更好地模拟热带气旋地中心气压和降水量。逐小时地时间分辨率使我们能够更精确地观察每日天气系统的演变。