日志
matlab处理自动站站点风场资料并给出比例尺
|
grads可以处理站点资料，但这里用matlab进行处理，哪个更方便些还请大家考量。比例尺是根据沙颖凯的方法来做的（函数Get_Autoscale.m，详情看连接http://bbs.06climate.com/forum.php?mod=viewthread&tid=21699&extra=&page=1），u，v风场根据风向和风速可以得到，即WS2UV.m函数可以得到风场，Obv_highwind_SCATTER.m是最终用来得到站点风场图的程序。
clc;
clear all;
% fileID1 = fopen('D:\gzxjun3\HVDOCUMENTS\vortexline\10分钟平均风速\整点前10分钟平均风速\201006192000.txt','r');
% fileID2 = fopen('D:\gzxjun3\HVDOCUMENTS\vortexline\10分钟平均风速\每10分钟平均风速\201006192030.txt','r');
% fileID3 = fopen('D:\gzxjun3\HVDOCUMENTS\vortexline\10分钟平均风速\整点前10分钟平均风速\201006192100.txt','r');
% fileID4 = fopen('D:\gzxjun3\HVDOCUMENTS\vortexline\10分钟平均风速\每10分钟平均风速\201006192130.txt','r');
% fileID5 = fopen('D:\gzxjun3\HVDOCUMENTS\vortexline\10分钟平均风速\整点前10分钟平均风速\201006192200.txt','r');
% fileID6 = fopen('D:\gzxjun3\HVDOCUMENTS\vortexline\10分钟平均风速\每10分钟平均风速\201006192230.txt','r');
% fileID7 = fopen('D:\gzxjun3\HVDOCUMENTS\vortexline\10分钟平均风速\整点前10分钟平均风速\201006192300.txt','r');
% fileID8 = fopen('D:\gzxjun3\HVDOCUMENTS\vortexline\10分钟平均风速\每10分钟平均风速\201006192330.txt','r');
%
% C1=textscan(fileID1,'%s %s %s %s','delimiter', ',','emptyvalue', NaN);
% C2=textscan(fileID2,'%s %s %s %s','delimiter', ',','emptyvalue', NaN);
% C3=textscan(fileID3,'%s %s %s %s','delimiter', ',','emptyvalue', NaN);
% C4=textscan(fileID4,'%s %s %s %s','delimiter', ',','emptyvalue', NaN);
% C5=textscan(fileID5,'%s %s %s %s','delimiter', ',','emptyvalue', NaN);
% C6=textscan(fileID6,'%s %s %s %s','delimiter', ',','emptyvalue', NaN);
% C7=textscan(fileID7,'%s %s %s %s','delimiter', ',','emptyvalue', NaN);
% C8=textscan(fileID8,'%s %s %s %s','delimiter', ',','emptyvalue', NaN);
%
% data1 = C1{3};
% data2 = C2{3};
% data3 = C3{3};
% data4 = C4{3};
% data5 = C5{3};
% data6 = C6{3};
% data7 = C7{3};
% data8 = C8{3};
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%上面的方法太笨拙，下面的通过uigetfile使程序简练不少
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);
