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

贫道敬孔

python在使用
from scipy.stats import t
计算t检验的值时，

value_confi = t.sf(value_t_test, df=time-1) 仅对正值（仅对目标value正值做检验，负值不做），
value_confi = t.cdf(value_t_test, df=time-1) 仅对负值（同上但是仅对目标value负值），
（即上述t.sf和t.cdf都是单边检验）
而想达到NCL中value_confi  = student_t( value_t_test, time-1 )的效果，
必须使用2 * (1 - t.cdf(np.abs(value_t_test), df=time-1)) 才可同时对目标value正值和负值都做t检验
注意，这里同时value同时有正负值时不能使用2 * (1 - t.sf(np.abs(value_t_test), df=time-1)) ！

贫道敬孔

画两点之间的剖面图并叠加地型阴影：NCL能跑出来且直观，PYTHON的还没弄出来，弄出来更新。

1、剖面图
核心思想：使用gc_latlon对指定两个点之间的经纬度进行插值，然后使用linint2_points对想画剖面的变量（使用已经插值好的经纬度）进行差值


start_lat = 12.5        ; 设置起始点经纬度
start_lon = 87.5
end_lat       = 35.0
end_lon       = 115.0
npts              = 21; 使差值后的经纬度的范围在1.25左右即可，可打印后修改，具体数值不太重要

lat     = w_mean&latitude
lon     = w_mean&longitude

dist    = gc_latlon( start_lat, start_lon, end_lat, end_lon, npts, 2 )  
printVarSummary(dist)
lat_new = dist@gclat
lon_new = dist@gclon
copy_VarAtts( lat, lat_new )
copy_VarAtts( lon, lon_new )
w_new     = linint2_points(lon, lat, w_mean, False, lon_new, lat_new, 0)    ;将要绘制的变量进行插值
uvd_v_new = linint2_points(lon, lat, uvd_v, False, lon_new, lat_new, 0) ;注这里用辐散风和w画剖面，更准确
topo_new  = linint2_points(lon, lat, topography, False, lon_new, lat_new, 0)

w_new          =  -1 * w_new * scale
w_new_verse    =  -1 * w_new   

w_new!0        = "level"
w_new&level    = w_mean_1958_2023&level
w_new!1        = "latitude"
w_new&latitude = lat_new
copy_VarCoords( w_new, q_new )
copy_VarCoords( w_new, w_new_verse )
copy_VarCoords( w_new, uvd_v_new )

print("lon_new is " + lon_new )
print("lat_new is " + lat_new )
topo_new!0          = "latitude"
topo_new&latitude   = lat_new

注意这里1不能使用linint2_points_Wrap计算，因为这样变量维度就不是目标的level和latitude了，所以要手动赋予坐标信息；2是地型topo_new不能简单地copy_VarCoords( w_new(:,0), topo_new )，这样会出错！


2、地型数据处理
核心思想：使用压高公式！
;;;;;; 用压高公式p=p0*(1-Τ*hgt/(273.15+15))^(g/(R*Τ))，把标准大气的高度转气压。然后画topo2的填色图。
;;;;;; 这里的T即Γ，即温度垂直递减率，每上升单位高度温度的降低量。单位：K/m（如国际标准大气中 Γ≈0.0065K/m）
topo2 = 1013.25 * (1 - topo_new*0.0065/288.15) ^ 5.25145
print(topo2)
copy_VarCoords( topo_new, topo2 )
printVarSummary(topo2)
print( "The max value of topo2 is " + max(topo2) )
print( "The min value of topo2 is " + min(topo2) )  

直接用，没问题的。

3、画图的话就用正常的剖面图如我用的gsn_csm_pres_hgt_streamline，但是要注意的有：

1）所有非地型数据一律要predraw，比如
res@cnFillDrawOrder        = "PreDraw" ;;;; 非常重要！不然后面地型的黑色不能显示   
res@stStreamlineDrawOrder  = "PreDraw" ;;;; 非常重要！不然后面地型的黑色不能显示
res_line@cnLineDrawOrder         = "PreDraw" ;;;; 非常重要！不然后面地型的黑色不能显示
res_line@cnLabelDrawOrder        = "PreDraw" ;;;; 非常重要！不然后面地型的黑色不能显示  
注意这里连等值线的标签都要predraw，不然标签又会在地型上！

2）想再叠加其他变量的等值线直接在前面和w_mean一样处理好相应的剖面数据后，使用plot_line = gsn_csm_contour(wks, q_new, res_line) 即可

3）地形画图
;;;;添加地形      
res_mask                        = True
res_mask@gsnDraw                = False
res_mask@gsnFrame               = False   
res_mask@gsnRightString         = " "   
res_mask@gsnLeftString          = " "
res_mask@xyLineThicknesses      = (/3.0/)              ; make line thicker
res_mask@xyLineColors           = (/"green"/)          ; change line color
res_mask@gsnYRefLine            = (/ min(topo2) , max(topo2) /)
res_mask@gsnBelowYRefLineColor  = (/ "gray", "gray"/)

plot_mask   = gsn_csm_xy( wks, topo2&latitude, topo2, res_mask )              
overlay( plot, plot_mask )

注意1是用xy画图！因为地型本来就是二维数据，再剖面后其实就剩一维序列了，只能用xy，不能用counter画！也不能用gsn_csm_pres_hgt画！这二者都要data是二维！
2是这里不用使用res_mask@xyCurveDrawOrder       = "PostDraw",因为前面已经使用过predraw了。

以上就是NCL画任意两点的剖面叠加地型的程序了，只要数据对绝对能跑通。python用metpy的cross函数总是报错，都不能直接取点，等能跑通的时候更新。


参考了以下帖子的代码和热心人回复中的思路和代码，统一致谢：
http://bbs.06climate.com/forum.p ... hlight=%C6%CA%C3%E6
http://bbs.06climate.com/forum.p ... hlight=%C6%CA%C3%E6
http://bbs.06climate.com/forum.p ... age=1&authorid=9306
http://bbs.06climate.com/forum.p ... %C6%CA%C3%E6&page=2
http://bbs.06climate.com/forum.p ... hlight=%C6%CA%C3%E6

贫道敬孔

NCL 画XY图（如gsn_csm_xy）

1、控制坐标绝对值是tr，而不是tmXB，因为后者“管的更宽”，前者只针对坐标值本身。
如一组序列范围在-8到8不等，如果只设置
res@tmXBMode              = "Explicit"   ; 注意不要用"Manual" ;
res@tmXBValues            = ispan(-8, 8, 2)  
res@tmXBLabels            = ispan(-8, 8, 2)  
而不设置tr的话，x的范围（图片的显示范围）会根据其实际大小而改变，即不受标签范围的控制

此时如想“精准控制”图片范围，还得设置
res@trXMinF               = -8
res@trXMaxF               = 8
;res@trYMinF                 = -10.0           ; min value on y-axis
;res@trYMaxF                 =  100             ; max value on y-axis

2、参考线既可以使用
res@gsnYRefLine             = 0.   
也可使用
lnres1                       = True      
lnres1@gsLineThicknessF      = 5         
lnres1@gsLineDashPattern     = 12  
lnres1@gsLineColor           = "gray35"
  
ref_plot_line1               = gsn_add_polyline( wks, plot(1), (/0, 0/), (/0, 1000/), lnres1 )

3、只使用gsn_csm_y画图的话，默认序列的“值”是Y轴，其属性（如level或time）是X轴，此时想翻转XY轴的显示的话，只能用gsn_csm_xy，并指定属性做Y轴
plot(0) = gsn_csm_xy(wks, q_mean, q_mean&level, res0)

注意ERA5的数据层数是默认125-1000hPa的，翻转数据不用在level维::-1，而是
res@trYReverse             = True
即可

贫道敬孔

NCL 给显示的title加平方：

~S~要加的数字或字母~N~，如~S~2~N~。注意S和N前后都要有波浪号~，否则报错！

PYTHON 给显示的title加平方：
$^要加的数字或字母$，如$^2$。注意这里似乎不能给两个不同的“属性”同时使用，如想要对-5做幂，NCL直接~S~-5~N~即可，而python只弄一次不行，必须要$^-$$^5$才可显示。

注：python显示psi的一瞥 ' 的命令：$\mathcal{\psi}$$^\prime$

贫道敬孔

这里注意对右下角的符号~B~**~N~后想要再连接数字/字母，要有两个波浪号，即~B~2R~N~~，而不是一个，如对-Q2R/L进行显示，则代码为"-~F34~a~F25~Q~B~2R~N~~F34~q~F25~/L"，而不是"-~F34~a~F25~Q~B~2R~N~F34~q~F25~/L"，否则就是个错。
可参考https://bbs.06climate.com/forum. ... hlight=%B7%FB%BA%C5
进行对比

贫道敬孔

NCL中的倒三角的符号可以有很多表达，比如F34的Q，即F34~Q；而三角或者说叫Δ，目前只在~F33~D中发现

贫道敬孔

NCL使用多个子图，想在已有的res情况下，对另外的子图赋值时，当使用
res_QL                       = res
res_QL@cnLevelSpacingF       = 0.03            
res_QL@cnMinLevelValF        = -0.3            
res_QL@cnMaxLevelValF        = 0.3

而出现fatal:["NclAtt.c":262]:Attribute assignment type mismatch.
这种错误时

应该
res_QL@cnLevelSpacingF       := 0.3            
res_QL@cnMinLevelValF        := -3            
res_QL@cnMaxLevelValF        := 3

即对cnLevelSpacingF等与原res“矛盾”的值进行“强赋值”即可

官网教程https://mailman.ucar.edu/pipermail/ncl-talk/2018-May/012144.html