【重磅】用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，所以要有一个单位的转换