【重磅】用NCL的思想理解Python【持续更新】

贫道敬孔

标题党了，其实就是想把NCL中常用的一些函数用python表示出来。

看似很简单，但其实作为新手，无论是搜索还是找例子还是比较难的，本贴是自己新学python的一点心得，和大家分享，不定期持续更新。

有大佬熟悉的欢迎一起分享呀：

1、复制变量属性
     NCL: copy_VarCoords( hgt, hgt_ltm )
     PYTHON：hgt_ltm  = xr.DataArray(hgt_ltm,  coords=hgt.coords, dims=hgt.dims)
     注：在python中目前没有方便的直接赋予“属性”的函数，要借助xr.DataArray。但是，Python中有一点好的是，只要一经这种赋予属性，后面遇到加减乘除乃至微分也不需要再次赋值了，但是众所周知NCL只要经过乘除类的计算，它的属性就自己没有了，还得重新赋予。这一点python要“坚强”的多。对了，SPYDER中可以在右上的Variable Explorer中方便地查看各变量的属性，如果有属性，就在type中显示core.dataarray.DataArray，否则就是简单的Array of float32。spyder是真的强。

2、读完(NC)数据后优雅地选择日期和层数
     NCL:
     TIME          = zfile->time
     YYYYMMDD = cd_calendar(TIME,-2)   
     ymStrt        = 20220501
     ymLast        = 20220630            
     iStrt            = ind(YYYYMMDD.eq.ymStrt)
     iLast           = ind(YYYYMMDD.eq.ymLast)

      hgt             =  zfile->z( iStrt:iLast, {300}, ::-1, : )

      PYTHON:
      time_start = '2022-05-01'
      time_last  = '2022-06-30'


      zfile = xr.open_dataset('F:/hgt.2022.nc')
      hgt   = zfile['z'].sel(time=slice(time_start, time_last), level=300)
      hgt   = hgt[:,::-1,:]
      注：NCL读取NC/GRIB等文件时间设置看似有点麻烦，其实是“一劳永逸”的；python只能用xr.open_dataset读取后用.sel来选取，看似简短了，但是貌似自己生成的NC文件不能用这种方法，所以这种NC数据的日期只能用Dataset读取后“数出”，不能用时间识别出，所以这个功能还是NCL占胜场。

3、扩充维度
     NCL: coslat_4d      = conform_dims( dimsizes(hgt), coslat, 2 )
     PYTHON：coslat_4d = np.expand_dims( coslat , axis = (0,1,3) )
     注1、：二者的思维是一样的，即确定需要【被】扩充的【自己的数组(coslat)】那一维，NCL是给出【目标数组】然后给出【自己的数组(coslat)】的数组所在的维数，python是直接将【自己的数组(coslat)】扩充，但是不需要【目标数组】。二者的不同是python不需要扩充成【目标数组】的具体维度，而是将【目标数组】的其他维扩充成1来“占位”，方便计算。个人感觉这种纯粹为了计算而产生的方法非常直接而简洁，也不占内存，感觉这个python占胜场。
     注2、：扩展维数还有一个np.broadcast_to。但是，np.broadcast_to只能“向前”复制扩展，比如(145)可以np.broadcast_to为(6,145)但是不能np.broadcast_to为(145,6)，所以为了一次性扩展成功，只能将原数组放到最后一维，而想要扩展的数组都放到前面！然后再reshape即可。综上，没有 np.expand_dims好用。

4、求平均
     NCL: psi_d_mean     = dim_avg_n_Wrap(psi_d, 0)
     PYTHON: psi_d_mean = np.mean(psi_d, axis=0)  
     注：这个PYTHON更容易。而且，如上所述，python在这种np下的计算中属性很多时候是保留的，NCL只能看是否有_Wrap的选项了，没有的话就没有属性。python胜。

5、计算cos_phi
     NCL: coslat            = cos( lat(:) * pi / 180. )
     PYTHON: coslat      = np.cos( lat*np.pi / 180. )
     注：涉及到计算，python中要经常用到np（关键python自带的函数不好用），可以看到python中连pi都有数值，NCL还要先赋予pi一个值，所以这个python方便。


6、计算微分
     NCL: dpsi_ddlon       = center_finite_diff_n( psi_d, lon*pi/180., True, 0, 3 )
     PYTHON: dpsi_ddlon = np.gradient(psi_d,dlon[1])[2] # 这里dlon=np.gradient(lon)*np.pi/180.0
     注1、：可以看出python的微分很简洁，也好懂，式子最后的[2]的意思是经度lon在第三维，NCL式子的最后3就是Lon是第四维了。二者等价。算出的值略有差异，但是差别不大，可以放心转换用之。
     注2、：不要用dpsi_ddlon = psi_d_mean.differentiate('longitude')这个来算微分！！算出来量级不对！！


今天就先更新到这里，不对的欢迎大家指出，下次有心得继续更新。

虹hhh

{:5_213:}{:5_213:}{:5_213:}

灭火器

关于 numpy 扩充维度，官方文档介绍了广播的具体规则：两个数组进行运算时形状从右往左对齐，缺失的维度长度先设为 1，然后逐维对齐长度。所以要将形如 (145,) 的数组 a 做 broadcast_to 操作，有两种方法

1. np.broadcast_to(a[:, np.newaxis], (len(a), 6))
   
2. a[:, np.newaxis] * np.ones(6, dtype=a.dtype)

复制代码

lbj2361

楼主总结的很好啊，赞一个

Swaddle

学习啦，很有收获

风叶ele

看起来感觉一旦已经对某种语言深入使用习惯了，再去改用另外一种还是挺费劲儿的

XiaoMaFenJu

灭火器 发表于 2024-1-10 10:04
关于 numpy 扩充维度，官方文档介绍了广播的具体规则：两个数组进行运算时形状从右往左对齐，缺失的维度长 ...

是的，用np.newaxis是比较方便的

XiaoMaFenJu

自己生成的nc其实也可以sel，主要看你的属性信息写的咋样，比如time维你要写入datetime64格式的数据才能和你写的那样去筛选。python去生成nc其实是比较麻烦的，那些coords和Attrs都要自己去写和分配，确实不如ncl我认为

XiaoMaFenJu

最后就是可以去看metpy这个包，常用的气象物理量计算 和 梯度，拉普拉斯算子等都包括了

Reality

赞一个，纠正一个小问题，python自己写得nc也可以通过.sel和.loc来定位，主要是时间time写入时变量类型必须是datetime64