Plot sunrise and sunset times in MATLAB or Octave

Winter is coming and so is summer (in the other hemisphere of the world). The length of the days are changing and combined with daylight saving time, you might want to know just when you get to see the sun again. This tutorial shows you how to make a graph that plots the sunrise and sunset times in MATLAB or Octave. You'll learn how to shade the area above, below and between curves, use date and time formatted axis tick labels and draw grid lines for each tick mark in a custom colour. The result will look like this:

This code works in Octave but it should be compatible with MATLAB, which I don't have access to at the moment.

For Octave users, fltk will mess up the plots in a lot of ways, so make sure you set the graphics toolkit to gnuplot:

% for Octave users
graphics_toolkit("gnuplot")

You can get sunrise and sunset data in a variety of ways:

• suncycle.m: computes the sunrise and sunset times (as well as solar altitude and radiation) for a set of dates. The times are reported in GMT and you will have to convert them to local time for the location you want. You also need to adjust for daylight saving time. Finally, you will need to convert the decimal hours that suncycle outputs into decimal days, i.e. for a sunrise time of 11:38AM, suncyle outputs 11.633 hours = 11.633 / 24 = 0.485 days.
• Sunrise, sunset, daylight in a graph: a website that allows you to save the times in a CSV file. You will need to parse the times and convert them to decimal days, i.e. 11:38 AM = (11 + 38 / 60) / 24 = 0.485 days.

If you are making the graphs for polar regions, you also need to take into account polar night and polar days (no sunrise or sunset) and adjust the code for nice graphs.

I don't want to complicate this tutorial so I will create some fake sunrise and sunset data for a non-polar city (Eutopia, United Democratic Republics of Freedomness):

clear all
 
% create a fake sunrise / sunset curve using cos.
year = 2099;
dates = [datenum(2099,1,1):datenum(2099,12,31)]';
scale = 2*pi / length(dates);
sunset_times = 2 * (cos(scale*(dates - floor(median(dates)))) + 1) + 17;
sunrise_times = -2 * (cos(scale*(dates - floor(median(dates)))) + 1) + 8;
daylength_times = sunset_times - sunrise_times;
% DST adjustments
dst_start_idx = find(dates == datenum(2099,3,10));
dst_end_idx = find(dates == datenum(2099,11,3));
sunset_times(dst_start_idx:dst_end_idx,:) = sunset_times(dst_start_idx:dst_end_idx,:) + 1;
sunrise_times(dst_start_idx:dst_end_idx,:) = sunrise_times(dst_start_idx:dst_end_idx,:) + 1;
city = 'Eutopia, UDRF';

The time from the code above and suncycle.m are both in decimal hours but we need the sunset and sunrise times to be in decimal days (or a valid date number) to make use of MATLAB's date formatting capabilities.

MATLAB represents dates and times with a date number, which is a decimal number of days since some internal reference date (January 0, 0000). To convert the hours to a date number, take the decimal hours [0,24] and normalize it to decimal days [0,1]:

% Convert the times to datenum by simply converting each time to a decimal day.
sunset_times = sunset_times / 24;
sunrise_times = sunrise_times / 24;
daylength_times = daylength_times / 24;
hours = [0:24] / 24; % hours of a day as decimal days datenum.

When 0 ⇐ date number ⇐ 1, MATLAB treats it as a date Jan 0, 0000 with only hours, minutes and seconds, which is what we want to be able to plot the sunrise and sunset time in hourly format.

If you think that the converted date numbers are odd (Octave, for example, which treats date any number that is [0,1] as Dec 31 -1 rather than January 0, 0000), you can add an arbitrary date number, e.g. datenum(2000,1,1) to all of the converted times above but this is not necessary:

% optional, not necessary
sunset_times = sunset_times + datenum(2000,1,1);
sunrise_times = sunrise_times + datenum(2000,1,1);
daylength_times = daylength_times + datenum(2000,1,1);
hours = hours + datenum(2000,1,1); % hours of a day as decimal days datenum.

The explanations in the rest of the tutorial assumes you did not add the arbitrary date number and that the times are [0,1], but the code works for both ways.

We want the X-axis to show tickmarks for every month at the 1st of each month. So construct a list of date numbers corresponding to the start of each month.

% Get the datenums corresponding to the start of each month.
month_starts = datenum(year, 1:12, 1);

Now let's start to build the graph. We want to shade the area between the sunrise and sunset curves with a light blue, representing day. We also want to shade the area below the sunrise curve and the area above the sunset curve with a dark blue, representing night. We do this by using building polygons and using patch to draw them on the curve.

% create a new figure.
figure;
 
% ---- Plot the shaded areas below, between and above the curves ----
 
hold on;
% fill between sunrise and sunset: day time. use flipud so the resulting vectors describe
% a properly ordered polygon.
patch([dates;flipud(dates)],[sunrise_times;flipud(sunset_times)],[0.7 0.9 1])
% fill between top of graph ( which is hours(end) ) and sunset time: night after sunset.
patch([dates;flipud(dates)],[sunset_times;repmat(hours(end),length(dates),1)],2*[0.05 0.1 0.3])
% fill between bottom of graph (hours(1)) and sunrise time: night before sunrise.
patch([dates;flipud(dates)],[sunrise_times;repmat(hours(1),length(dates),1)],2*[0.05 0.1 0.3])

Let's look at the first patch statement, the polygon between the sunrise and sunset curves. For patch, the first argument is the X coordinate of each point of the polygon and the second is the Y coordinate. We need to construct a polygon with the sunrise and sunset curve points in the right order. Consider the following:

C-------D
|       |
A-------B

If we have:

• Point A = Sunrise Start
• Point B = Sunrise End
• Point C = Sunset Start
• Point D = Sunset End

Then, we need to construct the polygon with points in the order A-B-D-C. This is why the arguments to patch are:

patch([dates;flipud(dates)],[sunrise_times;flipud(sunset_times)],[0.7 0.9 1])

dates and sunrise_times are A-B. Then we flip dates and sunset_times to represent D-C. Since dates and the different time arrays are column vectors, we use flipud to flip it up and down.

The other 2 patches (time before sunrise and time after sunset) are filled similarly. For the Y values, we use the sunset time and the top of the graph (hours(end)) or the sunrise time and the bottom of the graph (hours(1)).

Set up the X and Y limits of the graph:

% set limits to Jan 1 to Dec 31.
xlim([dates(1) dates(end)])
% set limits to time between 0:00 and 24:00 as decimal days.
ylim([hours(1) hours(end)])

Treating the dates and times as simply numbers makes this pretty straightforward.

Now, set up the tickmarks with date and time formatting:

h = gca; % get current axis.
 
% Set the XTicks to the values in month_starts.
set(h,'XTick',month_starts)
% Set the YTicks to the hours of the day.
set(h,'YTick',hours)
 
% Reformat the XTicks to 3 letter month names (Jan, Feb, etc) using 'mmm' format.
% Preserve the total number of ticks by using the 'keepticks' option and
% preserve the plot limits by using 'keeplimits'
datetick(h,'x','mmm','keepticks','keeplimits')
 
% Reformat YTicks to 24 hours format 'HH:MM'
datetick(h,'y','HH:MM','keepticks','keeplimits')

Here, we've set the X ticks to the date number values for the 1st of each month and the Y ticks to the date number values for each hour from 0:00 to 24:00 (which are decimal days, so [0,1]). Then, datetick formats the ticks to an actual date or time format, which you can specify using a format string(mmm or HH:MM).

We want vertical grid lines for each month tick and horizontal grid lines for each hour tick. Since we want to set the colour of the lines, we have to manually create each line to plot.

Create the vertical lines:

% vertical grid lines: one at each tick mark on the x-axis.
% Get the limits of the y axis, which are datenums representing [0:00,24:00]
ylimits= get(h,'Ylim');
% for each tick mark line, we have two points: (tickmark, 0:00) and (tickmark, 24:00).
% create X = [Jan-1st Feb-1st ... Dec-1st; Jan-1st Feb-1st ... Dec-1st] (as datenums)
X = repmat(month_starts,2,1); % x-values
% create Y = [0:00 0:00 ... 0:00; 24:00 24:00 ... 24:00] (as datenums)
Y = repmat(ylimits',1,size(month_starts,2)); % y-values.
% Plot the grid lines with gray colour.
plot(X,Y,'color',[0.5 0.5 0.5])

The above is pretty straightforward. X and Y represent the multiple grid lines, with each column representing one line. Each vertical grid line consists of two points, for example: (Jan 1st, 0:00) and (Jan 1st, 24:00), where both numbers are date numbers representing the date and time.

Create the horizontal grid lines:

% horizontal grid lines: one for each hour on the y-axis.
% Get limits of the x axis, which are datenums representing [Jan-1st,Dec-31st]
xlimits= get(h,'xlim');
% for each hour line, we have two points (as datenums): (Jan-1st, Hour) and (Dec-31st, Hour).
% create X = [Jan-1st Jan-1st ... Jan-1st; Dec-31st Dec-31st ... Dec-31st; ]
X2 = repmat(xlimits',1,size(hours,2)); % x-values
% create Y = [0:00 1:00 ... 24:00; 0:00 1:00 ... 24:00] (as datenums)
Y2 = repmat((hours),2,1); % y-values
% Plot the grid lines with gray colour.
plot(X2,Y2,'color',[0.5 0.5 0.5])

Again, X and Y represent the multiple grid lines, with each column representing one line. Each horizontal grid line consists of two points, for example: (Jan 1st, 0:00) and (Dec 31st, 0:00), where both numbers are date numbers representing the date and time.

Plot the sunrise time, sunset time and daylight length curves in a nice colour:

% finally, plot the sunset, sunrise and daylength times once so we get lines for them.
plot(dates,sunrise_times,'color',[1 0.5 0],'linewidth',3)
plot(dates,sunset_times,'color',[1 0.5 0],'linewidth',3)
plot(dates,daylength_times,'color',[1 0.9 0],'linewidth',3)

Next, add text labels for each curve:

% Text labels for each curve.
 
% get the index where the min / max occurs for each variable.
% median is just used to deal with curves where there are multiple extreme values.
 
sunrise_imin = floor(median(find(min(sunrise_times) == sunrise_times)));
daylength_imax = floor(median(find(max(daylength_times) == daylength_times)));
sunset_imax = floor(median(find(max(sunset_times) == sunset_times)));
 
% For all 3 labels, use the daylength max as the horizontal position.
% This could be adjusted for the southern hemisphere, since the 
% daylength minimum is in the middle of the year.
 
% sunrise time label: use the sunrise min as the vertical position.
text(dates(daylength_imax),sunrise_times(sunrise_imin),'Sunrise',...
	'VerticalAlignment','bottom',...
	'HorizontalAlignment','center',...
	'FontSize',10)
 
% sunset time label: use the sunset max as the vertical position.
text(dates(daylength_imax),sunset_times(sunset_imax),'Sunset',...
	'VerticalAlignment','top',...
	'HorizontalAlignment','center',...
	'FontSize',10)
 
% daylength label: use the daylength max as the verical position.
text(dates(daylength_imax),daylength_times(daylength_imax),'Daylight',...
	'VerticalAlignment','top',...
	'HorizontalAlignment','center',...
	'FontSize',10)

The text labels are positioned horizontally where the daylight hours is maximum and positioned vertically at either the minimum (for sunrise) or maximum (sunset, daylight hours) curves.

Finally add the title:

title(sprintf('Sunrise and Sunset Times\n%s (%d)', city, year))

The result, in Octave, looks like this:

I saved an SVG using Octave from sunrise and sunset data generated by suncycle.m for Toronto:

The above was exported from SVG to PNG using Inkscape and looks pretty good.

% Create a sunrise and sunset graph in MATLAB / Octave tutorial
% By Peter Yu
% http://www.peteryu.ca/tutorials/matlab/sunrise_sunset_graph
 
% for Octave users
graphics_toolkit("gnuplot")
 
clear all
 
% create a fake sunrise / sunset curve using cos.
year = 2099;
dates = [datenum(2099,1,1):datenum(2099,12,31)]';
scale = 2*pi / length(dates);
sunset_times = 2 * (cos(scale*(dates - floor(median(dates)))) + 1) + 17;
sunrise_times = -2 * (cos(scale*(dates - floor(median(dates)))) + 1) + 8;
daylength_times = sunset_times - sunrise_times;
% DST adjustments
dst_start_idx = find(dates == datenum(2099,3,10));
dst_end_idx = find(dates == datenum(2099,11,3));
sunset_times(dst_start_idx:dst_end_idx,:) = sunset_times(dst_start_idx:dst_end_idx,:) + 1;
sunrise_times(dst_start_idx:dst_end_idx,:) = sunrise_times(dst_start_idx:dst_end_idx,:) + 1;
city = 'Eutopia, UDRF';
 
% Convert the times to datenum by simply converting each time to a decimal day.
sunset_times = sunset_times / 24;
sunrise_times = sunrise_times / 24;
daylength_times = daylength_times / 24;
hours = [0:24] / 24; % hours of a day as decimal days datenum.
 
% Get the datenums corresponding to the start of each month.
month_starts = datenum(year, 1:12, 1);
 
% create a new figure.
figure;
 
% ---- Plot the shaded areas below, between and above the curves ----
 
hold on;
% fill between sunrise and sunset: day time.
patch([dates;flipud(dates)],[sunrise_times;flipud(sunset_times)],[0.7 0.9 1])
% fill between top of graph ( which is hours(end) ) and sunset time: night after sunset.
patch([dates;flipud(dates)],[sunset_times;repmat(hours(end),length(dates),1)],2*[0.05 0.1 0.3])
% fill between bottom of graph (hours(1)) and sunrise time: night before sunrise.
patch([dates;flipud(dates)],[sunrise_times;repmat(hours(1),length(dates),1)],2*[0.05 0.1 0.3])
 
% set limits to Jan 1 to Dec 31.
xlim([dates(1) dates(end)])
% set limits to time between 0:00 and 24:00 as decimal days.
ylim([hours(1) hours(end)])
 
h = gca; % get current axis.
 
% Set the XTicks to the values in month_starts.
set(h,'XTick',month_starts)
% Set the YTicks to the hours of the day.
set(h,'YTick',hours)
 
% Reformat the XTicks to 3 letter month names (Jan, Feb, etc) using 'mmm' format.
datetick(h,'x','mmm','keepticks','keeplimits')
 
% Reformat YTicks to 24 hours format 'HH:MM'
datetick(h,'y','HH:MM','keepticks','keeplimits')
 
% vertical grid lines: one at each tick mark on the x-axis.
% Get the limits of the y axis, which are datenums representing [0:00,24:00]
ylimits= get(h,'Ylim');
% for each tick mark line, we have two points: (tickmark, 0:00) and (tickmark, 24:00).
% create X = [Jan-1st Feb-1st ... Dec-1st; Jan-1st Feb-1st ... Dec-1st] (as datenums)
X = repmat(month_starts,2,1); % x-values
% create Y = [0:00 0:00 ... 0:00; 24:00 24:00 ... 24:00] (as datenums)
Y = repmat(ylimits',1,size(month_starts,2)); % y-values.
% Plot the grid lines with gray colour.
plot(X,Y,'color',[0.5 0.5 0.5])
 
% horizontal grid lines: one for each hour on the y-axis.
% Get limits of the x axis, which are datenums representing [Jan-1st,Dec-31st]
xlimits= get(h,'xlim');
% for each hour line, we have two points (as datenums): (Jan-1st, Hour) and (Dec-31st, Hour).
% create X = [Jan-1st Jan-1st ... Jan-1st; Dec-31st Dec-31st ... Dec-31st; ]
X2 = repmat(xlimits',1,size(hours,2)); % x-values
% create Y = [0:00 1:00 ... 24:00; 0:00 1:00 ... 24:00] (as datenums)
Y2 = repmat((hours),2,1); % y-values
% Plot the grid lines with gray colour.
plot(X2,Y2,'color',[0.5 0.5 0.5])
 
% finally, plot the sunset, sunrise and daylength times once so we get lines for them.
plot(dates,sunrise_times,'color',[1 0.5 0],'linewidth',3)
plot(dates,sunset_times,'color',[1 0.5 0],'linewidth',3)
plot(dates,daylength_times,'color',[1 0.9 0],'linewidth',3)
 
% Text labels for each curve.
 
% get the index where the min / max occurs for each variable.
% median is just used to deal with curves where there are multiple extreme values.
 
sunrise_imin = floor(median(find(min(sunrise_times) == sunrise_times)));
daylength_imax = floor(median(find(max(daylength_times) == daylength_times)));
sunset_imax = floor(median(find(max(sunset_times) == sunset_times)));
 
% For all 3 labels, use the daylength max as the horizontal position.
% This could be adjusted for the southern hemisphere, since the 
% daylength minimum is in the middle of the year.
 
% sunrise time label: use the sunrise min as the vertical position.
text(dates(daylength_imax),sunrise_times(sunrise_imin),'Sunrise',...
	'VerticalAlignment','bottom',...
	'HorizontalAlignment','center',...
	'FontSize',10)
 
% sunset time label: use the sunset max as the vertical position.
text(dates(daylength_imax),sunset_times(sunset_imax),'Sunset',...
	'VerticalAlignment','top',...
	'HorizontalAlignment','center',...
	'FontSize',10)
 
% daylength label: use the daylength max as the verical position.
text(dates(daylength_imax),daylength_times(daylength_imax),'Daylight',...
	'VerticalAlignment','top',...
	'HorizontalAlignment','center',...
	'FontSize',10)
 
title(sprintf('Sunrise and Sunset Times\n%s (%d)', city, year))