以下文章来源于MeteoAI ,作者学前班大队长
对于xarray之前已经介绍过两期了,分别是数据结构及数据读取和数据索引。 这一期要介绍的功能是插值与掩膜。 这两个方法在数据处理中会经常用到,实用等级☆☆☆☆☆。
xarray中对scipy的插值函数进行了进一步的封装,可以让我们方便的调用。
只需要对DataArray
,DataSet
使用interp()
函数就可以实现插值了,就像索引一样简单。不管是一维数据还是多维数据都可以轻松搞定。
下面是官方给出的例子,DataArray
的时间维度总共有四个值[0,1,2,3]。
da.sel(time=3)
,索引时间维的值为3(12行)。da.interp(time=2.5)
,将时间维原本不存在的2.5插值了出来(22行)。In [1]: da = xr.DataArray(np.sin(0.3 * np.arange(12).reshape(4, 3)),
...: [('time', np.arange(4)),
...: ('space', [0.1, 0.2, 0.3])])
...:
# label lookup
In [2]: da.sel(time=3)
Out[2]:
<xarray.DataArray (space: 3)>
array([ 0.42738 , 0.14112 , -0.157746])
Coordinates:
time int64 3
* space (space) float64 0.1 0.2 0.3
# interpolation
In [3]: da.interp(time=2.5)
Out[3]:
<xarray.DataArray (space: 3)>
array([0.700614, 0.502165, 0.258859])
Coordinates:
* space (space) float64 0.1 0.2 0.3
time float64 2.5
日常数据处理中经常要用到掩盖陆地或者海洋数据,一种方法就是在画图的时候最后添加地理信息,用地理信息掩盖掉海洋或者陆地的数据。这里主要想说的是另一种方法,直接对数据进行处理,把海洋或者陆地区域的值设为缺测。
DataArray
或者Dataset
创建一个新的坐标,将海陆数据附给他。def mask(ds, label='land'):
landsea = xr.open_dataset('landsea.nc')['LSMASK']
ds.coords['mask'] = (('latitude', 'longitude'), landsea.values)
if label == 'land':
ds = ds.where(ds.mask < 1)
elif label == 'ocean':
ds = ds.where(ds.mask > 0)
return ds
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
def create_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
def mask(ds, label='land'):
landsea = xr.open_dataset('landsea.nc')
landsea = landsea['LSMASK']
# --ds和地形数据分辨率不一致,需将地形数据插值
landsea = landsea.interp(lat=ds.latitude.values, lon=ds.longitude.values)
# --利用地形掩盖海陆数据
ds.coords['mask'] = (('latitude', 'longitude'), landsea.values)
print(ds.mask)
if label == 'land':
ds = ds.where(ds.mask < 0.8)
elif label == 'ocean':
ds = ds.where(ds.mask > 0.2)
return ds
if __name__ == '__main__':
# 数据读取及时间平均处理
fig, ax = create_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')
temp_mask = mask(temp_region, 'ocean')
# --画图
cbar_kwargs = {'label': '2m temperature (℃)',
'ticks': np.arange(-30, 30+5, 5)
本文转自MeteoAI微信公众号,点击文末阅读原文按钮即可跳转原文。