%% Grouped Calculations in Tables and Timetables
% Grouped calculations can help you interpret large datasets such as time-series data. In such calculations, you use a grouping variable to split a dataset into groups and apply a function to each group. A grouping variable contains values, such as time periods or station locations, that you can use to group other data values, such as temperature readings or atmospheric concentrations of a gas. In MATLAB®, you can store such data in tables or timetables. With grouped calculations in a table you can often calculate results in-place, in one table, instead of breaking data out into separate tables and merging results later.
% This example shows how to import nitrogen dioxide (NO2) data from the US Environmental Protection Agency (EPA) into a table and do grouped calculations on this data. NO2 is one of the Criteria Air Pollutants regulated under the US Clean Air Act. It is toxic by itself and is also a key component of photochemical smog that results in ground-level ozone production. NO2 is produced through high-temperature processes that can split nitrogen and oxygen gases and enable them to recombine. Natural processes contribute NO2 to the atmosphere, but so do human activities such as combustion in automobile engines and power plants, lightning, and biomass burning. The concentration of NO2 in the atmosphere is also influenced by the photochemical cycling between NO and NO2, atmospheric transport, and ultimately oxidation to nitric acid, causing acid rain. Different processes contribute NO2 to the atmosphere on different timescales, leading to daily (diurnal), weekly, and annual cycles in its atmospheric concentration. Time-series analysis of such data relies heavily on grouped calculations to examine different periodic behavior or to average the data over time to smooth out high-frequency variability and reveal long-term trends.
% The example first shows how to do preliminary data cleaning, including conversion of the table to a timetable. Then it shows simple ways to group the data by one grouping variable and calculate annual mean NO2 concentrations. It also shows how to group the NO2 data by two grouping variables together, time and location, enabling calculations that find locations exceeding EPA standards at various times. You can also group the NO2 data by time period to look for daily or yearly cycles. Finally it shows how to apply a function that requires inputs from multiple table variables to find the times at which the maximum NO2 concentrations occurred at each site.
%% Import NO2 Data to Table
% First, import NO2 data from the Air Quality System (AQS) database maintained by the EPA. This data consists of hourly measurements of NO2 concentrations from outdoor monitors across the United States, Puerto Rico, and the U.S. Virgin Islands. It is stored as a set of zipped spreadsheets, one for each year starting with 1980.
% Download hourly NO2 measurements for the years 1985–1989. You can download and unzip the compressed spreadsheets by using the unzip function. The result is set of files in your current folder with names such as hourly_42602_1985.csv. Here, 42602 is an EPA code for NO2. (Data from the US Environmental Protection Agency. Air Quality System Data Mart available via https://www.epa.gov/airdata. Accessed July 15, 2021.)

yrs = string(1985:1989);
urls = "https://aqs.epa.gov/aqsweb/airdata/hourly_42602_" + yrs + ".zip";
fnames = strings(numel(yrs),1);
for ii = 1:numel(yrs)
    fnames(ii) = unzip(urls(ii));
end
fnames

%
% Import data from the spreadsheets into a table. Start by creating an empty table. Then import data from the spreadsheets, one by one, by using the readtable function and adding it to the table.
% Create import options that help specify how readtable imports tabular data. To create import options based on the contents of the spreadsheets, use the detectImportOptions function. Read all the text data into table variables that store strings. You can also specify that only specified table variables have certain data types. To specify that only the TimeGMT and TimeLocal table variables store times as duration arrays, use the setvaropts function.

NO2data = table;
opts = detectImportOptions(fnames(1),"TextType","string");
opts = setvaropts(opts,["TimeGMT","TimeLocal"],"Type","duration","InputFormat","hh:mm");

%
% Import data from the spreadsheets by using the readtable function. You can vertically concatenate the tables you read in so that all the data is in one large table.
% The spreadsheets have column names, such as "Time GMT", that you cannot use as MATLAB identifiers. As the warning messages indicate, readtable converts these names into table variable names that are valid MATLAB identifiers, such as TimeGMT. When a table variable name is also a valid MATLAB identifier, it is easier to access the variable by using dot notation, as in NO2data.TimeGMT.

for ii = 1:numel(yrs)
    NO2data = [NO2data; readtable(fnames(ii),opts)];
end

%
% Display NO2data. It has 24 variables storing NO2 sample measurements, site locations, state names, times, and many other pieces of information.

NO2data

%
%% Clean NO2 Table and Convert to Timetable
% Next, prepare NO2data for analysis by cleaning the data. Data cleaning is the process of detecting and correcting (or removing) parts of the data set that are either corrupt, inaccurate, or irrelevant. You can also convert table variables so that they have data types that can be more convenient for analysis, such as categorical or datetime arrays.
% For example, the table variable SampleMeasurement has measurements of NO2 concentration. Concentrations below the method detection limit (MDL) are unreliable. To exclude them from analysis, find the rows where SampleMeasurement is below the MDL. Set those elements to NaN.

NO2data.SampleMeasurement(NO2data.SampleMeasurement < NO2data.MDL) = NaN;

%
% Create a table that contains only the subset of variables that are relevant to this example. You can use table subscripting to create a table that has all rows (specified by a colon) and only those variables that you name.

NO2data = NO2data(:,["DateLocal","TimeLocal","SampleMeasurement","StateName","CountyName","SiteNum","Latitude","Longitude"])

%
% Combine the local date and time into a single timestamp. The new Timestamp table variable is a datetime array. Delete the DateLocal and TimeLocal variables because they are now redundant.

NO2data.Timestamp = NO2data.DateLocal + NO2data.TimeLocal;
NO2data.Timestamp.Format = "default";
NO2data.DateLocal = [];
NO2data.TimeLocal = [];

%
% To categorize the data later, convert the StateName and CountyName variables to categorical arrays, first erasing space characters from the names. There are fixed sets of state and county names in the data, which makes it convenient to create categories based on them.

NO2data.StateName = categorical(erase(NO2data.StateName," "));
NO2data.CountyName = categorical(erase(NO2data.CountyName," "));

%
% Rename the SampleMeasurement variable to MeasuredNO2. One way to rename table variables is by using the VariableNames property of the table.

NO2data.Properties.VariableNames("SampleMeasurement") = "MeasuredNO2";

%
% Convert NO2data to a timetable. The datetime values in Timestamp are now row times that label the rows of the timetable. The dates and times of the original table were in separate variables. To put data like this data into a timetable, it is more convenient to import the data as a table, and then combine the separate date and time variables into one datetime variable. Then convert the modified table by using the table2timetable function.

NO2data = table2timetable(NO2data)

%
%% Simple Grouped Calculations by State
% Given the size of the timetable, it is obvious that there are many thousands of hourly measurements in every state. One way to calculate the number of measurements for each state is to sum the number of rows that have a particular state as a category. For example, calculate the number of measurements for Alaska, and then for Arizona.

numAlaska = sum(NO2data.StateName=="Alaska")
numArizona = sum(NO2data.StateName=="Arizona")

%
% It is tedious to perform this calculation multiple times or to store intermediate results in many variables or subtables. Instead, MATLAB provides functions that group data in tables and apply functions to each group in-place. For example, use the groupcounts function to group the data in NO2data by the states in StateName and count the rows in each group. Instead of calling sum many times, call groupcounts once.

NO2counts = groupcounts(NO2data,"StateName")

% To sort the results in a table or timetable, use the sortrows function. Sort gc on its GroupCount variable from highest to lowest value.

sortedNO2counts = sortrows(NO2counts,"GroupCount","descend")

% To calculate other statistics, use the groupsummary function. For example, find the maximum NO2 concentration measured in each state.

NO2max = groupsummary(NO2data,"StateName","max","MeasuredNO2");
sortedNO2max = sortrows(NO2max,"max_MeasuredNO2","descend")

% As an alternative, you can use the varfun function with the "GroupingVariables" name-value argument for grouped calculations. But the groupsummary function is simpler and performs most of the same grouped calculations as varfun.
%
%% Simple Grouped Calculations by Time
% Functions such as groupcounts, groupsummary, and varfun work equally well on tables and timetables. But timetables also provide the retime and synchronize functions, which can perform time-based calculations by using their row times. You can group timetable data by time and perform calculations on data within the time periods. The retime function is the best option for such cases.
% For example, group the data in NO2data into yearly time periods. Find the maximum NO2 concentration for each year.

yearlyMaxNO2 = retime(NO2data(:,"MeasuredNO2"),"yearly","max")

% This calculation is useful if you have one time series. In this case, the data in the MeasuredNO2 variable come from multiple sites. A more useful analysis is to group by both year and site.
%
%% Calculate Annual Means by Site
% The US EPA has two National Ambient Air Quality Standards (NAAQS) for NO2. A location is not in compliance with the NAAQS if either:
% The annual mean exceeds 53 ppb
% The 98th percentile of 1-hour daily maximum concentrations, averaged over 3 years, exceeds 100 parts-per-billion (ppb)
% Analyze data in NO2data to find locations that are not in compliance with the first standard, where the annual mean exceeded 53 ppb. There are three different ways to approach this analysis. What the three approaches have in common is that you can group the data by both time and site to calculate annual means by site.
%
%% Group by Multiple Grouping Variables
% To find sites that do not comply with the NAAQS, calculate the mean value for each site for each year. While NO2data does not include unique identifiers for the sites, you can use state names, county names, and site numbers together to uniquely identify air quality sites.
% The row times of NO2data are datetime values. Extract their year components and add a new variable to NO2data named Year. Calculate the annual means for each site by using groupsummary with StateName, CountyName, SiteNum, and Year as grouping variables.

NO2data.Year = year(NO2data.Timestamp);
meanNO2bySite = groupsummary(NO2data,["StateName","CountyName","SiteNum","Year"],"mean","MeasuredNO2")

% To find the sites that have the highest mean NO2, sort the timetable.

sortedMeanNO2bySite = sortrows(meanNO2bySite,"mean_MeasuredNO2","descend")

% You can create a table that includes only those sites exceeding 53 ppb by using logical indexing. Create a logical vector that indicates the rows where mean_MeasuredNO2 is greater than 53. Use that vector as a subscript to get matching rows from meanNO2bySite.

exceeded53ppb = meanNO2bySite.mean_MeasuredNO2 > 53;
sitesExceed53ppb = meanNO2bySite(exceeded53ppb,:)

%
%% Pivot to Find Relationships Between Grouping Variables
% Sometimes pivoting, or rearranging statistics calculated from tabular data, makes it easier to see and analyze results, particularly when you look at the relationship between two grouping variables. For example, you can create a pivot table for the annual mean NO2 by site. By pivoting, you can create a table where every site lists annual mean NO2 in its own table variable, showing the relationship between year and site. In MATLAB, you can create pivot tables by using the stack and unstack functions, which stack and unstack table variables into taller or wider formats.
% A complication in this case is that NO2data has three grouping variables that together uniquely identify sites: state name, county name, and site number. To create a pivot table, first combine these three table variables into one variable. Convert StateName, CountyName, and SiteNum into strings and add them together. Replace spaces and dashes with underscores, and erase periods and parentheses. The names in SiteID are unique site identifiers.

siteID = string(NO2data.StateName) + "_" + string(NO2data.CountyName) + "_" + string(NO2data.SiteNum);
siteID = replace(siteID,[" ","-"],"_");
siteID = erase(siteID,[".","(",")"]);

%
% Add SiteID to NO2data as a new table variable. Calculate annual means by using groupsummary, but this time use SiteID as a grouping variable.

NO2data.SiteID = categorical(siteID);
meanNO2bySiteID = groupsummary(NO2data,["SiteID","Year"],"mean","MeasuredNO2")

%
% To create a pivot table, use the unstack function. Each unique site in the SiteID variable of meanNO2bySiteID becomes the name of a separate table variable in the output, pivotedMeanNO2bySiteID, and has the annual means associated with that site. This unstacking operation is how you can create a pivot table in MATLAB.

pivotedMeanNO2bySiteID = unstack(meanNO2bySiteID,"mean_MeasuredNO2","SiteID","GroupingVariable","Year")

% This representation of the annual means by site has an advantage and a disadvantage.
% It is easier to look at the short five-year time series for each site. After unstacking, each site has its own variable in pivotedMeanNO2bySiteID. You can easily compare sites to each other.
% It is harder to sort and pick out the largest values across the whole pivoted table. After unstacking, pivotedMeanNO2bySiteID has 443 variables. The stacked version, meanNO2bySite, has only seven variables.
%
%% Group by Time and Another Grouping Variable
% To group data in NO2data by year and another grouping variable, it was necessary to add Year as an additional variable. Also, the output from groupsummary is a table even when the input is a timetable. But suppose you want to keep the results in a timetable instead. The retime function can also produce annual summaries. But it can group data only by time. To group data by site and by year, rearrange NO2data so that you can call retime on a timetable where the NO2 concentrations are already grouped by site.
% Group the raw data in NO2data by site by using the unstack function. The output timetable has a separate variable for each site. This timetable looks similar to a pivot table. But instead of having means or some other statistic, NO2bySite has all the raw data. It is just reorganized. For further convenience, sort the rows of the timetable by their row times so that the earliest timestamps come first.

NO2bySite = unstack(NO2data,"MeasuredNO2","SiteID","GroupingVariable","Timestamp");
NO2bySite = sortrows(NO2bySite)

%
% In this format you can easily plot the raw data by using the stackedplot function. This plot shows NO2 concentrations for each site as a function of time.

stackedplot(NO2bySite)

%
% To create a timetable that is also a pivot table, use retime to calculate annual means.

meanNO2bySiteTT = retime(NO2bySite,"yearly","mean")

%
% Create a stacked plot of the annual means.

stackedplot(meanNO2bySiteTT)

%
% You might want to preserve information from the original timetable NO2data in this timetable of results. For example, you might want to add the latitudes and longitudes of the sites to NO2bySite. They were stored for each timestamp in NO2data. But to store them more compactly in this timetable, add them as per-variable custom properties to NO2bySite.

LatLon = groupsummary(NO2data,"SiteID","mode",["Latitude","Longitude"]);
NO2bySite = addprop(NO2bySite,["Latitude","Longitude"],["variable","variable"]);
NO2bySite.Properties.CustomProperties.Latitude(string(LatLon.SiteID)) = LatLon.mode_Latitude';
NO2bySite.Properties.CustomProperties.Longitude(string(LatLon.SiteID)) = LatLon.mode_Longitude';

%
%% Moving Means for Data Grouped by Site
% To calculate compliance with the second NAAQS standard for NO2 requires a sequence of grouped calculations. By the second standard, a location is out of compliance if the 98th percentile of the 1-hour daily maximum concentrations of NO2, averaged over 3 years, exceeds 100 ppb.
% Start with the hourly concentrations of NO2 by site. To find the daily maximum for each site, use the retime function, specifying "max" as the method to find the maximum concentration for each day's worth of data. Then find the 98th percentiles of the daily maximums in each year's worth of data, calling retime a second time. To calculate percentiles, use the findPrctile supporting function referred to in this example.

dailyMax = retime(NO2bySite,"daily","max");
yearlyP98 = retime(dailyMax,"yearly",@(x)findPrctile(x,98))

%
% Next calculate a moving mean for each site, specifying a three-year window for the moving mean. The smoothdata enables you to apply the movmean function to each variable in yearlyP98.

moving3yearAvg = smoothdata(yearlyP98,"movmean",[years(3) 0])

%
% Display sites that are out of compliance. First specify a time range starting in 1987, the first year for which the moving three-year window has three full years of data.

full3years = timerange("1987-01-01","1989-01-01","closed")

%
% Next find sites that exceeded the standard during any year in the 1987–1989 time period. The vector exceed is a vector of logical values whose values are 1 (true) where the corresponding variables of moving3yearAvg have values that exceed the standard. You can use logical arrays to index into tables. In this case, index into moving3yearAvg by using exceed to display only the variables for the sites that exceed the standard.

exceed = any(moving3yearAvg{full3years,:}>100,1);
moving3yearAvg(full3years,exceed)

%
% The NO2bySite timetable has latitudes and longitudes for the sites, saved as custom properties. The latitudes and longitudes are associated with the variables of the timetable. You can mark the sites that exceeded the standard on a map by using the geoscatter function. The sites that exceeded the standard in 1987 are marked in yellow. (Here, moving3yearAvg{1,:} accesses the first row of the timetable, corresponding to the year 1987.) The sites out of compliance in 1987 were typically large cities such as Los Angeles. Due to the Clean Air Act, a similar analysis using the latest data for the years 2015–2019 would show no sites out of compliance with the NO2 standard.

geoscatter(moving3yearAvg.Properties.CustomProperties.Latitude,...
           moving3yearAvg.Properties.CustomProperties.Longitude,...
           moving3yearAvg{1,:},moving3yearAvg{1,:}>100,'filled')

%
%% Group by Time Periods
% Another way to group by time is by using periodic time units to look for things like seasonality or daily cycles. For example, consider the average daily pattern of NO2 concentrations at each site. It is likely that fossil fuel combustion emissions from human activity and atmospheric photochemistry driven by the sun contribute to a daily cycle in NO2 concentrations. You cannot calculate the mean daily cycle using retime. One approach is to use the varfun function by adding a grouping variable based on the time of day of each timestamp. The timeofday function returns the time of day, or length of time since midnight, as a duration. This code sample shows the approach by using varfun.
% NO2bySite.Hour = timeofday(NO2bySite.Timestamp);
% NO2bySite.Hour.Format = "hh:mm";
% meanDailyCycleNO2 = varfun(@mean,NO2bySite,"GroupingVariable","Hour")
% For common calendar periods such as hour of day or month of year, there is a simpler method. These common calendar periods are options of the groupsummary function. Call groupsummary with the "hourofday" option.

meanDailyCycleNO2 = groupsummary(NO2bySite,"Timestamp","hourofday","mean")

% To visualize the daily cycle for the first few sites, use stackedplot. Compare the pattern of a single peak in remote sites in Alaska and Arizona to the morning and evening peaks during rush hour at urban sites in Arizona and California.

meanDailyCycleNO2.GroupCount = [];
meanDailyCycleNO2.hourofday_Timestamp = hours(double(meanDailyCycleNO2.hourofday_Timestamp));
stackedplot(meanDailyCycleNO2,'XVariable','hourofday_Timestamp')

%
%% Grouped Calculations with Functions That Require Multiple Inputs
% The retime, groupsummary, and varfun functions all apply functions separately to each table variable. But sometimes you have functions that use more than one table variable as inputs. For example, you might want to find the time or index at which some condition occurred within each group of data values. In such cases, use the rowfun function. It enables you to apply functions that require multiple inputs.
% For example, determine when the maximum NO2 concentration occurred at each site. This determination requires a function such as the findMax supporting function referred to in this example. The findMax function requires both timestamps and data values as input arguments. It returns the maximum value with the time at which the maximum value occurred.
% To group the data in NO2data by SiteID and find the times when the maximum NO2 concentration occurred at each site, use rowfun. Specify that the inputs to findMax are Timestamp and MeasuredNO2 from NO2data. Convert NO2data to a table so that rowfun returns a table.

NO2data = timetable2table(NO2data);
rowfun(@findMax,NO2data,"GroupingVariable","SiteID","InputVariables",["Timestamp","MeasuredNO2"], ...
    "OutputVariableNames",["MaxMeasuredNO2","MaxOccurrenceTime"])

% With these results you could extend your analysis to find out why the NO2 concentrations were particularly high on those dates.
%
%% Supporting Functions
% Supporting local functions are defined below.
%

function [maxVal,maxTime] = findMax(times,vals)
    % Return time at which maximum element of vals occurred
    [maxVal,maxIndex] = max(vals);
    if ~isnan(maxVal)
        maxTime = times(maxIndex);
    else
        maxTime = NaT;
    end
end
function y = findPrctile(x,p)
    % Return data point nearest to percentile p, without interpolation
    xs = sort(x);
    n = sum(~isnan(x)); % use non-NaN elements only
    k = p*n/100 + 0.5;  % index of data point that represents 100*(i-0.5)/n th percentile
    y = xs(round(k));   % data point nearest specified p
end

% Copyright 2021 The MathWorks, Inc.