<p><div><div><div><div data-width="100" data-show-faces="false" data-send="true" data-layout="button_count" data-height="21" data-href="https://www.r-bloggers.com/introducing-the-pwfslsmoke-package/"></div></div><div></div></div><div></div></div>
<div>
(This article was first published on <strong><a href="http://mazamascience.com/WorkingWithData/?p=1755"> R – Working With Data</a></strong>, and kindly contributed to <a rel="nofollow" href="https://www.r-bloggers.com/">R-bloggers)</a> </div>
<p>Mazama Science has just released the <a target="_blank" rel="nofollow" href="https://cran.r-project.org/package=PWFSLSmoke">PWFSLSmoke</a> package. Source code is available on <a target="_blank" rel="nofollow" href="https://github.com/MazamaScience/PWFSLSmoke">GitHub</a>. Here is the package description:</p>
<blockquote><p>Utilities for working with air quality monitoring data with a focus on small particulates (PM2.5) generated by wildfire smoke. Functions are provided for downloading available data from the United States Environmental Protection Agency (US EPA) and it’s AirNow air quality site. Additional sources of PM2.5 data made accessible by the package include: AIRSIS (password protected), the Western Regional Climate Center (WRCC) and the open source site OpenAQ.</p>
</blockquote>
<p>In this post we discuss the reasons for creating this package and provide examples of its use.<span></span></p>
<p>For several years now, Mazama Science has been working with the <a target="_blank" rel="nofollow" href="https://www.airfire.org/">AirFire</a> team at The USFS Pacific Wildland Fire Sciences Lab. The <strong>PWFSLSmoke</strong> R package is being developed for PWFSL to help modelers and scientists more easily work with PM2.5 data from monitoring locations across North America. The package makes it easier to obtain data, perform analyses and generate reports. It includes functionality to:</p>
<ul>
<li>download and easily work with regulatory PM2.5 data from AirNow</li>
<li>download and quality control raw monitoring data from other sources</li>
<li>convert between UTC and local timezones</li>
<li>apply various algorithms to the data (nowcast, rolling means, aggregation, <em>etc</em>.)</li>
<li>provide interactive timeseries and maps through RStudio’s Viewer pane</li>
<li>create a variety of publication ready maps and timeseries plots</li>
</ul>
<h2>Data Model</h2>
<p>The <strong>PWFSLSmoke</strong> package is designed to provide a compact, full-featured suite of utilities for working with PM 2.5 data used to monitor wildfire smoke. A uniform data model provides consistent data access across monitoring data available from different agencies. The core data model in this package is defined by the <em>ws_monitor</em> object used to store data associated with groups of individual monitors.</p>
<p>To work efficiently with the package it is important to understand the structure of this data object and which functions operate on it. Package functions that begin with <code>monitor_</code>, expect objects of class <em>ws_monitor</em> as their first argument. (<em>Note that ‘ws_’ stands for ‘wildfire smoke’.</em>)</p>
<p>Monitoring data will typically be obtained from an agency charged with archiving data acquired at monitoring sites. For wildifre smoke, the primary pollutant is PM 2.5 and the sites archiving this data include <a target="_blank" rel="nofollow" href="https://www.airnow.gov/">AirNow</a>, <a target="_blank" rel="nofollow" href="http://www.wrcc.dri.edu/">WRCC</a> and <a target="_blank" rel="nofollow" href="https://app.airsis.com/USFS">AIRSIS</a>.</p>
<p>The data model for monitoring data consists of an <strong>R</strong> <code>list</code> with two dataframes: <code>data</code> and <code>meta</code>.</p>
<p>The <code>data</code> dataframe contains all hourly measurements organized with rows (the ‘unlimited’ dimension) as unique timesteps and columns as unique monitors. The very first column is always named <code>datetime</code> and contains the <code>POSIXct</code> datetime in Coordinated Universal Time (UTC).</p>
<p>The <code>meta</code> dataframe contains all metadata associated with monitoring sites and is organized with rows as unique sites and columns as site attributes. The following columns are guaranteed to exist in the <code>meta</code>dataframe:</p>
<ul>
<li><code>monitorID</code> – unique ID for each monitoring site (<em>i.e.</em> each instrument deployment)</li>
<li><code>longitude</code> – decimal degrees East</li>
<li><code>latitude</code> – decimal degrees North</li>
<li><code>elevation</code> – meters above sea level</li>
<li><code>timezone</code> – Olson timezone</li>
<li><code>countryCode</code> – ISO 3166-1 alpha-2 code</li>
<li><code>stateCode</code> – ISO 3166-2 alpha-2 code</li>
</ul>
<p><em>(The <a target="_blank" rel="nofollow" href="https://github.com/MazamaScience/MazamaSpatialUtils">MazamaSpatialUtils</a> package is used to assign timezones and state and country codes.)</em></p>
<p>Additional columns may be available in the <code>meta</code> dataframe and these will depend on the source of the data.</p>
<p>It is important to note that the <code>monitorID</code> acts as a unique key that connects data with metadata. The <code>monitorID</code> is used for column names in the <code>data</code> dataframe and for row names in the <code>meta</code> dataframe. So the following will always be true:</p>
<pre>rownames(ws_monitor$meta) == ws_monitor$meta$monitorID
colnames(ws_monitor$data) == c('datetime',ws_monitor$meta$monitorID)</pre>

<h2>File Size Reduction</h2>
<p>Most providers of monitoring data store the data as it comes in, retaining station metadata with every single datum. However, individual monitors rarely move. Thus, replication of metadata vastly inflates the size of data in this “as-it-comes-in” data management strategy. The <em>ws_monitor</em> data model reflects the reality of monitors that rarely move.</p>
<p>Our <a target="_blank" rel="nofollow" href="http://mazamascience.com/Presentations/Big_Data.pdf">BigData</a> presentation describes how we used this data model to convert unwieldy EPA annual datasets into bite sized .RData files without losing any information. In one example we reduced a 665 Mb ASCII CSV file into a 3.3 Mb .RData file — 1/200<sup>th</sup> the size of the original.</p>
<p>Reduced file sizes make working with data from multiple monitors much easier and the remainder of this post demonstrates the package capabilities by looking at smoke generated from the 2015 wildfire season in the Pacific Northwest</p>
<h2>Pacific Northwest 2015 Wildfires</h2>
<p>In the summer of 2015 Washington state had several catastrophic wildfires that led to heavy smoke in eastern Washington and northern Idaho for quite a few days. We will show how the mapping and timeseries plotting functions in the <strong>PWFSLSmoke</strong> package can help us visualize the spatial and temporal extent of wildfire smoke during these events.</p>
<p>To begin, let’s have a broader look at AirNow ambient monitoring data for the Pacific Northwest – Washington, Oregon and Idaho – from June 1 through October 31, 2015. First, we create a 24-hr rolling mean for each monitor:</p>
<pre>library(PWFSLSmoke)
PacNW &lt;- Northwest_Megafires
# To work with AirNow data directly, uncomment the next two lines
#N_M &lt;- airnow_load(startdate=20150531, enddate=20151101)
#PacNW &lt;- monitor_subset(airnow, stateCodes=c(&quot;WA&quot;, &quot;OR&quot;, &quot;ID&quot;))
PacNW_24 &lt;- monitor_rollingMean(PacNW, width=24)</pre>
<p>Now we can create a map where each monitor is color coded by the maximum value of this 24-hr rolling mean. By default, AQI colors and labels are used but arguments to <code>monitorMap()</code> and <code>addAQILegend()</code> allow users to specify their own.</p>
<pre>monitorMap(PacNW_24, slice=max)
addAQILegend(title=&quot;Max AQI&quot;, cex=0.7)</pre>
<p><a target="_blank" rel="nofollow" href="https://i1.wp.com/mazamascience.com/WorkingWithData/wp-content/uploads/2017/03/monitorMap.png"><img alt="" class="wp-image-1749" src="https://i1.wp.com/mazamascience.com/WorkingWithData/wp-content/uploads/2017/03/monitorMap.png?w=450" data-recalc-dims="1"></a></p>
<p>The map shows that many areas of the Pacific NW had days with unhealthy air but a cluster of sites in Idaho were particularly bad.</p>
<p>(<em><strong>Note</strong> that this is not the regulatory midnight-to-midnight AQI but a continuous 24-hr rolling mean.</em>)</p>
<p>We can us an interactive “leaflet” map to zoom in and get more information:</p>
<pre># Commented out for the blog post
#monitorLeaflet(PacNW_24, slice=max)</pre>
<p>The <code>monitorMap()</code> and <code>monitorLeaflet()</code> plots show us pretty much the same information, except that the leaflet plot allows you to get monitor-specific metadata by clicking on a monitor. In this manner, we can assemble a list of monitorIDs in and around the Nez Perce Reservation in northern Idaho and generate a multi-monitor timeseries plot showing the terrible smoke in late August.</p>
<pre>NezPerceIDs &lt;- c(&quot;160571012&quot;,&quot;160690012&quot;,&quot;160690013&quot;,&quot;160690014&quot;,&quot;160490003&quot;,&quot;160491012&quot;)
NezPerce &lt;- monitor_subset(PacNW, monitorIDs=NezPerceIDs)
monitorPlot_timeseries(NezPerce, style='gnats')</pre>
<p><a target="_blank" rel="nofollow" href="https://i1.wp.com/mazamascience.com/WorkingWithData/wp-content/uploads/2017/03/NezPerce_gnats.png"><img alt="" class="wp-image-1751" src="https://i1.wp.com/mazamascience.com/WorkingWithData/wp-content/uploads/2017/03/NezPerce_gnats.png?w=450" data-recalc-dims="1"></a></p>
<p>At this point it is clear that August is the month of interest so we’ll subset all of our existing <code>ws_monitor</code> objects to cover the month of August with full days according to West Coast time.</p>
<pre>PacNW &lt;- monitor_subset(PacNW, tlim=c(20150801,20150831), timezone=&quot;America/Los_Angeles&quot;)
PacNW_24 &lt;- monitor_subset(PacNW_24, tlim=c(20150801,20150831), timezone=&quot;America/Los_Angeles&quot;)
NezPerce &lt;- monitor_subset(NezPerce, tlim=c(20150801,20150831), timezone=&quot;America/Los_Angeles&quot;)</pre>
<p>We can use the <code>monitorPlot_dailyBarplot()</code> function to look at official, midnight-to-midnight AQI levels for each monitor during the month of August:</p>
<pre>layout(matrix(seq(6)))
par(mar=c(1,1,1,1))
for (monitorID in NezPerceIDs) {
  siteName &lt;- NezPerce$meta[monitorID,'siteName']
  monitorPlot_dailyBarplot(NezPerce, monitorID=monitorID, main=siteName, axes=FALSE) 
}
par(mar=c(5,4,4,2)+.1)
layout(1)</pre>
<p><a target="_blank" rel="nofollow" href="https://i2.wp.com/mazamascience.com/WorkingWithData/wp-content/uploads/2017/03/NezPerce_daily.png"><img alt="" class="wp-image-1750" src="https://i2.wp.com/mazamascience.com/WorkingWithData/wp-content/uploads/2017/03/NezPerce_daily.png?w=450" data-recalc-dims="1"></a></p>
<p>We could also take a more automated approach and directly calculate the location with the worst acute smoke (worst hourly value) during this time period:</p>
<pre>data &lt;- PacNW$data[,-1] # omit 'datetime' column
maxPM25 &lt;- apply(data, 2, max, na.rm=TRUE)
worstAcute &lt;- names(sort(maxPM25, decreasing=TRUE))[1:6]
intersect(worstAcute, NezPerceIDs)
## [1] &quot;160571012&quot; &quot;160490003&quot; &quot;160491012&quot;
PacNW$meta[worstAcute[1],c('siteName','countyName','stateCode')]
##                      siteName countyName stateCode
## 160571012 Juliaetta E-sampler      LATAH        ID</pre>
<p>We see that three of the sites we identified “by hand” also had the worst acute smoke of any sites in the Pacific NW. The site with the highest measured value of PM 2.5 for the month of august was Julietta in Latah county, Idaho.</p>
<p>Let’s do a similar analysis for chronic smoke. In this case we will calculate the number of days at or above the AQI “unhealthy” level:</p>
<pre>PacNW_dailyAvg &lt;- monitor_dailyStatistic(PacNW, FUN=mean, minHours=20)
## Warning in monitor_dailyStatistic(PacNW, FUN = mean, minHours = 20): Found
## 2 timezones. Only the first will be used
data &lt;- PacNW_dailyAvg$data[,-1]
unhealthyDays &lt;- apply(data, 2, function(x){ sum(x >= AQI$breaks_24[4], na.rm=TRUE) })
worstChronic &lt;- names(sort(unhealthyDays, decreasing=TRUE))[1:6]
intersect(worstChronic, NezPerceIDs)
## [1] &quot;160490003&quot; &quot;160571012&quot; &quot;160690014&quot;
PacNW$meta[worstChronic[1],c('siteName','countyName','stateCode')]
##            siteName countyName stateCode
## 160490003 Kamiah-ID      IDAHO        ID</pre>
<p>The area around the Nez Perce Reservation again had three of the worst sites for chronic smoke with the worst being Kamiah in Idaho county, Idaho.</p>
<p>Latah and Idaho counties were unfortunately downwind of several extraordinarily large <a target="_blank" rel="nofollow" href="https://en.wikipedia.org/wiki/2015_Washington_wildfires">Washington state wildfires in 2015</a> with the following ignition dates:</p>
<ul>
<li>Aug 11 – Kettle Complex</li>
<li>Aug 13 – Grizzly Bear Complex</li>
<li>Aug 13 – North Star</li>
<li>Aug 14 – Chelan Complex</li>
<li>Aug 15 – Okanogan Complex</li>
</ul>
<p>A little googling and we can obtain a set of coordinates for these fires to use in a new “Google” map of the August daily average maxima with fire icons at fire locations:</p>
<pre>fireLons &lt;- c(-118.461,-117.679,-120.039,-119.002,-119.662)
fireLats &lt;- c(48.756,46.11,47.814,48.338,48.519)
gmap &lt;- monitorGoogleMap(PacNW_dailyAvg, zoom=7, centerLon=-118, centerLat=47, slice=max)
addIcon('redFlame', fireLons, fireLats, map=gmap, expansion=0.2)
addAQILegend(cex=0.8)
title(&quot;August, 2015&quot;, line=-1.5, cex.main=2)</pre>
<p><a target="_blank" rel="nofollow" href="https://i1.wp.com/mazamascience.com/WorkingWithData/wp-content/uploads/2017/03/WA_megafires.png"><img alt="" class="wp-image-1754" src="https://i1.wp.com/mazamascience.com/WorkingWithData/wp-content/uploads/2017/03/WA_megafires.png?w=450" data-recalc-dims="1"></a></p>
<p>From this map we can see that the monitor in Omak on the Colleville reservation also registered a daily AQI level of “extreme” as it was in-between the two largest fires. The overall spatial pattern, however, shows that the worst impact from these large fires resulted from smoke drifting SE with the prevailing winds and bunching up in the valleys and plains just upwind of the Bitterroot mountains of northern Idaho.</p>
<p>We will now take a closer look at two tribal monitors: Omak for the Colleville tribe and Kamiah for the Nez Perce.</p>
<pre>Omak &lt;- monitor_subset(PacNW, monitorIDs=&quot;530470013&quot;)
Kamiah &lt;- monitor_subset(PacNW, monitorIDs=&quot;160490003&quot;)
layout(matrix(seq(2)))
monitorPlot_dailyBarplot(Omak, main=&quot;August 2015 Daiy AQI -- Omak, WA&quot;,
                         labels_x_nudge=0.8, labels_y_nudge=250)
monitorPlot_dailyBarplot(Kamiah, main=&quot;August Daily AQI -- Kamiah, ID&quot;,
                         labels_x_nudge=0.8, labels_y_nudge=250)
layout(1)</pre>
<p><a target="_blank" rel="nofollow" href="https://i1.wp.com/mazamascience.com/WorkingWithData/wp-content/uploads/2017/03/Omak_Kamiah_daily.png"><img alt="" class="wp-image-1752" src="https://i1.wp.com/mazamascience.com/WorkingWithData/wp-content/uploads/2017/03/Omak_Kamiah_daily.png?w=450" data-recalc-dims="1"></a></p>
<p>We can also examine the diurnal cycle during the very worst days:</p>
<pre>layout(matrix(seq(2)))
par(mar=c(3,4,4,2))
monitorPlot_timeOfDaySpaghetti(Omak, title=&quot;Aug 23-26 Diurnal Smoke -- Omak, WA&quot;,
                               xlab='', ylab='',
                               tlim=c(20150823,20150826))
monitorPlot_timeOfDaySpaghetti(Kamiah, title=&quot;Aug 23-26 Diurnal Smoke -- Kamiah, ID&quot;,
                               xlab='', ylab='',
                               tlim=c(20150823,20150826))
par(mar=c(5,4,4,2)+.1)
layout(1)</pre>
<p><a target="_blank" rel="nofollow" href="https://i0.wp.com/mazamascience.com/WorkingWithData/wp-content/uploads/2017/03/Omak_Kamiah_diurnal.png"><img alt="" class="wp-image-1753" src="https://i0.wp.com/mazamascience.com/WorkingWithData/wp-content/uploads/2017/03/Omak_Kamiah_diurnal.png?w=450" data-recalc-dims="1"></a></p>
<p>While both sites had horrible air all day, in Omak it was somewhat less horrible in the evenings. This is in contrast to Kamiah where the least horrible conditions were encountered in the mornings.</p>
<p>This ends the spatio-temporal exploration of smoke from Pacific NW mega-fire in 2015. We hope this inspires you to harness the mapping plotting functionality available in the <strong>PWFSLSmoke</strong> package.</p>
<p><em>Best Wishes for Healthy Air!</em></p>
<p> </p>
<p> </p>
<p>
<div>
<div>To <strong>leave a comment</strong> for the author, please follow the link and comment on their blog: <strong><a href="http://mazamascience.com/WorkingWithData/?p=1755"> R – Working With Data</a></strong>.</div>
<hr>
<a rel="nofollow" href="https://www.r-bloggers.com/">R-bloggers.com</a> offers <strong><a rel="nofollow" href="https://feedburner.google.com/fb/a/mailverify?uri=RBloggers">daily e-mail updates</a></strong> about <a title="The R Project for Statistical Computing" rel="nofollow" href="https://www.r-project.org/">R</a> news and <a title="R tutorials" rel="nofollow" href="https://www.r-bloggers.com/search/tutorial">tutorials</a> on topics such as: <a title="Data science" rel="nofollow" href="https://www.r-bloggers.com/search/data%20science">Data science</a>, <a title="Big Data" rel="nofollow" href="https://www.r-bloggers.com/search/Big%20Data">Big Data, <a title="R jobs" rel="nofollow" href="https://www.r-users.com/">R jobs</a>, visualization (<a title="ggplot and ggplot2 tutorials" rel="nofollow" href="https://www.r-bloggers.com/search/ggplot2">ggplot2</a>, <a title="Boxplots using lattice and ggplot2 tutorials" rel="nofollow" href="https://www.r-bloggers.com/search/boxplot">Boxplots</a>, <a title="Maps and gis" rel="nofollow" href="https://www.r-bloggers.com/search/map">maps</a>, <a title="Animation in R" rel="nofollow" href="https://www.r-bloggers.com/search/animation">animation</a>), programming (<a title="RStudio IDE for R" rel="nofollow" href="https://www.r-bloggers.com/search/RStudio">RStudio</a>, <a title="Sweave and literate programming" rel="nofollow" href="https://www.r-bloggers.com/search/sweave">Sweave</a>, <a title="LaTeX in R" rel="nofollow" href="https://www.r-bloggers.com/search/LaTeX">LaTeX</a>, <a title="SQL and databases" rel="nofollow" href="https://www.r-bloggers.com/search/SQL">SQL</a>, <a title="Eclipse IDE for R" rel="nofollow" href="https://www.r-bloggers.com/search/eclipse">Eclipse</a>, <a title="git and github, Version Control System" rel="nofollow" href="https://www.r-bloggers.com/search/git">git</a>, <a title="Large data in R using Hadoop" rel="nofollow" href="https://www.r-bloggers.com/search/hadoop">hadoop</a>, <a title="Web Scraping of google, facebook, yahoo, twitter and more using R" rel="nofollow" href="https://www.r-bloggers.com/search/Web+Scraping">Web Scraping</a>) statistics (<a title="Regressions and ANOVA analysis tutorials" rel="nofollow" href="https://www.r-bloggers.com/search/regression">regression</a>, <a title="principal component analysis tutorial" rel="nofollow" href="https://www.r-bloggers.com/search/PCA">PCA</a>, <a title="Time series" rel="nofollow" href="https://www.r-bloggers.com/search/time+series">time series</a>, <a title="finance trading" rel="nofollow" href="https://www.r-bloggers.com/search/trading">trading</a>) and more...
<img height="1" alt="" width="1" src="http://feeds.feedburner.com/~r/RBloggers/~4/itg3_iWtOdM"></a></div>