- 积分
- 19044
- 贡献
-
- 精华
- 在线时间
- 小时
- 注册时间
- 2015-4-23
- 最后登录
- 1970-1-1
|

楼主 |
发表于 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
|
|