以下文章来源于MeteoAI ,作者学前班大队长
xarray专题再次开讲,错过第一部分的可以先去补个课从xarray走向netCDF处理(一):数据结构及数据读取。 今天要介绍的就是xarray的索引功能,通过索引你可以对数据进行切片,从整体中提取你所关注的区域、高度或者时间。
在xarray的官方文档中给出了如下几种索引方式
对如下数据进行索引演示:名为ds的DataSet,名为temp的DataArray,数据链接在文末。
位置索引是最直接也是最简单的索引方式,但是位置索引只对DataArray
有效,对DataSet
无效。下面用两种不同方法获取相同的值。
1. 通过数字索引
>>>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. 通过标签索引
>>>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
通过维度的名字就可以不必按照指定的维度顺序对数据进行切片了。对DataArray
和DataSet
都有效,且方法一致。
1.通过数字索引
>>>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.通过标签索引
>>>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...
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函数,该函数主要用来创建地图相关的信息,后面画图的时候可以随时调用。
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
对数据中感兴趣的区域进行提取并简单的可视化。
# 生成地图
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()
本文转自MeteoAI微信公众号,点击文末阅读原文按钮即可跳转原文。