从xarray走向netCDF处理(二):数据索引

2019-07-24 16:47:08 浏览数 (1)

xarray专题再次开讲,错过第一部分的可以先去补个课从xarray走向netCDF处理(一):数据结构及数据读取。 今天要介绍的就是xarray的索引功能,通过索引你可以对数据进行切片,从整体中提取你所关注的区域、高度或者时间。

索引核心方法

在xarray的官方文档中给出了如下几种索引方式

索引演示

代码语言:javascript复制
对如下数据进行索引演示:名为ds的DataSet,名为temp的DataArray,数据链接在文末。
根据位置索引

位置索引是最直接也是最简单的索引方式,但是位置索引只对DataArray有效,对DataSet无效。下面用两种不同方法获取相同的值。

1. 通过数字索引

代码语言:javascript复制
>>>temp[:,1,1]
<xarray.DataArray 't2m' (time: 12)>
array([249.14844, 256.4179 , 247.45125, 254.26143, 267.49988, 273.03592,
   273.91003, 273.30096, 267.9147 , 264.65695, 255.07956, 251.31808],
   dtype=float32)
Coordinates:
 longitude  float32 0.75
 latitude   float32 89.25
* time     (time) datetime64[ns] 2018-01-01 2018-02-01 ... 2018-12-01
Attributes:
 units:     K
 long_name:  2 metre temperature

2. 通过标签索引

代码语言:javascript复制
>>>temp.loc[:, 89.25, 0.75]
<xarray.DataArray 't2m' (time: 12)>
array([249.14844, 256.4179 , 247.45125, 254.26143, 267.49988, 273.03592,
   273.91003, 273.30096, 267.9147 , 264.65695, 255.07956, 251.31808],
   dtype=float32)
Coordinates:
 longitude  float32 0.75
 latitude   float32 89.25
* time     (time) datetime64[ns] 2018-01-01 2018-02-01 ... 2018-12-01
Attributes:
 units:     K
 long_name:  2 metre temperature
根据维度名字索引

通过维度的名字就可以不必按照指定的维度顺序对数据进行切片了。对DataArrayDataSet都有效,且方法一致。

1.通过数字索引

代码语言:javascript复制
>>>ds.isel(longitude=1, time=0)
<xarray.Dataset>
Dimensions:    (latitude: 241)
Coordinates:
 longitude  float32 0.75
* latitude   (latitude) float32 90.0 89.25 88.5 87.75 ... -88.5 -89.25 -90.0
 time      datetime64[ns] 2018-01-01
Data variables:
 u10       (latitude) float32 ...
 v10       (latitude) float32 ...
 t2m       (latitude) float32 ...
Attributes:
 Conventions:  CF-1.6
 history:     2019-03-28 02:03:39 GMT by grib_to_netcdf-2.12.0: grib_to_n...

22.通过标签索引

代码语言:javascript复制
>>>ds.sel(longitude=0.75, time='2018-01-01')
<xarray.Dataset>
Dimensions:    (latitude: 241)
Coordinates:
 longitude  float32 0.75
* latitude  (latitude) float32 90.0 89.25 88.5 87.75 ... -88.5 -89.25 -90.0
 time      datetime64[ns] 2018-01-01
Data variables:
 u10       (latitude) float32 ...
 v10       (latitude) float32 ...
 t2m       (latitude) float32 ...
Attributes:
 Conventions:  CF-1.6
 history:     2019-03-28 02:03:39 GMT by grib_to_netcdf-2.12.0: grib_to_n...

索引及可视化实战

代码语言:javascript复制
import arrow
import numpy as np
import xarray as xr
import cartopy.crs as ccrs
import cartopy.feature as cfeat
from cartopy.mpl.gridliner import LONGITUDE_FORMATTER, LATITUDE_FORMATTER
import matplotlib.pyplot as plt

定义一个map函数,该函数主要用来创建地图相关的信息,后面画图的时候可以随时调用。

代码语言:javascript复制
def map():
  # 创建画图空间
  proj = ccrs.PlateCarree()  #创建坐标系
  fig = plt.figure()  #创建页面
  ax = fig.subplots(1, 1, subplot_kw={'projection': proj})  #创建子图
  # 设置地图属性
  ax.add_feature(cfeat.BORDERS.with_scale('50m'), linewidth=0.8)  #加载分辨率为50的国界线
  ax.add_feature(cfeat.COASTLINE.with_scale('50m'), linewidth=0.6)  #加载分辨率为50的海岸线
  ax.add_feature(cfeat.RIVERS.with_scale('50m'))  #加载分辨率为50的河流
  ax.add_feature(cfeat.LAKES.with_scale('50m'))  #加载分辨率为50的湖泊
  # 设置网格点属性
  gl = ax.gridlines(crs=ccrs.PlateCarree(), draw_labels=True,
                   linewidth=1.2, color='k', alpha=0.5, linestyle='--')
  gl.xlabels_top = False  #关闭顶端的经纬度标签
  gl.ylabels_right = False  #关闭右侧的经纬度标签
  gl.xformatter = LONGITUDE_FORMATTER  #x轴设为经度的格式
  gl.yformatter = LATITUDE_FORMATTER  #y轴设为纬度的格式
  return fig, ax

对数据中感兴趣的区域进行提取并简单的可视化。

代码语言:javascript复制
# 生成地图
fig, ax = map()
# 数据读取及时间平均处理
ds = xr.open_dataset('EC-Interim_monthly_2018.nc')
lat = ds.latitude
lon = ds.longitude
time = ds.time
temp = (ds['t2m'] - 273.15) # 把温度转换为℃
# 区域选择
lon_range = lon[(lon>70) & (lon<140)]
lat_range = lat[(lat>0) & (lat<60)]
temp_region = temp.sel(longitude=lon_range, latitude=lat_range, time='2018-02-01')
# 画图
cbar_kwargs = {
   'label': '2m temperature (℃)',
   'ticks': np.arange(-30, 30 5, 5)
}
levels = np.arange(-30, 30 1, 1)
temp_region.plot.contourf(ax=ax, levels=levels, cmap='Spectral_r', extend='both',
    cbar_kwargs=cbar_kwargs, transform=ccrs.PlateCarree())
fig.show()

0 人点赞