Determining months of each year

I have calculated the average rainfall for each month from 2005-2015 and have obtained a graph.However the graph doesn't differentiate which year the monthly average belong to? Is there anyway for me to label the months according to marker colours?

this is the code i used, and the graph generated

<<

>>

%% Site 525 rainfall data
% Simple timeseries plot

% Load the data
R= load ('../pdata/MIDAS_srcid00525_daily_rainfall.txt');
rain = (R(:,3)); % raw (not interpolated/filled)
rain(rain<0) = NaN; % identify the missing data and set to NaN
% Datenum: point in time as the number of days from January 0, 0000
dates = datenum(R(:,2));
% subset data (2005-2015)
sub_idx = find(dates>=datenum(2005,1,1) & dates<=datenum(2015,12,31));
dates = dates(sub_idx); %dates from 2005-2015
rainfall = rain(sub_idx); %rainfall data for 2005-2015
%% Monthly Average 
[Y,M] = datevec(dates);
%
[a,~,c] = unique([Y,M],'rows');
%
rain_data = [a, accumarray(c,rainfall,[],@(x)mean(x))];
rain_ave = (rain_data(:,3));
date_plot = (rain_data(:,2));

%% Make a figure and plot for max min
figure('units','normalized','outerposition',[0 0 1 1]);
% Set the font size for the figure
font_size = 22;
set(0, 'DefaultAxesFontSize', font_size);
% Set the figure background to be white
set(gcf,'color','w');

% Plot data

p1 = plot(date_plot,rain_ave, '+');
p1.LineWidth = 1.2;
p1.MarkerSize = 10;
p1.Color = [0 0 1];

dynamicDateTicks;

ylabel('Total rainfall (mm)');
xlabel('Year (2005-2015)')
grid on
box on
xlim([min(date_plot) max(date_plot)]);

Answer 1

[Akira Agata]

I would recommend converting your data set to 'timetable', and applying retime function. Assuming dates and rainfall are datenum and numerical vector with the same size, say N-by-1, the following code can calculate monthly average and plot the result.

dates = datetime(dates,'ConvertFrom','datenum');
rainfall = rand(size(dates));
TT = timetable(dates,rainfall);
TT = retime(TT,'monthly','mean');

% Plot monthly average
plot(TT.dates,TT.rainfall,'o');