学习笔记3 | 高维数据处理——Xarray

目录

一、数据结构

1.DataArray

（1）DataArray的创建

（2）DataArray的属性及常用方法

2.DataSet

（1）DataSet的创建

（2）DataSet的属性和常用方法

二、数据的读取

1.读取nc文件

2.读取grib文件

3.读取多个文件并 合并

三、数据的索引

1.通过位置索引

2.通过名字索引

四、数据的坐标系统

1.修改坐标

2.增加/删除坐标

五、数据统计与计算

1.基础运算

2.降维

3.分组与聚合

4.滑窗

5.插值

六、高维数据可视化

1.时间序列可视化

2.空间数据可视化

3.多子图

七、数据的掩膜

1.数值掩膜

2.空间位置掩膜

---

一、数据结构

1.DataArray

（1）DataArray的创建

1.  
   import xarray as xr
2.  
   import numpy as np
3.  
   import pandas as pd
4.  
    
5.  
   data = np.random.randn(4,3,2) #随机生成三维数组
6.  
   xr.DataArray(data) #简单的dataarray
7.  
    
8.  
   times = pd.date_range('2000-01-01',periods=4)
9.  
   lat = [0,1,2]
10.  
    lon = [10,20]
11.  
    da = xr.DataArray(
12.  
    data,
13.  
    coords={
14.  
    'time'=times,
15.  
    'lat'=lat,
16.  
    'lon'=lon
17.  
    },
18.  
    dims=['time','lat','lon']
19.  
    )

（2）DataArray的属性及常用方法

1.  
   da.dims #获取维度
2.  
   da.coords #获取坐标
3.  
   da.values #获取数值
4.  
   type(da.values)
5.  
   da.attrs #获取属性
6.  
   da.attrs['name']='precitation' #设置属性
7.  
   da.attrs['units']='mm'
8.  
   da.rename('rain')

2.DataSet

可以理解为多个DataArray的集合

（1）DataSet的创建

1.  
   #1.直接创建
2.  
   temp = 15 8*np.random.randn(4,3,2)
3.  
   precip = 10*np.random.randn(4,3,2)
4.  
   xr.Dataset(
5.  
   {
6.  
   'temperature'=(['time','lat','lon'],temp)
7.  
   'precipitation'=(['time','lat','lon'],precip)
8.  
   },
9.  
   coords={
10.  
    'lon'=lon
11.  
    'lat'=lat
12.  
    'time'=times
13.  
    }
14.  
    )
15.  
     
16.  
    #2.使用DataArray创建
17.  
    xr.Dataset({'precip':da})
18.  
    xr.Dataset({
19.  
    'precip':da,
20.  
    'temp':da 10
21.  
    })

（2）DataSet的属性和常用方法

1.  
   ds.data_vars #显示有哪些变量
2.  
   'temp' in ds #判断变量是否在Dataset中
3.  
   ds.coords
4.  
   ds._coord_names #获取坐标名
5.  
   ds._dims #获取维度
6.  
    
7.  
   ds.drop_vars('temp') #剔除变量
8.  
   ds.drop_dims('time') #剔除维度

二、数据的读取

1. nc数据：再分析数据、气候系统模拟
2. grib数据：数值预报中用的较多

1.读取nc文件

1.  
   #xarray.open_dataset() engine=netcdf4/cfgrib
2.  
   import xarray as xr
3.  
   ds = xr.open_dataset('filename') #engine默认netcdf4
4.  
   ds.longitude
5.  
   ds.longitude.values
6.  
   ds.time
7.  
   type(ds.time.values)
8.  
   ds['latitude']
9.  
   ds['t2m']
10.  
    ds['t2m'].values
11.  
     
12.  
    #xarray.open_dataarray #只能打开含有一个变量
13.  
    da = xr.open_dataset('filename',drop_variables=['u10','v10','swh'])
14.  
    da = xr.open_dataarray('filename',drop_variables=['u10','v10','swh'])
15.  
    da = xr.open_dataarray('filename')

2.读取grib文件

1.  
   xr.open_dataset('fd.grib2',engine='cfgrib') #打开后发现缺少t2m类似其他坐标系统里的变量1
2.  
   xr.open_dataset('fd.grib2',engine='cfgrib'，
3.  
   backend_kwargs={
4.  
   'filter_by_keys':{
5.  
   'typeOfLevel':'heightAboveGround',
6.  
   'level':2
7.  
   }})
8.  
   xr.open_dataset('fd.grib2',engine='cfgrib'，
9.  
   backend_kwargs={
10.  
    'filter_by_keys':{
11.  
    'typeOfLevel':'surface',
12.  
    'level':0
13.  
    }})

3.读取多个文件并 合并

1.  
   # xarray.open_mfdataset(paths,engine=' ')
2.  
   xr.open_dataset('filename')
3.  
   xr.open_mfdataset('../*.nc') #mfdataset也可以读取一个文件

三、数据的索引

1.通过位置索引

：在指定维度位置上输入相应的数值来索引，只能在dataArray中使用

1.  
   #1.通过数字查找
2.  
   import xarray as xr
3.  
   xr.open_dataset('filename')
4.  
   t = ds['t2m']
5.  
   t.values
6.  
   t.values.shape
7.  
   t[0,0,0] #得到的是dataArray
8.  
   t[:,1,1]
9.  
    
10.  
    #2.标签查找
11.  
    t.loc[:,89.75,70.25]
12.  
    t.loc[:,89.75:80,70.25:71] #维度从大到小排列，经度从小到大排列

2.通过名字索引

：索引时输入相关的维度名称和相应的数值

1.  
   #1.数字查找
2.  
   ds.isel(longitude=1,time=0) #1是位置，没有查找的坐标全部都取
3.  
    
4.  
   #2.标签查找
5.  
   ds.sel(longitude=70.25,latitude=50,time='2020-01-01') #70.25是值
6.  
    
7.  
   #3.时空的切片索引
8.  
   #空间索引
9.  
   ds.sel(latitude=slice(80,50),longitude=slice(80,155))
10.  
    ds.sortby('latitude') #按照latitude排序
11.  
    ds.sel(latitude=slice(50,80),longitude=slice(80,155))
12.  
    #时间索引
13.  
    ds.sel(time=slice('2020-02-03','2020-05')
14.  
    ds.sel(time='2020-01-01 06')
15.  
    ds.sel(time=(ds.time.dt.day==15)) #指定15号
16.  
    ds.sel(time=(ds.time.dt.month.isin([1,3,5]))) #提取1，3，5月的数据

四、数据的坐标系统

1.修改坐标

1.  
   import xarray as xr
2.  
   ds = xr.open_dataset('filename')
3.  
    
4.  
   #1.修改坐标名字
5.  
   ds.rename({
6.  
   'longitude':'x',
7.  
   'latitude':'y'})
8.  
    
9.  
   #2.修改坐标数值
10.  
    ds.assign_coords(longitude=ds['longitude'] 180)
11.  
     
12.  
    #3.时间坐标拆分
13.  
    yearv= = ds.time.dt.year
14.  
    monthv = ds.time.dt.month
15.  
    dayv = ds.time.dt.day
16.  
    hourv = ds.time.dt.hour
17.  
    ds_time = ds.assign_coords(year=yearv,month=monthv,day=dayv,hour=hourv) #只定义进去了但没有和变量关联起来
18.  
    ds1 = ds_time.set_index(time_new=('year','month','day','hour')) #增加坐标到dimension里
19.  
    ds1.sel(time_new=(2020,1,15,[0,12]))
20.  
    ds1.drop_dims('time').unstack('time_new') #但time是和变量关联的，删除之后变量也没了，所以要修改
21.  
    #修改上面三步
22.  
    ds1 = ds_time.set_index(time=('year','month','day','hour'))
23.  
    ds1.sel(time=(2020,1,15,[0,12]))
24.  
    time = ds1.unstack('time_new')
25.  
    time.sel(day=15,hour=18)

2.增加/删除坐标

1.  
   #1.增加坐标
2.  
   t = ds['t2m'].expand_dims(height=1) #但还没有关联到变量里,还没有赋值
3.  
   t.assign_coords(height=[2])
4.  
    
5.  
   #2.删除坐标
6.  
   t.squeeze('height') #只能把维度是1的坐标压缩（应该意思是只有一个值的维度吧）

五、数据统计与计算

1.基础运算

1.  
   import xarray as xr
2.  
   #1.四则运算
3.  
   ds = xr.open_dataset('....nc')
4.  
   t = ds['t2m']
5.  
   t.values
6.  
   t-273.15
7.  
   u = ds['u10']
8.  
   v = ds['v10']
9.  
   u ** 2 v ** 2 #风速的平方
10.  
     
11.  
    #2.常用函数
12.  
    t.max()
13.  
    t.min()
14.  
    t.mean()
15.  
    t.std()
16.  
     
17.  
    import numpy as np
18.  
    np.max(t) #和上面结果一样
19.  
    np.sqrt(u**2 v**2)
20.  
    abs(u)

2.降维

遇到缺测值可以自动忽略

1.  
   #1.时间降维
2.  
   t.mean()
3.  
   t.mean(dim='time') #年平均气温
4.  
   t-t.mean(dim='time') #距平气温
5.  
   #2.空间降维
6.  
   t.mean(dim=['latitude','longitude'])

3.分组与聚合

时间分辨率有：year、month、day、hour、minute、second、microsecond、season、dayofyear、dayofweek、days_in_month

1.  
   #xarray.DataArray.groupby(group,squeeze= ,restore_coord_dims= )
2.  
   t.groupby('time.month')
3.  
   t.groupby('time.season')
4.  
    
5.  
   tmonth = t.groupby('time.month')
6.  
   tmonth.groups
7.  
   for i,data in t_month:
8.  
   print(i,data['time'])
9.  
   tmonth.mean() #时间维度变成month

4.滑窗

：是在简单平均数基础上，通过顺序逐期增减新旧数据求算滑动平均，借以消除偶然变动因素，找出事物发展趋势，并进行预测的方法。

1.  
   #xarray.DataArray.rolling(dim=None,center=False,min_periods=None)
2.  
   t_time = t.mean(['longitude','latitude'])
3.  
   r = t_time.rolling({'time':5})
4.  
   r.mean() #前五个是nan
5.  
    
6.  
   r = t_time.rolling(time=5,center=True)
7.  
   r.mean() #前2个和最后2个是nan
8.  
   r.max() #其他函数也可以用
9.  
   #空间平滑
10.  
    t_region = t.mean('time')
11.  
    r = t_region.rolling(longitude=5,latitude=5,center=True)
12.  
    r.mean()

5.插值

1.  
   #Xarray.DataArray.interp(coords= ,method= ,...)
2.  
   #算法有：多维：linear、nearest；一维：linear、nearest、zero、slinear、qadratic、cubic
3.  
    
4.  
   #1.一维插值
5.  
   t_oneday = t.sel(time='2020-01-01') #有四个时间 0，6，12，18
6.  
   t_oneday.time
7.  
   import pandas as pd
8.  
   new_time = pd.data_range(t_oneday.time.values[0],t_oneday.time.values[-1],freq='3H')
9.  
   t_oneday.interp(time=new_time) #插值成功
10.  
     
11.  
    #2.多维插值
12.  
    t_region = t_t.sel(longitude=slice(100,101),latitude=slice(41,40)) #slice左右都是闭的
13.  
    import numpy as np
14.  
    nlat = np.linspace(41,40,11)
15.  
    nlon = np.linspace(100,101,11)
16.  
    t_region.interp(longitude=nlon,latitude=nlat)

六、高维数据可视化

画图对象：dataarray，使用matplotlib库

一维数据：line（）；二维数据：pcolormesh（）；其他：hist（）

1.时间序列可视化

1.  
   import xarray as xr
2.  
   da = xr.open_dataset('....nc')
3.  
   t = da['t2m']-273.15
4.  
   t_point = t.isel(latitude=10,longitude=10) #10是位置
5.  
   t_point.plot()
6.  
    
7.  
   t_month = t_point.groupby('time.month').mean()
8.  
   t_month.plot()
9.  
    
10.  
    #多个点的时间序列
11.  
    tp = t.isel(latitude=[10,100,200],longitude=10)
12.  
    tp.plot() #结果是堆叠的hist
13.  
    tp.plot(hue='latitude') #得到不同latitude的折线图

2.空间数据可视化

1.  
   t_season = t.sel(longtitude=slice(90,120),latitude=slice(50,20)).groupby('time.season').mean()
2.  
   t_season.iseal(season = 0).plot()

3.多子图

1.  
   t_season = t.sel(longtitude=slice(90,120),latitude=slice(50,20)).groupby('time.season').mean()
2.  
   g = t_season.plot(col='season',col_wrap=2) #col_wrap是两行的意思
3.  
   g.axes
4.  
    
5.  
   tp = t.isel(latitude=[10,100,200],longitude=[10,100,200])
6.  
   tp.plot(row='latitude',col='longtitude')

七、数据的掩膜

把不要的数据隐藏起来的最简便的方法就是设为nan，计算和可视化过程都不参与

1. 数值掩膜
2. 空间位置掩膜

1.数值掩膜

1.  
   #xarray.where(cond,x,y,keep_attrs=None)
2.  
   import xarray as xr
3.  
   import numpy as np
4.  
   ds = xr.open_dataset('...nc')
5.  
   u = ds['u10']
6.  
   v = ds['v10']
7.  
   ws = np.sqrt(u**2 v**2)
8.  
   ws_region = ws.mean('time')
9.  
   ws_region.plot()
10.  
    xr.where(ws_region >5,ws_region,np.nan).plot()

2.空间位置掩膜

1.shp文件 ；2.salem掩膜工具

1.  
   #1.陆地掩膜
2.  
   import geopandas as gpd
3.  
   import salem
4.  
   land_shp = gpd.read_file('filename')
5.  
   land = ws_region.salem.roi(shape = land_shp) #region of interest
6.  
   land.plot()
7.  
    
8.  
   #2.海洋掩膜
9.  
   ocean_shp = gpd.read_file('filename')
10.  
    ocean = ws_region.salem.roi(shape = ocean_shp)
11.  
    ocean.plot()
12.  
     
13.  
    #3.行政区划掩膜
14.  
    china_shp = gpd.read_file('filename')
15.  
    china = ws_region.salem.roi(shape = china_shp)
16.  
    select = china_shp[china_shp['省'].isin(['广东省','河南省'])]

这篇好文章是转载于：学新通技术网