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

zxj841025 · 发表于 2016-10-6 16:29:13

登录后查看更多精彩内容~

您需要登录才可以下载或查看，没有帐号？立即注册

x

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 上传

点击文件名下载附件

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

复制代码

二爷名声在外 · 发表于 2016-10-6 18:54:59

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

Andromedaaaaaa · 发表于 2017-2-23 20:17:48

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

albert007zh · 发表于 2017-11-22 18:49:42

谢谢楼主分享

xuezhanbin · 发表于 2018-3-26 20:20:15

感谢楼主分享！！！

zxj841025 · 发表于 2018-3-27 09:09:14

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

myeye · 发表于 2018-6-29 13:05:22

感谢分享 ~

		自动登录	找回密码
密码			立即注册

[源程序] matlab处理自动站站点风场资料并给出比例尺

登录后查看更多精彩内容~