Hourly variation of wind for each month, using Quiver

I have three vectors: date (TIEMPO), Wind speed (VV) and Wind direction (DV).

I would like to combine Contour Plot and Quiver to display the anual variation of wind, yet, given that's a single point instead of a map, I don't know how to use "quiver" function without coordinates.

Here's my code:

 %%TEST
 clear all
 load Viento_DMC

 %separates time in diferent variables (year, month, day and hour)
 [ANHO, MES, DIA, HORA, mn, s] = datevec(TIEMPO);
 %delete useless variables
 clearvars mn s fecha
 %creates a datenum based only on year and month
 TIEMPO2=datenum(ANHO,MES,1,0,0,0);
 %gives a single date for each month of each year (in order to get how many
 %months are in the dataset)
 TIEMPO3=unique(TIEMPO2);
 %max number of months in datafile
 meses_max=numel(TIEMPO3);

 %creates V and U components
 U=VV.*sind(DV);
 V=VV.*cosd(DV);

 %Creates a daly cicle for wind speed (mean of every hour of each month) 
 cda=grpstats(VV, {hour(TIEMPO), TIEMPO2},'mean'); 
 %creates a matrix for daly cicle of wind speed (24 hours x 60 months)
 cda2d=reshape(cda, [meses_max,24]); %24 hours

 %Creates a daly cicle U and V components 
 cdaU=grpstats(U, {hour(TIEMPO), TIEMPO2},'mean');
 cdaV=grpstats(V, {hour(TIEMPO), TIEMPO2},'mean');
 cdaDir=atan2d(cdaU,cdaV); 
 %creates a matrix for U and V components
 cda2dDir=reshape(cdaDir, [meses_max,24]);

 % create x/y grid 24x60 
 [x,y] = meshgrid(1:24,1:60) ;
 % plot Wind speed and wind direction matrix's on xy grid with quiver
  quiver( x,y,cda2d,cda2dDir)

This is how it's supposed to look:

The data it's attached. I would really apreciate any help you could give me. Thanks!

[jonas]

Undefined function 'hour' for input arguments of type 'double'.

Error in Prueba_Quiver (line 17)
cda=grpstats(vel, {hour(TIEMPO), mesunico},'mean'); %wind speed

Here you're using the duration function hour on a datevec? It's difficult to follow your coding so please fix any unnecessary bugs in the script so that we can focus on the actual issue, which is about plotting.

[Emanuel Valdes]

I just update the code, trying to be as clear as possible. Regarding the function hour, it's used on a datenum, not a datavec.

[jonas]

My point was that you cannot (as far as I know) use the function hour on a double (datevec/datenum, doesn't matter). You need a datetime format. I'm still getting the same error. Can you run the code yourself?

Answer 1

[jonas]

This is my attempt with normalized/scaled arrows. See attachment.

load Viento_DMC
T=datetime(datevec(TIEMPO));

%%Velocity and direction
U=VV.*sind(DV);
V=VV.*cosd(DV);
Dir=atan2d(U,V);
TT=timetable(T,VV,Dir,U,V);
%%Take hourly mean (probably not necessary)
TT2=retime(TT,'hourly','mean');

%%Calculate hourly average per unique month/year (not sure if this is what you want, but easy to adjust)
[~,~,G]=unique([year(TT2.T) month(TT2.T) hour(TT2.T)],'rows');
VVm = splitapply(@nanmean,TT2.VV,G);
Um = splitapply(@nanmean,TT2.U,G);
Vm = splitapply(@nanmean,TT2.V,G);
VVm=reshape(VVm,24,numel(VVm)/24);

%%Normalize Um and Vm (UNCOMMENT TO NORMALIZE VECTORS)
for i=1:length(Um)
    Umn(i)=Um(i)./norm([Um(i) Vm(i)]);
    Vmn(i)=Vm(i)./norm([Um(i) Vm(i)]);
end
Um=Umn;
Vm=Vmn;
Um=reshape(Um,24,numel(Um)/24);
Vm=reshape(Vm,24,numel(Vm)/24);

%%create x/y grid 24x60
[x,y] = meshgrid(linspace(0,1,24),linspace(0,1,60)) ;

%%plot Wind speed and wind direction matrix's on xy grid with quiver
axis([0 1 0 1])
contourf(x,y,VVm');hold on
quiver(x,y,Um',Vm',0.5,'color','k')
axis([0 1 0.27 1])

%%Manipulate Labels
nLabelsY=10; %(60/nLabelsY must be integer NOT float)
[~,id,~]=unique([year(TT2.T) month(TT2.T)],'rows');
yticklabs=cellstr(datestr(TT2.T(id),'yyyy-mmm'))
set(gca,'ytick',linspace(0,1,nLabelsY),'yticklabels',{yticklabs{1:60/10:end}})
set(gca,'xtick',linspace(0,1,24),'xticklabels',sprintfc('%g',1:24))

[Chad Greene]

Nice!

[Emanuel Valdes]

Excellent!!! Thanks a lot!!!!

[jonas]

Thanks, happy to help. I will fix the labels and call it a day.

[jonas]

@Emanual, I found a quite serious error with the code. From your original code, I figured that the input for quiver was velocity/angle. However, after doing some research, I think that it's actually the U/V components that you want to throw as input. If you do so, then you get this figure:

This is how a quiver plot normally looks like, where the arrows are proportional to the length of the vectors (code attached). If you want the same length on all vectors, then you can normalize, however the plot turns out quite ugly. This is because the velocity is almost zero at many places (the blue areas), so the direction is kind of random. I've attached the .m-file to this comment. Uncomment the relevant section if you want to normalize.

I've also decided to plot on a normalized axes scale. I believe the intended use of quiver is to plot vectors on a spatial grid, not direction over time. That's why it's a bit tricky to get it right. If you have different scales on the axes, then you have to re-scale all the quivers.

[Emanuel Valdes]

@Jonas, I have a last inquery... the plot you made is showing from where the wind is coming from, which is mathematically accurate, but meteorologically wrong. It should show where the wind is going, in other words, the arrows should go in the other direction. How could I fix that?

[jonas]

If my geometry is not all wrong then you just multiply both U and V by -1

[Emanuel Valdes]

I guess it's conceptually more accurate to add 180 to DV. It's the same result any way

Answer 2

[Chad Greene]

A few points:

1. Instead of this:

[ANHO, MES, DIA, HORA, mn, s] = datevec(TIEMPO);
%delete useless variables
clearvars mn s fecha

You can simply do

[~, ~, ~, HORA, ~, ~] = datevec(TIEMPO);

which never writes the unnecessary variables.

Here's the closest I can get to what you want. It uses my xyz2grid function from File Exchange.

load Viento_DMC

%separates time in diferent variables (year, month, day and hour)
[~, ~, ~, HORA, ~, ~] = datevec(TIEMPO);

%creates V and U components
U=VV.*sind(DV);
V=VV.*cosd(DV);

% Grid by day and hour: (UG and VG indicate "gridded")
[H,D,UG] = xyz2grid(HORA,round(TIEMPO),U); 
VG = xyz2grid(HORA,round(TIEMPO),V); 

% Get 31 day moving mean by operating down the rows:
UG = movmean(UG,31,1,'omitnan'); 
VG = movmean(VG,31,1,'omitnan'); 

figure
pcolor(H,D,hypot(UG,VG))
shading interp

cb = colorbar; 
ylabel(cb,'velocidad media m s^{-1}')
caxis([0 4]) 

hold on
sz = [48 24]; 
Hr = imresize(H,sz); 
Dr = imresize(D,sz); 

% Scale by speed so all vectors are same length.
UGr = imresize(UG./hypot(UG,VG),sz); 
VGr = imresize(VG./hypot(UG,VG),sz); 

da = daspect; 
quiver(Hr,Dr,UGr,VGr*da(2)/da(1),0.2,'k')

datetick('y','mmmyy','keeplimits')

I'm not sure how to make the arrows prettier. I've tried to account for the data aspect ratio (because the units of x and y axes are not the same), but the arrows are still a bit squashed. Perhaps one of the alternate quiver functions on File Exchange will solve the problem for you.

Answer 3

[Emanuel Valdes]

Guys, I think I solved it. The problem was that I didn't make a matrix for U and V components, I was trying to plot speed and direction.

Here's the code fixed:

 %%TEST
 clear all
 load Viento_DMC

 %separates time in diferent variables (year, month, day and hour)
 [ANHO, MES, DIA, HORA, mn, s] = datevec(TIEMPO);
 %delete useless variables
 clearvars mn s fecha
 %creates a datenum based only in year and month
 TIEMPO2=datenum(ANHO,MES,1,0,0,0);
 %gives a single date for each month of each year (in order to get how many
 %months are in the dataset)
 TIEMPO3=unique(TIEMPO2);
 %max number of months in datafile
 meses_max=numel(TIEMPO3);
 mesunico=ANHO*100+MES;
 meses=datestr(TIEMPO3,'mmm-yy','local');
 meses_max=numel(TIEMPO3);

 %creates V and U components
 U=VV.*sind(DV);
 V=VV.*cosd(DV);

 %Creates a daly cicle for wind speed (mean of every hour of each month) 
 cda=grpstats(VV, {hour(TIEMPO), TIEMPO2},'mean'); 
 %creates a matrix for daly cicle of wind speed (24 hours x 60 months)
 cda2d=reshape(cda, [meses_max,24]); %24 hours  

 %Creates a daly cicle U and V components 
 cdaU=grpstats(U, {hour(TIEMPO), TIEMPO2},'mean');
 cdaV=grpstats(V, {hour(TIEMPO), TIEMPO2},'mean');
 %create a matrix for U and V 
 cdaU2=reshape(cdaU, [meses_max,24]);
 cdaV2=reshape(cdaV, [meses_max,24]);
 % create x/y grid 24x60 
 [x,y] = meshgrid(1:24,1:meses_max) ;
 % put it on quiver
 quiver(x,y,cdaU2,cdaV2) 

 %formato del gráfico
 horas=0:1:23;
 set(gca,'xtick',1:24,'xticklabel',horas);
 set(gca,'ytick',1:meses_max,'yticklabel',meses); 

 set(gcf, 'Position', [100, 100, 600, 900]);
 ylabel('Mes')
 xlabel('Hora')
 colorbar; 
 colormap jet

[jonas]

OK! Happy it worked out for you. I'm still getting the same error.

Undefined function 'hour' for input arguments of type 'double'.

Error in Prueba_Quiver (line 25) cda=grpstats(VV, {hour(TIEMPO), TIEMPO2},'mean');

From the doc: h = hour(t) returns the hour numbers of the datetime values in t. The h output is a double array the same size as t and contains integer values from 0 to 23.