展示/绘制xarray DataArray时内核重启/崩溃的原因排查

处理含zg变量的CMIP6 CESM2 piControl试验数据（维度：time、levels、lat、lon）时，选取500hPa层级后，执行以下操作会导致内核冻结并重启：

• 计算时间平均并绘图：GPT_500Hpa.mean(dim='time').plot()
• 计算空间平均生成时间序列并绘图：GPT_500Hpa.mean(dim=('lon','lat')).plot()
• 查看数据值：GPT_500Hpa.mean(dim=('lon','lat')).values
• 保存数据到NetCDF文件：GPT_500Hpa.to_netcdf('/CMIP.NCAR.CESM2.piControl.Amon.gn_500Hpa.nc')

但绘制单时间片空间图（GPT_500Hpa[0].plot()）完全正常，此前处理其他CMIP6模型数据未出现该问题，已尝试升级matplotlib、Python、pip，降级freetype，均无效。

代码示例 #

import pandas as pd
import matplotlib.pyplot as plt
import xarray as xr
import numpy as np
import dask
xr.set_options(display_style='html')
import intake
from xmip.preprocessing import rename_cmip6,promote_empty_dims
import cftime

cat_url = "https://storage.googleapis.com/cmip6/pangeo-cmip6.json" #not a local but object storage
col = intake.open_esm_datastore(cat_url)

#Using Geopotential hts for index calculation:
cat  = col.search(experiment_id=['piControl'], table_id=['Amon'], variable_id=['zg'],source_id=['CESM2'],member_id = ['r1i1p1f1'], grid_label=['gn'])

z_kwargs = {'consolidated': True, 'use_cftime':True}
dset_dict = cat.to_dataset_dict(zarr_kwargs=z_kwargs)

#selecting the 500HPa level from the dataset:
dset_dict['CMIP.NCAR.CESM2.piControl.Amon.gn'] = dset_dict['CMIP.NCAR.CESM2.piControl.Amon.gn'].squeeze() #remove empty dims

GPT_500Hpa = dset_dict['CMIP.NCAR.CESM2.piControl.Amon.gn'].zg[:,5,:,:]
GPT_500Hpa[0].plot() #this works!

#Following don't work and kill the kernel:
GPT_500Hpa.mean(dim='time').plot()
GPT_500Hpa.mean(dim=('lon','lat')).plot()

#Not even the .value works:
GPT_500Hpa.mean(dim=('lon','lat')).value

#It doesn't even let me save it to a new netcdf file:
GPT_500Hpa.to_netcdf('/CMIP.NCAR.CESM2.piControl.Amon.gn_500Hpa.nc')

选取500hPa层级后的DataArray结构 #

<xarray.DataArray 'zg' (time: 14400, lat: 192, lon: 288)>
dask.array<rechunk-merge, shape=(14400, 192, 288), dtype=float32, chunksize=(500, 192, 288), chunktype=numpy.ndarray>
Coordinates:
  * lat        (lat) float64 -90.0 -89.06 -88.12 -87.17 ... 88.12 89.06 90.0
  * lon        (lon) float64 0.0 1.25 2.5 3.75 5.0 ... 355.0 356.2 357.5 358.8
    plev       float64 5e+04
  * time       (time) object 0001-01-15 12:00:00 ... 1200-12-15 12:00:00
    member_id  <U8 'r1i1p1f1'
Attributes:
    cell_measures:  area: areacella
    cell_methods:   time: mean
    comment:        Geopotential is the sum of the specific gravitational pot...
    description:    Geopotential is the sum of the specific gravitational pot...
    frequency:      mon
    id:             zg
    long_name:      Geopotential Height
    mipTable:       Amon
    out_name:       zg
    prov:           Amon ((isd.003))
    realm:          atmos
    standard_name:  geopotential_height
    time:           time
    time_label:     time-mean
    time_title:     Temporal mean
    title:          Geopotential Height
    type:           real
    units:          m
    variable_id:    zg

空间平均后的DataArray结构（触发崩溃的对象） #

<xarray.DataArray 'zg' (time: 14400)>
dask.array<mean_agg-aggregate, shape=(14400,), dtype=float32, chunksize=(500,), chunktype=numpy.ndarray>
Coordinates:
    plev       float64 5e+04
  * time       (time) object 0001-01-15 12:00:00 ... 1200-12-15 12:00:00
    member_id  <U8 'r1i1p1f1'
Attributes:
    cell_measures:  area: areacella
    cell_methods:   time: mean
    comment:        Geopotential is the sum of the specific gravitational pot...
    description:    Geopotential is the sum of the specific gravitational pot...
    frequency:      mon
    id:             zg
    long_name:      Geopotential Height
    mipTable:       Amon
    out_name:       zg
    prov:           Amon ((isd.003))
    realm:          atmos
    standard_name:  geopotential_height
    time:           time
    time_label:     time-mean
    time_title:     Temporal mean
    title:          Geopotential Height
    type:           real
    units:          m
    variable_id:    zg

内容的提问来源于stack exchange，提问作者Shreya Trivedi