- 积分
- 205
- 贡献
-
- 精华
- 在线时间
- 小时
- 注册时间
- 2014-8-21
- 最后登录
- 1970-1-1
|
发表于 2015-1-25 12:47:27
|
显示全部楼层
用MATLAB画风向杆，全免费啦！
之前自己写了一个画风旗的程序，后来在网上发现了国外某个大侠的程序，我略加修改一下，贴出来给大家参考
下面的绿色代码保存为windbarbm.m
function [] = windbarbm(ha,lat,lon,u,v,varargin)
%WINDBARBM Project wind barbs onto map axes
%
% WINDBARBM(lat,lon,u,v) projects two dimensional wind barbs onto the
% current map axes. The vector components (u,v) are in units of knots and
% are specified at the points (lat,lon). It handles winds up to 130 knots.
% Winds exceeding 130 knots will appear as 130 knots.
%
% WINDBARBM(lat,lon,u,v,s) uses the input s to scale the vectors after
% they have been automatically scaled to fit within the grid. If omitted,
% s = 0.9 is assumed.
%
% WINDBARBM(lat,lon,u,v,'PropertyName',PropertyValue,...) and
% WINDBARBM(lat,lon,u,v,s,'PropertyName',PropertyValue,...) uses the
% windbarbm object properties specified to display the windbarb objects.
% The properties supported by windbarbm are the same as the properties
% supported by linem.
%
%
% MFILE: windbarbm.m
% MATLAB: 7.8.0 (R2009a)
% VERSION: 1.3 (28 November 2011)
% AUTHOR: Nick Siler
% CONTACT: siler@atmos.washington.edu
%Argument tests (from quiverm.m)
%==============================================
axis(ha);
%==============================================
% return;
if any([ndims(lat) ndims(lon) ...
ndims(u) ndims(v) ] > 2)
error(['map:' mfilename ':inputContainsPages'], ...
'Input data can not contain pages.')
elseif length(lat) == 1 && size(lat,1) ~= size(u,1)
error(['map:' mfilename ':invalidLat'], ...
'Lat vector input must have row dimension of u.')
elseif length(lon) == 1 && size(lon,1) ~= size(u,2)
error(['map:' mfilename ':invalidLon'], ...
'Lon vector input must have column dimension of u.')
elseif ~isequal(size(lat),size(lon),size(u),size(v))
error(['map:' mfilename ':inconsistentDims'], ...
'Inconsistent dimensions for inputs.')
end
%check for scale and wind barb property specification
wbproperties = '''color'',''b'''; %default wind barb color is blue.
switch length(varargin)
case 1
if ischar(varargin{1})
error(['map:' mfilename ':invalidScale'], ...
'Invalid scale factor.')
end
scale = varargin{1};
case 0
scale = .9;
otherwise
%for an odd number of arguments, the first will be the scale factor
if rem(length(varargin),2)==1
if ischar(varargin{1})
error(['map:' mfilename ':invalidScale'], ...
'Invalid scale factor.')
end
scale = varargin{1};
nn = 2;
else
% for an even number of arguments, no scale factor is specified
scale = .9;
nn = 1;
end
for ii = nn:length(varargin)
if ischar(varargin{ii})
wbproperties = [wbproperties,',''',varargin{ii},''''];
else
wbproperties = [wbproperties,',',num2str(varargin{ii})];
end
end
end
umag = sqrt(u.^2+v.^2); %wind speed
%find theta; add pi to atan(v/u) when u<0
dummy = (u<0)*pi;
theta = atan(v./u)+dummy;
[a,b] = size(umag);
%create 18 logical matrices for 18 possible barbs. Non-zero when the barb
%is called for at that gridpoint.
g1 = umag > 7.5 & umag <= 47.5;
g2 = umag > 17.5 & umag <= 47.5;
g3 = umag > 27.5;
g4 = (umag > 37.5 & umag <= 47.5) | (umag > 57.5 & umag <= 97.5);
g5 = umag > 67.5;
g6 = (umag > 77.5 & umag < 97.5) | umag > 107.5;
g7 = umag > 87.5 & umag < 97.5 | umag > 117.5;
g8 = umag > 127.5;
g9 = (umag > 2.5 & umag <= 7.5) | (umag > 12.5 & umag <= 17.5);
g10 = umag > 22.5 & umag <= 27.5;
g11 = (umag > 32.5 & umag <= 37.5) | (umag > 52.5 & umag <= 57.5);
g12 = (umag > 42.5 & umag <= 47.5) | (umag > 62.5 & umag <= 67.5);
g13 = (umag > 72.5 & umag <= 77.5) | (umag > 102.5 & umag <= 107.5);
g14 = (umag > 82.5 & umag <= 87.5) | (umag > 112.5 & umag <= 117.5);
g15 = (umag > 92.5 & umag <= 97.5) | (umag > 122.5 & umag <= 127.5);
g16 = umag > 47.5;
g17 = umag > 97.5;
g18 = true(a,b);
%position of each barb relative to grid point: [x0 y0; x1 y1]
c1 = [-1 0;-1.125 .325];
c2 = [-.875 0; -1 .325];
c3 = [-.75 0; -.875 .325];
c4 = [-.625 0; -.75 .325];
c5 = [-.5 0; -.625 .325];
c6 = [-.375 0; -.5 .325];
c7 = [-.25 0; -.375 .325];
c8 = [-.125 0; -.25 .325];
c9 = [-.875 0; -.9375 .1625];
c10 = [-.75 0; -.8125 .1625];
c11 = [-.625 0; -.6875 .1625];
c12 = [-.5 0; -.5625 .1625];
c13 = [-.3750 0; -.4375 .1625];
c14 = [-.25 0; -.3125 .1625];
c15 = [-.125 0; -.1875 .1625];
c16 = [-1 0; -.875 .325];
c17 = [-.75 0; -.625 .325];
c18 = [0 0; -1 0];
%set scale based on average latitude spacing
[m,n]=size(lat);
scale2 = scale*(max(max(lon))-min(min(lon)))/n;
%% fang
lonlim=getm(ha,'maplonlimit');
scale2 = scale*(lonlim(2)-lonlim(1))/n;
%draw the barbs
for nn = 1:18
eval(['dummy = reshape(g',int2str(nn),',1,a*b);']);
count = sum(dummy); % number of barbs to draw
if count == 0
continue
end
%rotation operations
eval(['x1 = c',int2str(nn),'(1,1)*cos(theta)-c',int2str(nn),...
'(1,2)*sin(theta);']);
eval(['y1 = c',int2str(nn),'(1,1)*sin(theta)+c',int2str(nn),...
'(1,2)*cos(theta);']);
eval(['x2 = c',int2str(nn),'(2,1)*cos(theta)-c',int2str(nn),...
'(2,2)*sin(theta);']);
eval(['y2 = c',int2str(nn),'(2,1)*sin(theta)+c',int2str(nn),...
'(2,2)*cos(theta);']);
x1 = x1*scale2+lon;
x2 = x2*scale2+lon;
%multiply y1 and y2 by cos(lat) to compensate for the closer spacing of
%meridians.
y1 = y1*scale2.*cos(lat*pi/180)+lat;
y2 = y2*scale2.*cos(lat*pi/180)+lat;
x = [reshape(x1(dummy),1,count);reshape(x2(dummy),1,count)];
y = [reshape(y1(dummy),1,count);reshape(y2(dummy),1,count)];
eval(['linem(y,x,',wbproperties,')']);
end
下面蓝色代码是测试程序，保存为tryWindbarbm.m
wdata=load('windUpload.data','w');
latout=wdata(:,1); % 纬度
lonout=wdata(:,2); % 经度
uktout=wdata(:,3); % u分量 (KT为单位)
vktout=wdata(:,4); % v分量 (KT为单位)
figure
set(gcf,'outerposition', get(0,'screensize')); % 全屏图像
hold on
set(gca,'position',[0.02 0.02 0.96 0.96]);
axis off
%% ===============================================
% 设置地图坐标属性
% setm(hca);%查看当前可以设置的所有图形坐标轴(map axes)的属性
latlim=[25 37];
lonlim=[100 130];
hca=axesm ('MapProjection','eqdcylin',...
'MapLatLimit',latlim,...
'MapLonLimit',lonlim,...
'Frame','on','grid','on');
setm(hca,'Frame','on');%使框架可见
setm(hca,'Grid','on');%打开网格
intcint=2; % 网格间距
setm(hca,'mlinelocation',intcint);
setm(hca,'plinelocation',intcint);
lonn=getm(gca,'maplonlimit');
latt=getm(gca,'maplatlimit');
setm(hca,'MLabelLocation',[min(lonn)-intk :intcint :max(lonn)]);%标上经度刻度标签,每隔10度
setm(hca,'MeridianLabel','on');%设置经度刻度标签可见
setm(hca,'MLabelParallel','south');%将x轴标识，也就是180w， 150w 等标在-90S（南维90）上
setm(hca,'PLabelLocation',[min(latt)-intk:intcint:max(latt)]);%标上经度刻度标签,[-90:30:90]
setm(hca,'ParallelLabel','on');%设置经度刻度标签可见
setm(hca,'PLabelMeridian','west');%将纬度标签放置在经度为90度的地方
% % ==============================================
load coast;
axis(hca);
plotm(lat,long,'k-')
widthwindbar=0.03;
windbarbm(hca,latout,lonout,uktout,vktout,widthwindbar,'color','b','linewidth',1.2);
下面红色的是测试数据，保存为windUpload.data
36.3833 115.0167 54.4647 7.6545
27.9667 115.8000 35.7317 4.3873
28.7000 116.5333 40.6944 4.9966
28.9500 122.3167 48.0997 9.3496
36.7500 115.2167 58.2736 9.2296
30.1833 106.8500 30.5290 -5.3831
29.4667 106.6000 28.9823 1.0121
29.0833 106.6000 21.9464 -1.5346
28.3167 115.9333 37.0000 0
28.9833 116.1333 44.8904 -3.1390
28.3000 114.6667 33.5814 5.3188
29.4500 116.8333 42.9411 -2.2504
30.0333 117.2667 47.8831 -3.3483
30.3833 117.2833 48.8806 -3.4181
36.3500 123.3667 67.9731 -31.6964
28.1667 113.5000 30.8302 3.2404
30.7500 117.3000 48.8806 -3.4181
31.1000 117.3167 51.9683 -1.8148
31.8167 117.3500 44.4326 11.9057
35.5667 123.3833 52.1301 -25.4255
34.8000 123.3667 50.3325 -24.5488
28.1167 112.3167 27.8935 -2.4404
36.6833 118.7667 60.9628 2.1289
36.8000 117.5667 62.9904 1.0995
32.1833 117.3333 42.7644 4.4947
32.5333 117.3167 40.9938 -0.7155
32.8833 117.3167 44.9948 -9.5639
28.0333 111.1667 25.4318 -5.4057
27.7500 110.0333 26.8521 -2.8223
25.0167 115.9000 33.8137 3.5540
25.4833 116.5667 33.9793 1.1866
36.8000 116.3333 56.6878 5.9581
25.9667 117.2167 32.9950 -0.5759
26.4500 117.8667 27.6553 4.3802
33.6000 117.2833 49.0000 0
27.0000 118.9833 26.8973 2.3532
27.3167 119.7333 32.9548 -1.7271
样图为：
|
-
|
[复制链接]
窥视卡
雷达卡
显身卡