[原创]p坐标地形垂直剖面的MATLAB实现

沙颖凯

分享一个在p坐标下绘制变量的垂直剖面时叠加地形的方法填补MATLAB处理此类问题的空白。

处理地形剖面最核心的问题是如何在p坐标下确定海拔高度的对应值，因为大气的气压是不断变化的，因此以气压为参照的p坐标其实也没有完全固定的值，但是大气气压的变化可以视为均值叠加扰动的结果，在取平均值得过程中，扰动能够被略去。在此基础上能够得到相对准确的海拔高度对应的p坐标值。

按上面的思路这里用NCEP/NCAR在1981-2010年的气候态等压面位势高度和etopo的地形计算p坐标地形高度，具体做法是把再分析数据插值到与地形数据相同的分辨率，然后用数据中每个等压面层次（1000, 850, ... 100 hPa）和其对应的位势高度（m）构建一组关系，用线性回归估计地形数据中海拔高度（m）所对应的p坐标值（hPa）。这里提供一份1×1 deg的计算结果 Topo_Barotropic_Low.mat (230.02 KB, 下载次数: 214)

2014-6-14 11:54 上传

点击文件名下载附件 。本来还有一份高分辨率的（0.07×0.07 deg），但是好像太大上传不了，就不提供了。还有就是海洋/湖泊等的区域全部被我定义为缺测了，使用的时候按情况处理。

有了p坐标下的海拔高度对应值后，在处理垂直剖面时加入地形就很方便了，思路和GrADS画剖面一样，对三维数据固定一维，保留垂直方向的变化再操作就行了。这里另外提供一个函数，功能是从数据中提取特定经纬度范围的结果，这样能提取特定范围的数据减少计算量：
GetLatLonInfo.m (1.7 KB, 下载次数: 196)

2014-6-14 11:55 上传

点击文件名下载附件

[rowLIM colLIM]=GetLatLonInfo(longitude, lonlim, latitude, latlim);
具体用法下面的例子中能见到。

为了体现出地形的效果，这里以某西南涡过程沿30N, 95-115E的相对涡度剖面为例， Profile_curlz_Sichuan.mat (5.56 KB, 下载次数: 151)

2014-6-14 11:56 上传

点击文件名下载附件 数据包含算好的涡度剖面结果（来自FNL 1×1 deg），对应的经度和p坐标高度还有一个专用的色标：

1. clear all
   
2. close all
   
3. %% Import data
   
4. load Profile_curlz_Sichuan
   
5. load Topo_Barotropic_Low
   
6. %% Label
   
7. label=[-14 -12 -10 -8 -6 -4 -2 2 4 6 8 10 12 14];
   
8. %% Geo-Info.
   
9. [rowLIM colLIM]=GetLatLonInfo(longitude_topo, [90 120], latitude_topo, [30 30]);
   
10. topo_Sichuan=Topo_Barotropic_Low(rowLIM(1), colLIM(1):colLIM(2));
    
11. longitude_topo=longitude_topo(rowLIM(1), colLIM(1):colLIM(2));
    
12. topo_Sichuan(isnan(topo_Sichuan))=1013.25;
    
13. %% Plot
    
14. fig=figure('Color', 'w');
    
15. hold on
    
16. contourf(lonDIM, hPaDIM, curlz_Sichuan, label, 'LineStyle', ':', 'EdgeColor', [0.5 0.5 0.5])
    
17. contour(lonDIM, hPaDIM, curlz_Sichuan, label(label>0), 'LineStyle', '-', 'EdgeColor', [0.5 0.5 0.5])
    
18. area(longitude_topo, topo_Sichuan, 'EdgeColor', 'none', 'FaceColor', 'k', 'BaseValue', 1013.25)
    
19. hold off
    
20. box on; grid off
    
21. caxis([min(label) max(label)]')
    
22. set(gca, 'YDir', 'reverse', 'FontSize', 8, 'FontWeight', 'bold', 'LineWidth', 1.5,...
    
23.     'XTickLabel', ['95' char(176) 'E' '|' '100' char(176) 'E' '|' '105' char(176) 'E' '|' '110' char(176) 'E' '|' '115' char(176) 'E'])
    
24. %    'YTickLabel', '100 hPa|200 hPa|300 hPa|400 hPa|500 hPa|600 hPa|700 hPa|800 hPa|900 hPa|1000 hPa')
    
25. %xlabel('\fontsize{8}Longitude (^\circE )', 'FontWeight', 'bold')
    
26. ylabel('\fontsize{8}Pressure level ( hPa )', 'FontWeight', 'bold')
    
27. title('\fontsize{10}Jul 21^{st} 2008 08:00 CST', 'FontWeight', 'bold')
    
28. colormap(CMap)
    
29. Handle=colorbar('EastOutside');
    
30. set(get(Handle, 'YLabel'), 'string', '\fontsize{8}Relative Vorticity ( 10^{-5} s^{-1} )', 'FontWeight', 'bold')
    
31. set(Handle, 'YTick', label, 'TickLength', [0.05 0.10], 'LineWidth', 1.0, 'FontWeight', 'bold', 'FontSize', 8)

复制代码

结果大概是这样，代码包含了一些图形句柄的操作，不要也无所谓，不过注意set(gca, 'YDir', 'reverse')要有，填充地形时area的'BaseValue'要修改。

freekiller

画得不错，谢谢分享

lqouc

嗯嗯，很不错啊~~~话说地形上的等值线要是去掉会不会更好呢？

river

感谢分享。楼主辛苦了

逆手一文字

必须赞一个

sfy0717

赞一个   楼主辛苦

quentin011099

如果该1*1地形数据能够用于GrADS作图就好了，请问楼主有吗。先谢谢了

沙颖凯

quentin011099 发表于 2014-7-16 20:56
如果该1*1地形数据能够用于GrADS作图就好了，请问楼主有吗。先谢谢了

想用的话把数据读入MATLAB然后二进制输出就可以了，不过我看GrADS应该也有自己的方法，比如http://bbs.06climate.com/forum.php?mod=viewthread&tid=22506

风儿飘飘

matlab绘制的图形还是比较流畅的

quentin011099

沙颖凯 发表于 2014-7-17 10:30
想用的话把数据读入MATLAB然后二进制输出就可以了，不过我看GrADS应该也有自己的方法，比如http://bbs.06 ...

诶，多谢楼主，我去捣鼓捣鼓