python 利用Xarray处理nc数据
注意数据的格式(经纬度范围、分辨率、时间格式、变量等基本信息）
import xarray as xr
f = xr.open_dataset('D:\download\hgt.mon.mean.nc')
print(f)
#查看维度信息
lat=f.lat
lon=f.lon
time=f.time
print(lat)
时间格式为：***datetime64[ns]***时最方便，直接索引即可。
eg：
* time (time) datetime64[ns] 1948-01-01 1948-02-01 ... 2019-01-01
需要选取2005年到2010年冬季的数据
#选取2005年到2010年冬季的数据
hgt = f['hgt'].loc[f.time.dt.month.isin([12,1,2])].loc ['2005-12-01':'2011-02-01']
print(hgt)
#选取夏季数据
hgt_JJA= f['hgt'].loc[f.time.dt.season=='DJF']
print(hgt_JJA)
#叠加时间和范围的选取
#a = tmax.loc[:,40:55,115:135].loc[tmax.time.dt.season=='JJA']
当格式为object
import xarray as xr
f = xr.open_dataset('/data/home/zenggang/yxy/CMIP6/SSP245/CanESM5/tas/t2m.nc')
f= f.assign_coords(time = f.indexes['time'].to_datetimeindex())
tas = f['tas'].loc[f.time.dt.month.isin([12,1,2])].loc ['2015-12-01':'2100-02-28']
有的模式的时间是0时(比如说 1979-01-01 00:00:00)，而有的模式的时间是12时(比如说 1979-01-01 12:00:00)，如果需要统一的话，只需要搭配datetime库进行调整即可，比如说将所有的12时变为0时：
from datetime import datetime, timedelta
f = f.assign_coords(time = (f.indexes['time'].to_datetimeindex()- timedelta(hours=12)))
#在得到了一个DataArray之后，用于画图时，需要获取他的经度、纬度、高度等
f.coords['time']
f['time']
print(f['time'])
Xarray官方提供了三种方法用来索引数据：
1.利用下标索引(index)
2.利用坐标值索引(coords)
3.利用标签索引(labels)
da指DataArray;ds指Dataset
#method 1
print(hgt[:2])
#method 2
print(hgt[:, [2, 1]])
#method 3
print(f['hgt'].loc['2000-01':'2003-01'].lat[45:90])
h=f['hgt'] .loc['time1':'time2',lat1:lat2,lon1:lon2]
print(h)
#索引出来的数据直接可以修改或者计算
#数据缺测处理
print(hgt.where(hgt.lat < 60))
首先建立一个带有缺测的dataarray
x = xr.DataArray([0, 1, np.nan, np.nan, 2], dims=['x'])
print(x)
判断是否缺测,还有notnull()函数
print(x.isnull())
判断有几个未缺测
print(x.count())
去除缺测
print(x.dropna(dim='x'))
将缺测填为-1
print(x.fillna(-1))
.interp()函数
插值
比如说气象上需要求逐周，逐旬，逐候的数据，处理起来很麻烦，我们则可以使用这个方法先将其聚块，再进行平均。
print(da.coarsen(time=7).mean())
coarsen(time=7)将时间轴每7天固定在一起，通过平均函数，得到了逐周的平均结果。
但是这里364是可以整除7的，如果不能刚好整除的话，还提供了一些参数来处理边界结果。
boundary='trim'，或者boundary='pad'.
trim是指末尾多余的不参与计算，或者叫修建多余条目。
pad指填充nan到不足的地方，补全为可以整除的更大数据再进行聚块。
还可以通过coord_func函数返回不同的坐标信息。
比如说默认下，7天平均，结果的时间是第4天，就是中间那天，但是通过coord_func={'time': 'min'}可以指定时间为最小的那天。
split·将数据分为多个独立的组。
apply·对各个组进行操作。
combine·将各个组合并为一个数据对象。
读取nc文件
ds = xr.open_dataset('ds.nc')
xr.open_mfdataset('my/files/*.nc', parallel=True)
生成nc文件
ds = xr.Dataset({'prec': (('xy', 'time'), np.random.rand(4, 5))},coords={'lat': ('xy', [15, 25, 35, 45]),'lon': ('xy', [15, 25, 35, 45]),'time': pd.date_range('2000-01-01', periods=5),})
print(ds)
ds.to_netcdf('output.nc')
data = pd.read_csv("cslist.txt",sep=',',header=None, names=['a','b','c','d','e','f','g','h','i','j','k'])
读取txt的同时，对每列赋予了一个列名，通过data.a可以直接按列名调用相应数据。
对于较复杂的.txt文件，仍可通过该函数读取
grib文件可通过pygrib库读取
import pygrib
f = pygrib.open('xxx.grb')
常用的统计方法
import numpy as np
平均值
在计算气候态，区域平均时均要使用到求均值函数，对应NCL中的dim_average函数，在python中通常使用** np.mean()**函数
💦气象数据中常有缺测，在NCL中，使用求均值函数会自动略过，而在python中，当任意一数与缺测(np.nan)计算的结果均为np.nan，比如求[1,2,3,4，np.nan]的平均值，结果为np.nan
因此，当数据存在缺测数据时，通常使用**np.nanmean()函数，用法同上，此时[1,2,3,4，np.nan]的平均值为(1+2+3+4)/4 = 2.5
同样的，求某数组最大最小值时也有np.nanmax(), np.nanmin()**函数来补充np.max(), np.min()的不足。
其他很多np的计算函数也可以通过在前边加‘nan’来使用。
标准差
标准化
x = (x - np.mean(x)) / np.std(x)
相关系数
皮尔逊相关系数：
numpy函数中的np.corrcoef(x, y)，可以实现相关计算
scipy.stats中的函数
from scipy.stats import pearsonr
r,p = pearsonr(x, y)
优点是可以直接返回相关系数R及其P值，这避免了我们进一步计算置信度。而缺点则是该函数只支持两个一维数组的计算，也就是说当我们需要计算一个场和一个序列的相关时，我们需要循环来实现。
r = np.zeros((a.shape[1],a.shape[2]))
p = np.zeros((a.shape[1],a.shape[2]))
for i in range(sic.shape[1]):
for j in range(sic.shape[2]):
r[i,j], p[i,j] = pearsonr(b , a[:,i,j])
其中a[time,lat,lon]，b[time]
线性回归系数
(NCL中为regcoef()函数)
同样推荐Scipy库中的stats.linregress(x,y)函数：
slope, intercept, r_value, p_value, std_err = stats.linregress(x,y)
slop: 回归斜率
intercept：回归截距
r_value： 相关系数
p_value： P值
std_err： 估计标准误差
直接可以输出P值，同样省去了做置信度检验的过程，遗憾的是仍需同相关系数一样循环计算。
💫来源：简书-摸鱼咯大佬！！！！✍️