从xarray走向netCDF处理(三)：插值与掩膜

2022-11-02 10:21:55 浏览数 (1)

以下文章来源于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行)。
代码语言：javascript复制
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
数据掩膜

日常数据处理中经常要用到掩盖陆地或者海洋数据，一种方法就是在画图的时候最后添加地理信息，用地理信息掩盖掉海洋或者陆地的数据。这里主要想说的是另一种方法，直接对数据进行处理，把海洋或者陆地区域的值设为缺测。

  1. 首先要下载海陆分布数据landsea.nc：https://apps.ecmwf.int/datasets/data/interim-full-invariant/选择 Land-sea mask下载。
  2. 对任意的DataArray或者Dataset创建一个新的坐标，将海陆数据附给他。
  3. 根据海陆分布数据中海洋或者陆地的值来提取掩膜数据。
代码语言：javascript复制
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
实例
代码语言：javascript复制
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微信公众号，点击文末阅读原文按钮即可跳转原文。

编程算法 数据处理

0 人点赞