爱气象,爱气象家园! 

气象家园

 找回密码
 立即注册

QQ登录

只需一步,快速开始

搜索
楼主: 贫道敬孔

[经验总结] 【重磅】用NCL的思想理解Python【持续更新】

  [复制链接]
 楼主| 发表于 2025-4-4 13:53:03 | 显示全部楼层
本帖最后由 贫道敬孔 于 2025-4-6 09:33 编辑

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)) !
密码修改失败请联系微信:mofangbao
回复 支持 反对

使用道具 举报

 楼主| 发表于 2025-4-8 20:55:24 | 显示全部楼层
本帖最后由 贫道敬孔 于 2025-4-8 20:58 编辑

画两点之间的剖面图并叠加地型阴影: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



密码修改失败请联系微信:mofangbao
回复 支持 反对

使用道具 举报

 楼主| 发表于 2025-4-9 20:24:18 | 显示全部楼层
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
即可
密码修改失败请联系微信:mofangbao
回复 支持 反对

使用道具 举报

 楼主| 发表于 2025-4-9 20:26:47 | 显示全部楼层
本帖最后由 贫道敬孔 于 2025-4-10 09:53 编辑

NCL 给显示的title加平方:

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

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

注:python显示psi的一瞥 ' 的命令:$\mathcal{\psi}$$^\prime$
密码修改失败请联系微信:mofangbao
回复 支持 反对

使用道具 举报

 楼主| 发表于 5 天前 | 显示全部楼层
这里注意对右下角的符号~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
进行对比
密码修改失败请联系微信:mofangbao
回复 支持 反对

使用道具 举报

 楼主| 发表于 2 小时前 | 显示全部楼层
NCL中的倒三角的符号可以有很多表达,比如F34的Q,即F34~Q;而三角或者说叫Δ,目前只在~F33~D中发现
密码修改失败请联系微信:mofangbao
回复 支持 反对

使用道具 举报

您需要登录后才可以回帖 登录 | 立即注册

本版积分规则

Copyright ©2011-2014 bbs.06climate.com All Rights Reserved.  Powered by Discuz! (京ICP-10201084)

本站信息均由会员发表,不代表气象家园立场,禁止在本站发表与国家法律相抵触言论

快速回复 返回顶部 返回列表