Python版anatool编写记录（一）

第N+N次计划

计划编写一个python版本的anatool，用作模式结果分析，能做到以下几点：

功能要求

1.查看数据内容（包括GrADS格式，grib2格式和二进制bin格式）；
2.根据选择获取相应的要素；
3.对提取出来的要素数据进行处理加工（包括不仅限于，集合平均，离散度，谱分解等等）；
4.绘制相应的图形
5.待增加

性能要求

1.面向对象，需要设置、选择或输入的尽量少；
2.运行速度快；
3.图形美观；
4.有可能的话做成UI交互。

GrADS格式数据信息显示和读取

使用了miniufo大大发布在github上的工具xgrads，传送门
https://github.com/miniufo/xgrads
注意：
1.xgrads的Examples里面调用方法在实际中应该是；
2.xgrads使用了Dask包，可以有效地减少数据读取时间，Dask官方描述自己为:

Dask is a flexible library for parallel computing in Python.
Dask是一个十分灵活的用于Python的并行计算库。

所以在读取数据时需要按照dask.array的要求来读取（还没弄的很清楚）
3.设计的ReadBin包结构如下（注意__init__.py必须有，可以是空文件）：
ReadBin/
├── getvar.py
├── __init__.py
├── __pycache__
└── showctl.py

其中showctl.py为读取CTL描述文件，并显示：

# -*- coding: utf-8 -*-
from xgrads import CtlDescriptor

class showctl:
    '输出CTL文件内容信息'
    def __init__(self,ctlname):
        self.ctl = ctlname

    def ctlinfo(self):
        ctlinfo = CtlDescriptor(file=self.ctl)
        print(ctlinfo)

其中getvar.py为读取GrADS数据，从中提取某要素的单层数据：

# -*- coding: utf-8 -*-
from xgrads import open_CtlDataset

class Grads_data:
    'Extract data from GrADS file with CTL file.'
    def __init__(self,ctlname):
        self.ctl = ctlname

    def getall(self):
        dset = open_CtlDataset(self.ctl)
        return dset

    def getsinglelevel(self,varname,ntime,nlev):
        dset = open_CtlDataset(self.ctl)
        return dset[varname][ntime,nlev,:,:]

编写了一个测试脚本，读取最底层温度，并绘图

# -*- coding: utf-8 -*-
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
import ReadBin.getvar as gv

ctlname = "./database/post.ctl_202006020004400"
ds = gv.Grads_data(ctlname)
ntim = 0
nlev = 25 #这里按照常理应该是0，但是查询读取之后的数据发现垂直层似乎交换了，存疑？
varn = "t"
labeln = "temperature"
t = ds.getsinglelevel(varn,ntim,nlev).compute() - 273.15
lons=t.lon.data
lats=t.lat.data

temp = xr.DataArray(t.data, coords=[lats,lons], dims=['latitude','longitude'])

# 创建画图空间
proj = ccrs.PlateCarree() #创建投影
fig = plt.figure(figsize=(16,9)) #创建页面
ax = fig.subplots(1, 1, subplot_kw={'projection': proj}) #子图
# 设置地图属性:加载国界、海岸线、河流、湖泊
ax.add_feature(cfeat.BORDERS.with_scale('50m'), linewidth=0.8, zorder=1)
ax.add_feature(cfeat.COASTLINE.with_scale('50m'), linewidth=0.6, zorder=1)
ax.add_feature(cfeat.RIVERS.with_scale('50m'), zorder=1)
ax.add_feature(cfeat.LAKES.with_scale('50m'), zorder=1)
# 设置网格点属性
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轴设为纬度格式
# 设置colorbar
cbar_kwargs = {
'orientation': 'horizontal',
'label': labeln,
'shrink': 0.8,
}
# 画图
levels1 = np.arange(0,45,1)
levels2 = np.arange(0,45,5)
a = temp.plot.contourf(ax=ax, levels=levels1, cmap='Spectral_r',cbar_kwargs=cbar_kwargs, transform=ccrs.PlateCarree())
b = temp.plot.contour(ax=ax,levels=levels2,colors=['black'],transform=ccrs.PlateCarree())
plt.clabel(b, inline=True, fontsize=10,fmt='%d')
plt.savefig(labeln+'.jpg')

display  test.jpg

显示图形如下：

temperature.jpg