MATLAB Examples

## Contents

```% ThingSpeak Weather Station Data Analysis % % This script includes examples of reading data from ThingSpeak channel and % performing various kind of visualization on the data. % % Retrieve the data first (see the first section) before executing other sections. % % ThingSpeak Support Toolbox on File Exchange is required to run this script. % It can be downloaded here: % http://www.mathworks.com/matlabcentral/fileexchange/52244-thingspeak-support-toolbox % ```

## Retrieve data from ThingSpeak channel

```% Channel ID to read data from readChannelID = 87179; % Specify date range dateRange = [datetime('March 7, 2016'),datetime('March 13, 2016')]; % Read data including the timestamp, and channel information. [data,time,channelInfo] = thingSpeakRead(readChannelID,'Fields',1:7,... 'DateRange',dateRange); % Create variables to store different sorts of data temperatureData = data(:,1); humidityData = data(:,2); pressureData = data(:,3); rainData = data(:,4); windSpeedData = data(:,5); windGustData = data(:,6); windDirectionData = data(:,7); ```

## Temperature, Humidity, Pressure, Rain, WindSpeed, WindDirection histogram

```% Create a figure to display plots figure % Temperature histogram subplot(2,3,1) % Create 2-by-3 axis on the same figure, and work on the first axis histogram(temperatureData); title(channelInfo.FieldDescriptions{1}); grid on % Humidity histogram subplot(2,3,2) histogram(humidityData); title(channelInfo.FieldDescriptions{2}); grid on % Pressure histogram subplot(2,3,3) histogram(pressureData); title(channelInfo.FieldDescriptions{3}); grid on % Rain fall histogram subplot(2,3,4) histogram(rainData); title(channelInfo.FieldDescriptions{4}); grid on % WindSpeed histogram subplot(2,3,5) histogram(windSpeedData); title(channelInfo.FieldDescriptions{5}); grid on % Wind Direction histogram rad = windDirectionData*pi/180; % Convert to radians rad = -rad+pi/2; % Adjust the wind direction data to match map compass, such that North is equal to 0 degree subplot(2,3,6) rose(rad,12) % Plot the wind direction histogram in a polar axis title(channelInfo.FieldDescriptions{7}) ax = gca; ax.View = [-90 90]; % Rotate axis 90 degrees counterclock-wise such that North is equal to 0 degree ```

## Smooth Temperature and Trend

```% Remove missing data from the temperature variable, in order to perform % the fitting methods idx = ~isnan(temperatureData); rawTemp = temperatureData(idx); newTime = time(idx); % Smooth the raw temperature data with local 60-point mean values smoothTemp = movmean(temperatureData(idx),60); % Fit the data for a trend line [p,~,mu]= polyfit(datenum(newTime),rawTemp,1); trend = polyval(p,datenum(newTime),[],mu); % Optional: fit the data by using the sum of eight sine functions % To use the fit function, the Curve Fitting Toolbox is required. % It can be downloaded here: http://www.mathworks.com/products/curvefitting/ f = fit(datenum(newTime),rawTemp,'sin8'); % Plot the raw data, smooth data, trend and fitting curve figure hold on plot(newTime,rawTemp,'b') % Plot the raw temperature plot(newTime,smoothTemp,'g','LineWidth',1.5) % Plot the smoothing curve plot(newTime,trend,'m','LineWidth',2) % Plot the trend line plot(f) % Plot the fitting curve hold off xlim([datenum(time(1)) datenum(time(end))]) % Adjust the x-axis to properly display the curves xlabel('Date') ylabel('Temperature (F)') legend({'Raw Data','Smooth Data','Trend','Fitting Curve'},'Location','NE') ```

## Daily Temperature Statistics

```% Remove missing data idx = ~isnan(temperatureData); rawTemp = temperatureData(idx); newTime = time(idx); % Plot temperature for each day. figure hold on % Use dateshift and findgroups functions to bin the temperature data in % terms of individual day timeShift = dateshift(newTime,'start','day'); % Shift hour, minute, and second to the start of each day, i.e. 00:00:00 dayGroup = findgroups(timeShift); % Group the shifted time in terms of day splitapply(@plot,second(newTime,'secondofday')/3600,rawTemp,dayGroup); % Apply the plot function to the temperature for each day by using the group defined above % Compute the absolute Max, Min, and Mean for each hour hourGroup = findgroups(hour(newTime)); % Create a vector of grouping number for each hour tMax = splitapply(@max,rawTemp,hourGroup); % Calculate the max for each hour tMin = splitapply(@min,rawTemp,hourGroup); % Calculate the min for each hour tMean = splitapply(@mean,rawTemp,hourGroup); % Calculate the mean for each hour % Prepare x and y data for the average and the variation area xHalf = [0,repelem(1:23,2),24]; % x coordinate at the end points for each interval (each hour) x = [xHalf, fliplr(xHalf)]; % x coordinate for variation area (including both bottom and top profiles) tMax = repelem(tMax,2); % y coordinate for max per hour tMin = repelem(tMin,2); % y coordinate for min per hour tMean = repelem(tMean,2);% y coordinate for mean per hour y = [tMin; flipud(tMax)]; % y coordinate for variation area (bottom and top) % Plot Statistics plot([0,24], [max(tMax),max(tMax)],'.--k', 'LineWidth',1.5) % Plot global Max as a horizontal line plot([0,24], [min(tMin),min(tMin)],'.--k', 'LineWidth',1.5) % Plot global Min as a horizontal line h1 = fill(x,y, 'b', 'LineStyle', 'none', 'FaceAlpha', 0.1); % Highlight the variation area per hour h2 = plot(xHalf,tMean,'--','LineWidth',2); % Plot average per hour hold off title('Daily temperature over the past week') ylabel('Temperature (F)') xlabel('Day Time') axis tight % Adjust axis to fit the plot ylim([20 80]) % Adjust to allow space on the top and at the bottom of the axis ax = gca; % Get the current axis object ax.XTick = 1:24; % Change the X-Tick to 24 hours legend([h1,h2],'Variation per hour','Average per hour') ```

## Compare with historical temperature data on March 7, 2011-2015

```% Load historical data rawData = load('March7'); % Load the data file named "March7" in the current folder. rawData = struct2cell(rawData); % Convert the structure rawData to cell, which allows numeric indexing without knowing the field name. % Calculate hourly Max, Min, and Mean for historical data m = zeros(24,length(rawData)); % Pre-allocate matrix for i = 1:length(rawData) % Loop over all years (2011-2015) histData = rawData{i}; % Get the temperature data for the i-th year hourGroup = findgroups(hour(histData.TimeEST)); % Group for each hour m(:,i) = splitapply(@mean,histData.TemperatureF,hourGroup); % Calculate the average temperature for each hour end histMax = max(m,[],2); % Max per hour histMin = min(m,[],2); % Min per hour histMean = mean(m,2); % Mean per hour % Temperature on March 7,2016 per hour temp = temperatureData(month(time)==3 & day(time)==7); s = second(time(month(time)==3 & day(time)==7),'secondofday'); % Coresponding time slot in seconds s = s/3600; % Convert to hours % Plot figure % In the first subplot, show the historical temperatures for each year subplot(2,1,1) plot(m, 'LineWidth',2) title('Daily Temperature on March 7, 2011-2015') ylabel('Temperature (F)') legend({'2011','2012', '2013', '2014','2015'}) ax = gca; ax.XTick = 1:24; % Change the X-Tick to 24 hours % In the second subplot, compare the current temperature with historical data subplot(2,1,2) hold on plot(s,temp,'g', 'LineWidth',2) % Plot temperature on March 7, 2016 plot(histMax,'--','LineWidth',2) % Plot historical max plot(histMin,'--','LineWidth',2) % Plot historical min plot(histMean,'--','LineWidth',2) % Plot historical mean hold off title('Daily Temperature on March 7 - Current v.s. History') ylabel('Temperature (F)') legend({'March 7, 2016','Historical Max', 'Historical Min', 'Historical Mean'}) ax = gca; ax.XTick = 1:24; % Change the X-Tick to 24 hours ```

## Temperature and Humidity 3D bar charts

```% Create a day range vector dayRange = day(dateRange(1):dateRange(2)); % Pre-allocate matrix weatherData = zeros(length(dayRange),24); % Generate temperature 3D bar chart % Get temperature per whole clock for each day for m = 1:length(dayRange) % Loop over all days for n = 1:24 % Loop over 24 hours if any(day(time)==dayRange(m) & hour(time)==n); % Check if data exist for this specific time hourlyData = temperatureData((day(time)==dayRange(m) & hour(time)==n)); % Pull out the hourly temperature from the matrix weatherData(m,n) = hourlyData(1); % Assign the temperature at the time closest to the whole clock end end end % Plot figure h = bar3(datenum(dateRange(1):dateRange(2)), weatherData); for k = 1:length(h) % Change the face color for each bar h(k).CData = h(k).ZData; h(k).FaceColor = 'interp'; end title('Temperature Distribution') xlabel('Hour of Day') ylabel('Date') datetick('y','mmm dd') % Change the Y-Tick to display specified date format ax = gca; ax.XTick = 1:24; % Change the X-Tick to 24 hours ax.YTickLabelRotation = 30; % Rotate label for better display colorbar % Add a color bar to indicate the scaling of color % Generate humidity 3D bar chart % Get humidity per whole clock for each day for m = 1:length(dayRange) % Loop over all days for n = 1:24 % Loop over 24 hours if any(day(time)==dayRange(m) & hour(time)==n); % Check if data exist for this specific time hourlyData = humidityData((day(time)==dayRange(m) & hour(time)==n)); % Pull out the hourly humidity from the matrix weatherData(m,n) = hourlyData(1); % Assign the humidity at the time closest to the whole clock end end end % Plot figure h = bar3(datenum(dateRange(1):dateRange(2)), weatherData); for k = 1:length(h) % Change the face color for each bar h(k).CData = h(k).ZData; h(k).FaceColor = 'interp'; end title('Humidity Distribution') xlabel('Hour of Day') ylabel('Date') datetick('y','mmm dd') % Change the Y-Tick to display specified date format ax = gca; ax.XTick = 1:24; % Change the X-Tick to 24 hours ax.YTickLabelRotation = 30; % Rotate label for better display colorbar % Add a color bar to indicate the scaling of color ```

## Interpolation and contour for Temperature, Humidity and Pressure

```% Replace missing data by interpolation, rather than removing the missing % data directly from the variable. This allows to keep the dimension of the % array being consistent xNew = linspace(1,size(data,1),100)'; % Create new x coordinates tNew = interp1(temperatureData(~isnan(temperatureData)),xNew,'linear','extrap'); % Temperature interpolation. Extrapolation is applied here in case that the last entry is NaN. hNew = interp1(humidityData(~isnan(humidityData)),xNew,'linear','extrap'); % Humidity interpolation pNew = interp1(pressureData(~isnan(pressureData)),xNew,'linear','extrap'); % Pressure interpolation % Find the index of the max pressure [pMax,idx] = max(pNew); % Create surface fitting data sf = fit([tNew,hNew],pNew,'linearinterp'); % Plot figure hsf = plot(sf,[tNew,hNew],pNew); % Plot the surface with nodes. This plot function is provided in Curve Fitting Toolbox. The output is an array of a surface object and a line object. hsf(1).EdgeColor = 'interp'; % Change face edge color of the surface hsf(1).FaceAlpha = 0.5; % Change the transparency of the surface xlabel('Temperature') ylabel('Humidity') zlabel('Pressure') title('Linear Interpolation Surface') % 2D View with the location of max pressure figure hsf = plot(sf); % Plot the surface only hsf.EdgeColor = 'interp'; % Change face edge color hold on plot3(tNew(idx),hNew(idx),pMax,'r.', 'MarkerSize',30) % Plot the location of max pressure text(tNew(idx)+2,hNew(idx)+2,pMax,['P= ',num2str(pMax),', T=',... num2str(tNew(idx)),', H=',num2str(hNew(idx))]) % Display the values at the location above title('Contour of the Pressure') xlabel('Temperature') ylabel('Humidity') grid off view(2) % Set the view to 2D, i.e., observing the plot from top to bottom along z-axis hold off ```

## Dew Point

```% Convert temperature from Fahrenheit to Celsius tempC = (5/9)*(temperatureData-32); % Calculate dew point. Refer to Wiki (https://en.wikipedia.org/wiki/Dew_point) for the formula and more details % Specify the constants for water vapor (b) and barometric (c) pressure. b = 17.67; c = 243.5; % Calculate the intermediate value 'gamma' gamma = log(humidityData/100) + b*tempC ./ (c+tempC); % Calculate dew point in Celsius dewPoint = c*gamma ./ (b-gamma); % Convert to dew point in Fahrenheit dewPointF = (dewPoint*1.8) + 32; % Plot figure hold on plot(time, dewPointF,'d') % Plot dew points xlabel('Time') ylabel('Dew Point') xlim([datenum(time(1)) datenum(time(end))]) % Adjust the x-axis limits fill([xlim fliplr(xlim)], [68 68 80 80], 'r', 'LineStyle', 'none', 'FaceAlpha', 0.1) % Highlight the uncomfortable zone text(0.7*datenum(time(1)) + 0.3*datenum(time(end)), 75, 'Uncomfortable', 'FontWeight','bold') % Add text for the zone fill([xlim fliplr(xlim)], [50 50 68 68], 'g', 'LineStyle', 'none', 'FaceAlpha', 0.1) % Highlight the comfortable zone text(0.7*datenum(time(1)) + 0.3*datenum(time(end)), 60, 'Comfortable', 'FontWeight','bold') % Add text for the zone fill([xlim fliplr(xlim)], [min(ylim) min(ylim) 50 50], 'y', 'LineStyle', 'none', 'FaceAlpha', 0.1) % Highlight the dry zone text(0.65*datenum(time(1)) + 0.35*datenum(time(end)), 20, 'Dry', 'FontWeight','bold') % Add text for the zone hold off ```

## Wind Compass and Feather

```% Specify the latest n+1 wind directions to be displayed n = 9; % Convert to radians rad = windDirectionData*pi/180; % Create a feather plot % Remove missing data and any wind speed with value 0 idx = (~isnan(rad)) & (~isnan(windSpeedData)) & (windSpeedData~=0); % Convert polar coordinates to Cartesian. Note that dividing by the maximum % wind speed allows to scale the length of each arrow by its relative wind % speed, rather than the wind direction. [x,y] = pol2cart(rad(idx),windSpeedData(idx)/max(windSpeedData(idx))); % Plot figure subplot(2,1,2) feather(x((end-n):end),y((end-n):end)) % Plot the feather xlim([0 n+2]) % Adjust the x-axis ylim([-1 1]) % Adjust the y-axis xlabel(['The last ',num2str(n+1),' wind direction']) % Add x lable title('Wind Direction Changes') grid on ax = gca; ax.YTickLabel = {}; % Hide the Y-Tick ax.XTick = 1:(n+1); % Adjust the X-Tick for n+1 data % Create a compass plot % Adjust the wind direction to match map compass, such that North is equal to 0 degree rad = -rad+pi/2; % Calculate the cosine component u = cos(rad) .* windSpeedData; % x coordinate of wind speed on circular plot % Calculate the sine component v = sin(rad) .* windSpeedData; % y coordinate of wind speed on circular plot % Plot subplot(2,1,1) compass(u((end-n):end),v((end-n):end)) % Plot compass title('Wind Compass') ax = gca; ax.View = [-90 90]; % Rotate axis 90 degrees counterclock-wise such that North is equal to 0 degree ```

## Instant Weather Flag

```% Convert to radians theta = windDirectionData(end)*pi/180; % Covert wind speed to a relative angle of the wind flag with z axis, by % assuming the max wind speed is 5mph. maxWindSpeed = 5; speedAngle = pi/2*windSpeedData(end)/maxWindSpeed; ytheta = pi-speedAngle; % Create the cone as a flag % Specify the radius of the cone ConeRadius = 0.1; % Create coordinate limits and rotation matrix xMax = 1.5; xMin = -xMax; yMax = 1.5; yMin = -yMax; zMax = 1.5; zMin = -0.5; yRotate = [cos(ytheta) 0 sin(ytheta); 0 1 0; -sin(ytheta) 0 cos(ytheta)]; zRotate = [cos(theta) -sin(theta) 0; sin(theta) cos(theta) 0; 0 0 1]; % Create the cone surface data t = [0;ConeRadius]; [X,Y,Z]=cylinder(t,50); % Adjust the direction of the cone flag such that the cone is pointing down % along z-axis XYZ = zRotate*(yRotate*[X(2,:);Y(2,:);Z(2,:)]); X(2,:) = XYZ(1,:); % Pull out x coordinate of the cone Y(2,:) = XYZ(2,:); % Pull out y coordinate of the cone Z(2,:) = XYZ(3,:); % Pull out z coordinate of the cone Z = Z + zMax; % Raise the flag to the top of the post % Plot figure c = jet; % Define the colormap % Plot the cone and use the instant temperature data to display its color h = surf(X,Y,Z,'FaceColor', c(round(temperatureData(end)/100*64),:),... 'LineStyle','none','FaceAlpha',0.8,'FaceLighting','gouraud'); % Add lights to adjust the color of the cone light('Position',[xMin 0 0],'Color',c(round(temperatureData(end)/100*64),:)) light('Position',[xMax 0 0],'Color',c(round(temperatureData(end)/100*64),:)) light('Position',[0 yMax 0],'Color',c(round(temperatureData(end)/100*64),:)) light('Position',[0 yMin 0],'Color',c(round(temperatureData(end)/100*64),:)) % Plot directions and the vertical post hold on plot3(0,0,zMax,'k.', 'MarkerSize',25) % Post plot3([0 0],[0 0],[0 zMax],'k','LineWidth',3) % Top node of the post plot3([xMin xMax],[0 0],[0 0],'k') % WE directions plot3([0 0],[yMin yMax],[0 0],'k') % NS directions % Plot auxiliary projection to help identify the direction where the flag % points to plot3([cos(theta)*sin(speedAngle) cos(theta)*sin(speedAngle)],... [0 sin(theta)*sin(speedAngle)],[0 0],'k--') plot3([0 cos(theta)*sin(speedAngle)],[sin(theta)*sin(speedAngle) sin(theta)*sin(speedAngle)],[0 0],'k--') plot3([cos(theta)*sin(speedAngle) cos(theta)*sin(speedAngle)],... [sin(theta)*sin(speedAngle) sin(theta)*sin(speedAngle)],[0 zMax-cos(speedAngle)],'k--') hold off % Add direction letters, color and annotations text([0,0,xMax,xMin,0], [yMax,yMin,0,0,0], [0.05,0,0,-0.05,zMax+0.15],... {'N','S','E','W',['Wind Speed = ',num2str(round(windSpeedData(end),1)),'mph']}) colormap('jet') cb = colorbar; % Get the colorbar object colorbar('Ticks',linspace(cb.Limits(1), cb.Limits(2),6),... 'TickLabels',{'Freeze','Cold','Cool','Neutral','Warm','Hot'}) % Label the color bar in terms of the human feeling about the temperature title('Weather Flag') axis equal axis off ```