NCEP风场求散度和涡度

wlzhongouc

闲着，继续玩，Matlab中有直接求散度和涡度的函数，散度命令：divergence，涡度命令：curl，以下为用NCEP风场数据做的北极上空的散度场和涡度场图：

上传图片的格式限制，就传了jpg的，看起来不是很好，有问题欢迎指出：）
代码见附件中 NCEP_wind.m (252.79 KB, 下载次数: 462)

2012-5-8 19:41 上传

点击文件名下载附件
代码
阅读权限: 20 ，希望大家多交流多贡献哦

补充一下：图中靠近0度经线的地方缺了一小块，不知道是不是NCEP数据分辨率的问题，用的双线性插值，希望听听大家的意见
更新：2012、5、8 今天听看了 @无径之林 指出的错误后，我修改了程序重新进行计算，原先错误的采用了经纬度网格点来进行计算，正确的应该是采用实际的距离来进行，这样量级在10e-7到10e-5左右，若是还有问题欢迎指出，不胜感激啦！
再次更新补充(2012、5、12)：为了验证计算的散度和涡度的正确性，听从了 @无径之林 的建议，在http://oceanwatch.pifsc.noaa.gov/las/servlets/index 下载了wind stress vector和wind stress curl的数据对程序计算的结果进行验证，对比结果如下，左图是网上涡度的数据，右图是采用我的方法计算的实际距离进而计算的涡度，能够看出大致趋势是一样的，由于风应力数据两个分量的缺省值不太匹配，所以出现计算结果出现很多空白区域


猛牛
！！！注意！！！！！！！！！！！！！！！！！！！！！！！！！！！！！！！！！！！！！！！！！！！！

补充说明：这种方法计算的散度和涡度在赤道附近准确度一般，在极地区域相当不准，没有很好地考虑到地球在极区曲率的变化，最好的办法是采用一个矩形网格，然后将矢量数据转化到矩形网格中，然后在矩形网格中计算出散度和涡度，画图的时候再对应到原来的经纬度上

补充时间：2013年8月17号 by 猛牛

声明：看到后来还有很多人来留言，我这里再更新一下说明吧，原帖子中提供的算法是不准确的。。。请查找一下在球面坐标系下的涡度和散度公式，用那个公式计算

补充时间：2016年11月2号 by 蒙牛

或者参考以下帖子，链接：
http://bbs.06climate.com/forum.php?mod=viewthread&tid=20665

ShaG

哈哈又见楼主，不过我觉气象上计算涡度和散度的时候是不可以用MATLAB自带的函数的，应该要把地球的曲率也考虑进去

% calculate curvature vorticity
R=6.3781e6; % earth's radius
[dx, ~]=gradient(longitude);
dx=dx*pi./180;
[~, dy]=gradient(latitude);
dy=dy*pi./180;
dy(1, :)=NaN;  
dy(73, :)=NaN; % 由于使用差分的缘故，把dy的头尾元素设定为缺测
[dv, ~]=gradient(v);
[~, du]=gradient(u.*cos(latitude*pi./180));
curlz=(dv./dx-du./dy)./(R.*cos(latitude*pi./180));

无径之林

wlzhongouc 发表于 2012-5-8 10:17

                               
登录/注册后可看大图

该问题后来解决了  呵呵  您说得很对  采用pcolor命令时，由于数据精度的问题，极地投影的时候总会缺一小 ...

很不错！你的程序我没看，但你的量级肯定不对，一般是1e-8 - 1e-5附近吧，具体记不清啦！
估计你是直接用lon lat做的，实际上应该用距离，curl等函数一样！
如下：
%%-------calculate global wind curl using tau data;
%------- 0.2 lat*1/3 lon----------
[m,n]=size(ua);
R0=6371004;
pi=3.1415926;
  %---calculate deltx delty dx dy,and x y--units m actual location
  deltx(1:m)=2*pi*R0*cosd(lat)/360*mean(diff(lon));  % x lon
  delty(1:n)=2*pi*R0/360*mean(diff(lat));            % y lat
  clear pi R0
  for i=1:m
      dx(i)=sum(deltx(1:i),2);
  end
  clear i m deltx
  for i=1:n
      dy(i)=sum(delty(1:i),2);
  end
  clear i n delty
  
  %y=repmat(deltx',[1,n]);
  %x=repmat(delty,[m,1]);
  [x,y]=meshgrid(dy,dx);
  [w,cav]=curl(x,y,ua,va);
  [sn,cav]=curl(x,y,va,ua);
  ss = divergence(x,y,va,ua);
  W(:,:) = (sn).^2 + (ss).^2 - w.^2 ;
  clear sn ss w cav dx dy x y

zane

请问诸位大侠，小生计算到这里，如何把散度场绘制出来？？
x=-5:0.1:5;
>> y=x;
>> u=x.*x;
>> v=-y.*y;
>> [xx,yy]=meshgrid(x,y);
>> [uu,vv]=meshgrid(u,v);
>> div=divergence(xx,yy,uu,vv);

谢过~~

子达

感谢楼主的分享

随缘

图看起来不错哦，学习了，O(∩_∩)O谢谢楼主分享

随缘

原来是猛牛作品啊，不错啊

戴戴

matlab怎么读 nc文件呀。我搞了几次 没搞好

topmad

图片很漂亮，支持一下

mofangbao

这个设置真有必要，现在只顾着索取的人太多了，而且提问之前自己不做足功课，看来论坛在某些方面需要改进了。再次感谢楼主的分享！

mofangbao

戴戴 发表于 2012-4-23 21:36

                               
登录/注册后可看大图

matlab怎么读 nc文件呀。我搞了几次 没搞好

那就继续搞吧...现有的教程加上自己的思考应该能解决的

wlzhongouc

随缘 发表于 2012-4-23 20:57

                               
登录/注册后可看大图

原来是猛牛作品啊，不错啊

图中靠近0度经线的地方缺了一小块 可能是数据分辨率的问题 也没找到办法解决。。

wlzhongouc

戴戴 发表于 2012-4-23 21:36

                               
登录/注册后可看大图

matlab怎么读 nc文件呀。我搞了几次 没搞好

在网上一搜吧 很多关于这方面的内容的 不同版本用的NC插件还不太一样

wlzhongouc

mofangbao 发表于 2012-4-23 22:11

                               
登录/注册后可看大图

这个设置真有必要，现在只顾着索取的人太多了，而且提问之前自己不做足功课，看来论坛在某些方面需要改进了 ...

呵呵 清风不用客气 论坛你贡献最大啊