前往小程序,Get更优阅读体验!
立即前往
首页
学习
活动
专区
工具
TVP
发布
社区首页 >专栏 >手把手带你科研入门系列 | PyAOS基础教程十:大数据文件

手把手带你科研入门系列 | PyAOS基础教程十:大数据文件

作者头像
气象学家
发布2021-07-27 15:19:58
1.2K0
发布2021-07-27 15:19:58
举报
文章被收录于专栏:气象学家

1、前言

文章解答以下疑问:

第一:如何在多CMIP6文件的场景下避免内存泄漏。

文章的目标

第一:了解netCDF数据块chunk的概念;

第二:导入dask库,并启动并行处理机制;

第三:计算并绘制高分辨率模型的最大日降雨量。

2、数据处理

首先看一下测试nc文件,总计7个文件,每个文件大约6.7G,是CNRM-CM6-1-HR模式按照25年的时间分开存储的。

由于模式数据非常巨大,一般pc的内存不够大,无法一次性处理如此大的文件,因此这里不再使用xarray库直接读取数据,而是先用glob库,通过glob库提供的方法将上述7个文件导入系统,但这个时候数据还未读取到系统内存

代码语言:javascript
复制
import glob

pr_files = glob.glob('data/pr*.nc')
pr_files.sort()
print(pr_files)

输出:

代码语言:javascript
复制
['/Users/irving/Desktop/data-carpentry/data/pr_day_CNRM-CM6-1-HR_historical_r1i1p1f2_gr_18500101-18741231.nc',
 '/Users/irving/Desktop/data-carpentry/data/pr_day_CNRM-CM6-1-HR_historical_r1i1p1f2_gr_18750101-18991231.nc',
 '/Users/irving/Desktop/data-carpentry/data/pr_day_CNRM-CM6-1-HR_historical_r1i1p1f2_gr_19000101-19241231.nc',
 '/Users/irving/Desktop/data-carpentry/data/pr_day_CNRM-CM6-1-HR_historical_r1i1p1f2_gr_19250101-19491231.nc',
 '/Users/irving/Desktop/data-carpentry/data/pr_day_CNRM-CM6-1-HR_historical_r1i1p1f2_gr_19500101-19741231.nc',
 '/Users/irving/Desktop/data-carpentry/data/pr_day_CNRM-CM6-1-HR_historical_r1i1p1f2_gr_19750101-19991231.nc',
 '/Users/irving/Desktop/data-carpentry/data/pr_day_CNRM-CM6-1-HR_historical_r1i1p1f2_gr_20000101-20141231.nc']

然后,用xarray读取数据,但是这里读取数据的方法,与前面的课程有非常明显的不同(前面用的是xarray.open_dataset来一次性读取nc文件到内存中),这里用到的是xarray.open_mfdataset函数分批读取数据,我们具体来看看它是如何读取数据的。

代码语言:javascript
复制
import xarray as xr

dset = xr.open_mfdataset(pr_files, chunks={'time': '500MB'})
print(dset)

pr_files即为上面的glob抓取的文件,虽说glob一次性抓取了7个nc文件,但是这里xarray读取依然类似于一个文件,参数chunks(数据块)是一个关键,这里的意思是在time维度上一次性读取500MB的数据块,实现按需读取数据。

输出:

代码语言:javascript
复制
<xarray.Dataset>
Dimensions:      (axis_nbounds: 2, lat: 360, lon: 720, time: 60265)
Coordinates:
  * lat          (lat) float64 -89.62 -89.12 -88.62 -88.13 ... 88.62 89.12 89.62
  * lon          (lon) float64 0.0 0.5 1.0 1.5 2.0 ... 358.0 358.5 359.0 359.5
  * time         (time) datetime64[ns] 1850-01-01T12:00:00 ... 2014-12-31T12:...
Dimensions without coordinates: axis_nbounds
Data variables:
    time_bounds  (time, axis_nbounds) datetime64[ns] dask.array<chunksize=(9131, 2), meta=np.ndarray>
    pr           (time, lat, lon) float32 dask.array<chunksize=(397, 360, 720), meta=np.ndarray>
Attributes:
    Conventions:            CF-1.7 CMIP-6.2
    creation_date:          2019-05-23T12:33:53Z
    description:            CMIP6 historical
    title:                  CNRM-CM6-1-HR model output prepared for CMIP6 and...
    activity_id:            CMIP
    contact:                contact.cmip@meteo.fr
    data_specs_version:     01.00.21
    dr2xml_version:         1.16
    experiment_id:          historical
    experiment:             all-forcing simulation of the recent past
    external_variables:     areacella
    forcing_index:          2
    frequency:              day
    further_info_url:       https://furtherinfo.es-doc.org/CMIP6.CNRM-CERFACS...
    grid:                   data regridded to a 359 gaussian grid (360x720 la...
    grid_label:             gr
    nominal_resolution:     50 km
    initialization_index:   1
    institution_id:         CNRM-CERFACS
    institution:            CNRM (Centre National de Recherches Meteorologiqu...
    license:                CMIP6 model data produced by CNRM-CERFACS is lice...
    mip_era:                CMIP6
    parent_experiment_id:   p i C o n t r o l
    parent_mip_era:         CMIP6
    parent_activity_id:     C M I P
    parent_source_id:       CNRM-CM6-1-HR
    parent_time_units:      days since 1850-01-01 00:00:00
    parent_variant_label:   r1i1p1f2
    branch_method:          standard
    branch_time_in_parent:  0.0
    branch_time_in_child:   0.0
    physics_index:          1
    product:                model-output
    realization_index:      1
    realm:                  atmos
    references:             http://www.umr-cnrm.fr/cmip6/references
    source:                 CNRM-CM6-1-HR (2017):  aerosol: prescribed monthl...
    source_id:              CNRM-CM6-1-HR
    source_type:            AOGCM
    sub_experiment_id:      none
    sub_experiment:         none
    table_id:               day
    variable_id:            pr
    variant_label:          r1i1p1f2
    EXPID:                  CNRM-CM6-1-HR_historical_r1i1p1f2
    CMIP6_CV_version:       cv=6.2.3.0-7-g2019642
    dr2xml_md5sum:          45d4369d889ddfb8149d771d8625e9ec
    xios_commit:            1442-shuffle
    nemo_gelato_commit:     84a9e3f161dade7_8250e198106a168
    arpege_minor_version:   6.3.3
    history:                none
    tracking_id:            hdl:21.14100/223fa794-73fe-4bb5-9209-8ff910f7dc40

从第1行输出信息来看,dset依然是xarray.Dataset类型的变量,请注意看第9和10行的变量中新增的dask.array对象下的chunksize属性,这是由于我们在读取dset数据时指定chunk参数的原因。按照chunk参数指定的500MB的大小,dask并非将7个nc文件的数据一次性读取到系统内存中,而是遵从一块一块数据读取的原则。当然dask也可以把这些chunks分发到不同的cpu核上进行处理。

那么多大的chunk比较合适呢?如果chunk太小,频繁的调度数据并处理数据将导致效率低下,整体耗时可能依然比较高;如果chunk太大,可能会导致系统运行缓慢,甚至内存泄漏。因此chunk既不能太大,也不能太小,dask的官方文档中给的推荐值是10MB-1GB,比如上面的例子中就是选用的中间值500MB的chunk。

接下来,看一下上面读取的数据大小

代码语言:javascript
复制
dset['pr'].data

输出:

代码语言:javascript
复制
           Array                 Chunk
Bytes     62.48 GB              499.74 MB
Shape     (60265, 360, 720)     (482, 360, 720)
Count     307 Tasks             150 Chunks
Type      float32               numpy.ndarray

最后,按照时间顺序计算日最大降雨量

代码语言:javascript
复制
pr_max = dset['pr'].max('time', keep_attrs=True)
print(pr_max)

输出:

代码语言:javascript
复制
<xarray.DataArray 'pr' (lat: 360, lon: 720)>
dask.array<nanmax-aggregate, shape=(360, 720), dtype=float32, chunksize=(360, 720), chunktype=numpy.ndarray>
Coordinates:
  * lat      (lat) float64 -89.62 -89.12 -88.62 -88.13 ... 88.62 89.12 89.62
  * lon      (lon) float64 0.0 0.5 1.0 1.5 2.0 ... 357.5 358.0 358.5 359.0 359.5
Attributes:
    long_name:           Precipitation
    units:               kg m-2 s-1
    online_operation:    average
    cell_methods:        area: time: mean
    interval_operation:  900 s
    interval_write:      1 d
    standard_name:       precipitation_flux
    description:         at surface; includes both liquid and solid phases fr...
    history:             none
    cell_measures:       area: areacella

上面的计算过程看上去是在很短的时间里就完成了,但实际上它依然是xarray懒人模式的一种,一般来说,xarray非必要的情况下不会计算,但是绘图或者写入netCDF文件则会发生计算操作。那么有没有办法强制xarray进行数据计算呢?办法当然是有的,computer函数就可以实现此目的。

代码语言:javascript
复制
%%time
pr_max.compute()

第一行代码的作用是打印当前cell的运行时间。

输出:

代码语言:javascript
复制
CPU times: user 4min 1s, sys: 54.2 s, total: 4min 55s
Wall time: 3min 44s

3、并行化

上面的例子中,所有的计算处理都是运行在单核上,而dask client可以把任务分发至不同的cpu核上,实现并行化处理。使用方法如下:

代码语言:javascript
复制
from dask.distributed import Client

client = Client()
client

输出:

代码语言:javascript
复制
Client                                    Cluster
Scheduler: tcp://127.0.0.1:59152          Workers: 4
Dashboard: http://127.0.0.1:8787/status   Cores: 4
                                          Memory: 17.18 GB

然后,我们再来调用一下computer函数,来看看数据处理花了多少时间,跟前面的单核场景进行对比:

代码语言:javascript
复制
%%time
pr_max.compute()

输出:

代码语言:javascript
复制
CPU times: user 10.2 s, sys: 1.12 s, total: 11.3 s
Wall time: 2min 33s

从这个结果中,可以看到,虽然是4个cpu核参加数据处理,整个cell的运行时间是2min33s,但跟前面单核处理时间3min44s,并没有减少75%的运行时间。说明在多核cpu之间进行系统调度也是耗费时间的,因此,多核cpu并行处理化场景可能不是最优解决方案,需要根据实际情况选择方案。

4、绘图

在完成了日最大降雨量的数据计算后,即可以完成画图工作。

代码语言:javascript
复制
import matplotlib.pyplot as plt
import cartopy.crs as ccrs
import numpy as np
import cmocean

pr_max.data = pr_max.data * 86400
pr_max.attrs['units'] = 'mm/day'

fig = plt.figure(figsize=[12,5])
ax = fig.add_subplot(111, projection=ccrs.PlateCarree(central_longitude=180))
pr_max.plot.contourf(ax=ax,
                     levels=np.arange(0, 450, 50),
                     extend='max',
                     transform=ccrs.PlateCarree(),
                     cbar_kwargs={'label': pr_max.units},
                     cmap=cmocean.cm.haline_r)
ax.coastlines()

model = dset.attrs['source_id']
title = f'Daily maximum precipitation, 1850-2014 ({model})'
plt.title(title)

plt.show()

5、总结

本文的主要知识点:

  • 学会用dask和xarray库让netCDF数据加载、处理和可视化等操作更加简单;
  • Dask可以通过并行加速数据处理,但需要特别注意数据分块大小。
本文参与 腾讯云自媒体同步曝光计划,分享自微信公众号。
原始发表:2021-07-25,如有侵权请联系 cloudcommunity@tencent.com 删除

本文分享自 气象学家 微信公众号,前往查看

如有侵权,请联系 cloudcommunity@tencent.com 删除。

本文参与 腾讯云自媒体同步曝光计划  ,欢迎热爱写作的你一起参与!

评论
登录后参与评论
0 条评论
热度
最新
推荐阅读
目录
  • 1、前言
  • 2、数据处理
  • 3、并行化
  • 4、绘图
  • 5、总结
领券
问题归档专栏文章快讯文章归档关键词归档开发者手册归档开发者手册 Section 归档