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

Lancelot

实际上xarray只是给numpy array加了一层包装。许多讨论不如直接把里面numpy array的部分抽出来，直接用矩阵计算处理，有时比xarray方便很多。

puck66

很好的主题，墙裂支持。。蹲住

puck66

贫道敬孔 发表于 2024-7-18 21:32
python的几大坑：

1、使用np.std时，一定要在括号的最后有, ddof=1，否则会出现后续计算错的离谱比如相 ...

matplotlib绘图这里有点没看明白，，那不叠加的情况是啥样子？

sensensen

收藏，感谢楼主。

贫道敬孔

继续更：

对于有“属性”的变量，比如已经读取的nc文件的u风场，变量名为u。现想要提取其中的level属性,命名为p

PYTHON有两种方法：
1、p = u.coords['level']
2、p = u.level
第二种是不是超级简洁！

NCL:
p  = u&level(::-1) ; 这里的(::-1)仅仅是颠倒层数的意思

总结：PYTHON的“属性”用点即. , 而NCL的属性用连字符即&

贫道敬孔

想获取数组某一维的长度：如想知道一个(time,level,lat,lon)的数组a纬度的格点数

NCL：
先用dimsizes(a)获得其真实维度的信息，然后想要a的哪一维度是哪一维
dimx   = dimsizes(a)
ntim   = dimx(0)
nlev   = dimx(1)
nlat   = dimx(2)
nlon   = dimx(3)  

PYTHON:
a.latitude[::-1]
nlat = len(lat)

二者均需得到所在维后再用各自的方法获得长度

贫道敬孔

垂直积分：

python中numpy和scipy的计算垂直积分的效果相同：

p = T['level'].metpy.quantify() * 100 * units.Pa  # 假设level单位为hPa（p一定要乘以100！化为帕！！）

Q1_integrated = np.trapz( Q1, p, axis=0 ) / g # 这里p即x=p的简写；与不乘以100即hPa的结果相同
Q2_integrated = np.trapz( Q2, p, axis=0 ) / g

与
from scipy import integrate
Q1_integrated = integrate.trapezoid(Q1, p, axis=0) / g # 这里p亦即x=p的简写；与不乘以100即hPa的结果相同
Q2_integrated = integrate.trapezoid(Q2, p, axis=0) / g
   
二者计算结果一模一样，当然，scipy还有一个integrate.simpson，与二者的计算结果近似，但是值稍微有点差别。

至于dataarray的直接积分，即
Q1 = xr.DataArray(Q1, coords=u.coords, dims=u.dims)
Q2 = xr.DataArray(Q2, coords=u.coords, dims=u.dims)
Q1_integrated = Q1.integrate("level") / g
Q2_integrated = Q2.integrate("level") / g

则算出的结果非常之小，要小两个量级，不知为何（可能是积分直接对百帕而不是帕？但是一个DataArray已经定义好了维度的名称（如level），又该如何修改呢？）。所以推荐前两个。

以方便程度来说推荐第一个即np.trapz

注：最新的numpy已将np.trapz经更新为np.trapezoid，二者调用和计算结果相同，注意版本即可。老一点的numpy兼容性高一点，np.trap好用的。

贫道敬孔

贫道敬孔 发表于 2025-3-10 20:47
垂直积分：

python中numpy和scipy的计算垂直积分的效果相同：

破案了，在这个问题中
至于dataarray的直接积分，即
Q1 = xr.DataArray(Q1, coords=u.coords, dims=u.dims)
Q2 = xr.DataArray(Q2, coords=u.coords, dims=u.dims)
Q1_integrated = Q1.integrate("level") / g
Q2_integrated = Q2.integrate("level") / g
则算出的结果非常之小，要小两个量级，不知为何（可能是积分直接对百帕而不是帕？但是一个DataArray已经定义好了维度的名称（如level），又该如何修改呢？）

直接给
Q1_integrated = Q1.integrate("level")*100 / g
Q2_integrated = Q2.integrate("level")*100 / g
即在积分完后给数值乘以100即可，原因应该就是dataarray对垂直积分直接是对层数即hPa，而不是对Pa，所以要有一个单位的转换

贫道敬孔

python用np.trapz和NCL用vibeta积分原理上结果都一样，但是要注意3点：

1、一律都不要化作dp，而是直接用气压层p去计算！！
对你没有看错，什么NCL的dp = dpres_plevel_Wrap(p, pbot, ptop, 0)
或者python仿照这个函数做的
def dpres_plevel(p_hpa):
    dp = (p_hpa[:-1] - p_hpa[1:]) * 100
    return dp

p_sorted = np.array([1000, 925, 850, 700, 600, 500, 400, 300, 250, 200, 150, 125, 100])
dp = dpres_plevel(p_sorted)
print(dp)
都不要！！一律都直接用气压层直接计算！

2、建议都直接用hPa计算，方便直观，不用乘以100！pattern也不会错！！

3、python直接计算即可，无所谓气压层的方向（1000-125或125-1000hPa），而NCL的vibeta中需要的变量必须都是从地面到高空！！这里NCEP的数据无所谓，ERA5的数据是从125-1000hPa的，必须转置方向，否则就是错的！
注意这里python如果数据本身是125-1000hPa，那p和数据都不用专门转置方向，直接用即可，要是专门转置算出来也是错的！

xzheng2024

这里的np.trapz使用的时候如果气压层方向是递减的为什么不用转置方向啊，我看np.trapz函数，如果维度是递减的，算出来结果是和递增结果相反