matlab处理自动站站点风场资料并给出比例尺

zxj841025

grads可以处理站点[color=rgb(51, 51, 51) !important]资料，但这里用matlab进行处理，哪个更方便些还请大家考量。比例尺是根据沙颖凯的方法来做的（函数Get_Autoscale.m，详情看连接http://bbs.06climate.com/forum.p ... p;extra=&page=1），u，v风场根据风向和风速可以得到，即WS2UV.m函数可以得到风场，Obv_highwind_SCATTER.m是最终用来得到站点风场图的程序。 Get_Autoscale.m (293 Bytes, 下载次数: 35)

2016-10-6 16:25 上传

点击文件名下载附件 obv_highwind_vector.m (5.17 KB, 下载次数: 14)

2016-10-6 16:25 上传

点击文件名下载附件 STRCELL2NUMMATR.m (379 Bytes, 下载次数: 15)

2016-10-6 16:25 上传

点击文件名下载附件 WS2UV.m (760 Bytes, 下载次数: 34)

2016-10-6 16:25 上传

点击文件名下载附件 2010061918.txt (28.95 KB, 下载次数: 14)

2016-10-6 16:27 上传

点击文件名下载附件

1. clc;
   
2. clear all;
   
3. 
   
4. num = 1;
   
5. fn=sprintf('%2imv',33);%note:这里是%2i，所以必须是33（两位数），如果改成3（即一位数），后面的print打印会出错
   
6. [s,e]=dos(['md ' fn]);%用dos命令创建一个文件
   
7. % highwind =zeros(num,1);
   
8. for i = 1:num
   
9.     [FileName,PathName] = uigetfile('*.txt','选择需要读入的数据文件');%这里的数据为桂自动站资料,时间为20:00-23:30（LST）
   
10.     fileID=fopen(strcat(PathName,FileName),'r');%strcat 连接字符串
    
11.     C=textscan(fileID,'%s %s %s %s','delimiter', ',','emptyvalue', NaN);
    
12.     data = C{3};
    
13.    [row,col] = size(data);
    
14.    higwnd = STRCELL2NUMMATR(data);
    
15.    data1 = C{4};
    
16.    wnddir = STRCELL2NUMMATR(data1);
    
17.    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    
18.    %获取经纬度信息%
    
19.     [FileName1,PathName1] = uigetfile('*.txt','选择需要读入的数据文件');%这里的数据为桂自动站资料,时间为20:00-23:30（LST）
    
20.     fileID=fopen(strcat(PathName1,FileName1),'r');%strcat 连接字符串
    
21.     C1=textscan(fileID,'%s %s %s %s','delimiter', ',','emptyvalue', NaN);
    
22.     data2=C1{2};
    
23.     data3=C1{3};
    
24.     lon = STRCELL2NUMMATR(data2);
    
25.     lat = STRCELL2NUMMATR(data3);
    
26.    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    
27.    
    
28.    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    
29.    %根据风速，风向求解出u，v风场，并根据经纬度信息画出经纬度对应点处的离散风场
    
30.    [u,v] = WS2UV(higwnd,wnddir);
    
31.    %% quiver label
    
32. labelx=[108.8 108.85]; labely=[24.1 24.15]; % location of label % <------- Use the actual grid distance of the data
    
33. labelu=[5 0]; labelv=[0 0]; % set 8m/s as a standard
    
34. scale_auto=Get_Autoscale(lon, lat, u, v);
    
35. scale_label=Get_Autoscale(labelx, labely,  labelu, labelv);
    
36. scale_factor=scale_label/scale_auto;
    
37. Handle=quiver(labelx, labely, labelu, labelv, 1.0);
    
38. set(Handle, 'Color', 'r', 'LineWidth', 2.5);
    
39. gtext([num2str(labelu(1),'%.2d') ' m/s'],...
    
40.     'FontWeight', 'bold', 'Color', 'r', 'FontName', 'Helvetica', 'FontSize', 12);
    
41. % text(labely(1)-0.1, labelx(2), [num2str(labelu(1),'%.2d') ' m/s'],...
    
42. %     'FontWeight', 'bold', 'Color', 'r', 'FontName', 'Helvetica', 'FontSize', 12);
    
43.    hold on
    
44.    Handle = quiver(lon,lat,u,v,1.0*scale_factor,'color','k');%<----- size: 2.0*scale_factor
    
45.    set(Handle, 'Color', [0.5 0.5 0.5], 'LineWidth', 1.5)
    
46.    
    
47.    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    
48.    %添加比例尺度，目的是让多附图有可比性
    
49. 
    
50.    guojie=shaperead('F:\usefull_skills\matlab-画中国地图\bou2_4l.shp');%读取国界shp文件的内容
    
51. %     guojie=shaperead('F:\bou2_4l.shp');%读取国界shp文件的内容
    
52.     bou1_4lx=[guojie(:).X];%提取经度信息
    
53.     bou1_4ly=[guojie(:).Y];%提取纬度信息
    
54.     plot(bou1_4lx,bou1_4ly,'-k','LineWidth',1.2)%绘国界
    
55. %     set(gca,'fontsize',19,'fontweight','bold','xlim',[105.958558576 112.01936941],'ylim',[22.952702714 27.045045072],'xtick',[106:1:112],'xticKlabel',{'106E','107E','108E',...
    
56. %     '109E','110E','111E','112E'},'ytick',[23:0.5:27],'YTICKLABEL',{'23N','','24N','','25N','','26N','','27N'},'linewidth',3);
    
57.     axis([min(lon(:)) max(lon(:)) min(lat(:)) max(lat(:))]);%设置显示的经纬度范围
    
58.     xlabel('Longitude','fontweight','bold','fontsize',25);
    
59.     ylabel('Latitude','fontweight','bold','fontsize',25);
    
60.     set(gca,'fontsize',15)
    
61.    % axes('position',[0.1 0.1 0.8 0.8])%axes是figure的子对象，定义一个画图区域
    
62. %    print ('-djpeg', 'D:\gzxjun3\HVDOCUMENTS\vortexline\10分钟平均风速\FileName.jpg')
    
63. hold off
    
64. %     print(gcf,'-djpeg',sprintf('%s\\%10d',fn,str2double('FileName(1:end-4)')));
    
65.     print(gcf,'-djpeg',sprintf('%s\\%10s',fn,FileName(1:end-4)));
    
66.    fclose(fileID);
    
67. end
    
68. % plot(highwind/10,'LineStyle','-','Marker','o','color','k','linewidth',4);
    
69. % %除以10是因为资料里的风速扩大了10倍
    
70. % % set(gca,'xlim',[1,8],'xtick',(1:1:8),'xticklabel',{'1700','','1900','','2100','','2300',''},'fontsize',15);
    
71. % set(gca,'xlim',[1,5],'xtick',(1:1:5),'xticklabel',{'1700','','1900','','2100'},'fontweight','bold','fontsize',25);
    
72. %
    
73. % legend('Wind speed','fontsize',25);
    

复制代码

二爷名声在外

个人感觉这种箭头显得有点杂乱
可能比较习惯windbarbm那种风向标

Andromedaaaaaa

如何选择的标准风速呢？为什么选作8m/s呢？

albert007zh

谢谢楼主分享

xuezhanbin

感谢楼主分享！！！

zxj841025

哈哈 没想到薛师弟你也看到这个帖子了

myeye

感谢分享  ~