<p><div><div><div><div data-show-faces="false" data-height="21" data-width="100" data-layout="button_count" data-send="true" data-href="https://www.r-bloggers.com/smoke-from-a-distant-fire/"></div></div><div></div></div><div></div></div>
<div>
(This article was first published on <strong><a href="http://bayesball.blogspot.com/2018/08/smoke-from-distant-fire.html"> Bayes Ball</a></strong>, and kindly contributed to <a rel="nofollow" href="https://www.r-bloggers.com/">R-bloggers)</a> </div>
<h4>Forest fires and air quality</h4>
<h4><em>August 31, 2018</em></h4>
<p>It was recently announced that during 2018, <a rel="nofollow" href="https://www.cbc.ca/news/canada/british-columbia/state-emergency-bc-wildfires-1.4803546" target="_blank">British Columbia has seen the most extensive forest fire season on record</a>. As I write this (2018-08-31) there are <a rel="nofollow" href="http://bcfireinfo.for.gov.bc.ca/hprScripts/WildfireNews/Fires.asp?Mode=normal&AllFires=0&FC=0" target="_blank">currently 442 wildfires burning in British Columbia</a>. These fires have a significant impact on people's lives - many areas are under evacuation order and evacuation alert, and there are reports that homes have been destroyed by the blazes.</p>
<p>The fires also create a significant amount of smoke, which has been pushed great distances by the shifting winds. This includes the large population centres of Vancouver and Victoria in British Columbia, as well as the Seattle metropolitan region and elsewhere in Washington. (Clifford Mass, Professor of Atmospheric Sciences at the Universtiy of Washington in Seattle, has written extensively about the smoke events in the region; see for example <a rel="nofollow" href="http://cliffmass.blogspot.com/2018/08/western-washington-smoke-darkest-before.html" target="_blank">Western Washington Smoke: Darkest Before the Dawn</a> from 2018-08-22.)</p>
<p>The Province of British Columbia has many air quality monitoring stations around the province, and makes the data available. The measure most used for monitoring the effects on health is PM25 or PM2.5, for fine particles with a diameter of 2.5 microns (millionths of a metre). The B.C. government has a <a rel="nofollow" href="http://www.env.gov.bc.ca/epd/bcairquality/readings/find-stations-map-PM25.html" target="_blank">Current Particulate Matter map</a> that colour codes the one hour average measures for all the testing stations around the province.</p>
<h3>The data file and a simple plot</h3>
<p>The DataBC Catalogue provides access to air quality data. There's <a rel="nofollow" href="https://catalogue.data.gov.bc.ca/dataset/air-quality-monitoring-verified-hourly-data" target="_blank">&quot;verifiedâ€� to the end of 2017</a>, and <a rel="nofollow" href="https://catalogue.data.gov.bc.ca/dataset/air-quality-monitoring-unverified-hourly-air-quality-and-meteorological-data" target="_blank">&quot;unverifiedâ€� for the past 30 days</a>. Since we want to see what happened this month, it's the latter we want. (The page with the <a rel="nofollow" target="_blank">links to the raw files is here</a>.)</p>
<p>The files are arranged by particulate or gas type; there's a table for ozone and another for sulpher dioxide, and others for the particulate matter. Note that the data are made available under the Province of B.C.'s Open Data license, and are in nice tidy form. And the date format is <a rel="nofollow" href="https://en.wikipedia.org/wiki/ISO_8601" target="_blank">ISO 8601</a>, which makes me happy.</p>
<p>To make sure we've got a reproducible version, I've saved the file I downloaded early this morning to my google drive. The link to the folder is <a rel="nofollow" href="https://drive.google.com/drive/folders/1zxhz9niF5ooLHw-_JCylCyRVbH8wPkfK?usp=sharing" target="_blank">here</a>.</p>
<p>For the first plot, let's look at the PM2.5 level for my hometown of Victoria, B.C. The code below loads the R packages we’ll use, reads the data, and generates the plot.</p>
<p></p>
<pre><code># tidyverse packages<br>library(tidyverse)<br>library(glue)<br></code></pre>
<pre><code><pre><code><br></code></pre>
<pre><code>PM25_data &lt;- readr::read_csv(&quot;PM25_2018-08-31.csv&quot;)</code></pre>
<div></div>
<p>filter(STATION_NAME == &quot;Victoria Topaz&quot;) %>% ggplot() + geom_line(aes(x = DATE_PST, y = REPORTED_VALUE)) + labs(x = &quot;date&quot;, title = glue(&quot;Air quality: Victoria Topaz&quot;), subtitle = &quot;one hour average, Âµg/m3 of PM2.5&quot;, caption = &quot;data: B.C. Ministry of Environment and Climate Change Strategy&quot;) </p>
<div><a rel="nofollow" href="https://i1.wp.com/3.bp.blogspot.com/-Kp9IY2x_4C4/W4lIYC4khfI/AAAAAAAAJQA/k-UuKU3wHAIu-7SXIoSrkU-7JtoSi-npwCLcBGAs/s1600/Victoria-1.png?ssl=1" target="_blank"><img border="0" src="https://i0.wp.com/3.bp.blogspot.com/-Kp9IY2x_4C4/W4lIYC4khfI/AAAAAAAAJQA/k-UuKU3wHAIu-7SXIoSrkU-7JtoSi-npwCLcBGAs/s320/Victoria-1.png?resize=450%2C961&ssl=1" data-original-width="450" width="450" data-original-height="961" data-recalc-dims="1"></a></div>
<p>There are 61 air quality monitoring stations around British Columbia. It would be interesting to see how the air quality was in other parts of the region - and since over half (54% in 2017) of the province's population lives in the Vancouver Census Metropolitan Area (CMA), let's plot the air quality there. There are multiple stations in the Vancouver CMA, so I chose the one at Burnaby Southâ€¦it's fairly central in the region.We can run this line of code to see a listing of all 61 stations (but we won't do that nowâ€¦)</p>
<pre><code># list all the air quality stations for which there is PM2.5 data<br>unique(PM25_data$STATION_NAME)</code></pre>
<p>And since we're going to be doing this often, let's wrap the code that filters for the location we want and runs the plot in a function. Note that we'll create a new variable <code>station_name</code> so all we need to do to change the plot is assign the name of the station we want, and off we go. Not only does this simplify our lives now, but is all-but-essential for a Shiny application.</p>
<pre><code><br></code></pre>
<pre><code># the air quality plot<br>PM25_plot &lt;- function(datafile, station_name){<br>  datafile %>%<br>  filter(STATION_NAME == station_name) %>%<br>  ggplot() +<br>  geom_line(aes(x = DATE_PST, y = REPORTED_VALUE)) +<br>  labs(x = &quot;date&quot;,<br>       title = glue(&quot;Air quality: &quot;, station_name),<br>       subtitle = &quot;one hour average, Âµg/m3 of PM2.5&quot;, <br>       caption = &quot;data: B.C. Ministry of Environment and Climate Change Strategy&quot;)<br>}</code></pre>
<p>Now that we’ve got the function, the code to create the plot for the Burnaby South station is significantly simplified: assign the station name, and call the function.</p>
<pre><code># our Burnaby plot<br>station_name &lt;- &quot;Burnaby South&quot;<br><br>PM25_plot(PM25_data, station_name)</code></pre>
<div><a rel="nofollow" href="https://i0.wp.com/2.bp.blogspot.com/-DoRTZyhw5Oc/W4lIkBQHe2I/AAAAAAAAJQE/gZbUHLPlBDcCrgEhfkBwFKGQ_RCNvUI5ACLcBGAs/s1600/Burnaby-1.png?ssl=1" target="_blank"><img border="0" src="https://i0.wp.com/2.bp.blogspot.com/-DoRTZyhw5Oc/W4lIkBQHe2I/AAAAAAAAJQE/gZbUHLPlBDcCrgEhfkBwFKGQ_RCNvUI5ACLcBGAs/s320/Burnaby-1.png?resize=450%2C961&ssl=1" data-original-width="450" width="450" data-original-height="961" data-recalc-dims="1"></a></div>
<p>And what about the towns that are the closest to the fires? While there are fires burning across the province, the fires that are burning the forests of the Nechako Plateau have understandably received a lot of attention. You may have seen the news stories and images from Prince George like <a rel="nofollow" href="https://globalnews.ca/news/4393740/prince-george-wildfire-smoke/" target="_blank">this</a> and <a rel="nofollow" href="https://globalnews.ca/news/4395799/fort-st-james-smoke-timelapse/" target="_blank">this</a>, or the <a rel="nofollow" href="https://worldview.earthdata.nasa.gov/?p=geographic&l=VIIRS_SNPP_CorrectedReflectance_TrueColor(hidden),MODIS_Aqua_CorrectedReflectance_TrueColor(hidden),MODIS_Terra_CorrectedReflectance_TrueColor,Reference_Labels,Reference_Features(hidden),Coastlines&t=2018-08-17-T00%3A00%3A00Z&z=3&v=-129.8449397281824,48.691966947332034,-115.07931472818241,57.067943509832034" target="_blank">images of the smoke plume from the NASA Worldview site</a>.</p>
<div>Prince George is east of many major fires, downwind of the prevailing westerly winds. So what has the air quality in Prince George been like?
<pre><code>station_name &lt;- &quot;Prince George Plaza 400&quot;<br>PM25_plot(PM25_data, station_name)</code></pre>
<p></p>
<div></div>
<p></p>
<div><a rel="nofollow" href="https://i1.wp.com/3.bp.blogspot.com/-rgwmoGNRx5M/W4lIq4oFmPI/AAAAAAAAJQM/ENyfndjIskofYlVxCAhD3roUuQI-D9ZXwCLcBGAs/s1600/PrinceGeorge-1.png?ssl=1" target="_blank"><img border="0" src="https://i0.wp.com/3.bp.blogspot.com/-rgwmoGNRx5M/W4lIq4oFmPI/AAAAAAAAJQM/ENyfndjIskofYlVxCAhD3roUuQI-D9ZXwCLcBGAs/s320/PrinceGeorge-1.png?resize=450%2C961&ssl=1" data-original-width="450" width="450" data-original-height="961" data-recalc-dims="1"></a></div>
<p>Or still closer to the fires, the town of Burns Lake.</p>
<pre><code>station_name &lt;- &quot;Burns Lake Fire Centre&quot;<br>PM25_plot(PM25_data, station_name)</code></pre>
<p></p>
<div><a rel="nofollow" href="https://i2.wp.com/3.bp.blogspot.com/-_fNeBOY8XS4/W4lIqwvIOKI/AAAAAAAAJQI/0Ay7F5PhYjEzrLm-bVyPadIYYEoM9PQNgCEwYBhgL/s1600/BurnsLake-1.png?ssl=1" target="_blank"><img border="0" src="https://i1.wp.com/3.bp.blogspot.com/-_fNeBOY8XS4/W4lIqwvIOKI/AAAAAAAAJQI/0Ay7F5PhYjEzrLm-bVyPadIYYEoM9PQNgCEwYBhgL/s320/BurnsLake-1.png?resize=450%2C961&ssl=1" data-original-width="450" width="450" data-original-height="961" data-recalc-dims="1"></a></div>
<p>The town of Smithers is west of the fires that are burning on the Nechako Plateau and producing all the smoke experienced in Burns Lake and Prince George. The residents of Smithers have had a very different experience, only seeing smoke in the sky when the winds shifted to become easterly.<br><code><br></code><code>station_name &lt;- &quot;Smithers St Josephs&quot; </code><br><code>PM25_plot(PM25_data, station_name)</code><br><code><br></code></p>
<div><a rel="nofollow" href="https://i0.wp.com/4.bp.blogspot.com/-xFN3GIPtqHY/W4lIq3fjhgI/AAAAAAAAJQk/IvFfadiy7AUK_RBeaN8H8ULRf6-6XktagCEwYBhgL/s1600/Smithers-1.png?ssl=1" target="_blank"><img border="0" src="https://i1.wp.com/4.bp.blogspot.com/-xFN3GIPtqHY/W4lIq3fjhgI/AAAAAAAAJQk/IvFfadiy7AUK_RBeaN8H8ULRf6-6XktagCEwYBhgL/s320/Smithers-1.png?resize=450%2C961&ssl=1" data-original-width="450" width="450" data-original-height="961" data-recalc-dims="1"></a></div>
<p><code><br></code></p>
<h3>multiple stations in one plot</h3>
<p>You may have noticed that the Y axis on the plots can be quite different - for example, Victoria reaches 300, Smithers gets to 400, and Prince George is double that at 800, and Burns Lake is more than double again. There are two ways we can compare multiple stations: a single plot, or faceted plots.</p>
<h4>a line plot with four stations</h4>
<pre><code>station_name &lt;- c(&quot;Burns Lake Fire Centre&quot;, &quot;Prince George Plaza 400&quot;, </code></pre>
<pre><code>                  &quot;Smithers St Josephs&quot;, &quot;Victoria Topaz&quot;)<br><br>PM25_data %>%<br>  filter(STATION_NAME %in% station_name) %>%<br>  ggplot() +<br>  geom_line(aes(x = DATE_PST, y = REPORTED_VALUE, colour = STATION_NAME)) +<br>  labs(x = &quot;date&quot;,<br>       title = glue(&quot;Air quality: Burns Lake, Prince George, Smithers, Victoria&quot;),<br>       subtitle = &quot;one hour average, Âµg/m3 of PM2.5&quot;, <br>       caption = &quot;data: Ministry of Environment and Climate Change Strategy&quot;)</code></pre>
<p></p>
<div><a rel="nofollow" href="https://i0.wp.com/2.bp.blogspot.com/-8yylRry3GM0/W4lIsobKIFI/AAAAAAAAJQs/VYOXXn3dNVs3oLTtNRTtI5qFnU6CgJNbgCEwYBhgL/s1600/four_plot-1.png?ssl=1" target="_blank"><img border="0" src="https://i2.wp.com/2.bp.blogspot.com/-8yylRry3GM0/W4lIsobKIFI/AAAAAAAAJQs/VYOXXn3dNVs3oLTtNRTtI5qFnU6CgJNbgCEwYBhgL/s320/four_plot-1.png?resize=450%2C961&ssl=1" data-original-width="450" width="450" data-original-height="961" data-recalc-dims="1"></a></div>
<p>With four complex lines as we have here, it can be hard to discern which line is which. </div>
<div></div>
<div>
<h4>Use facets to plot the four stations separately</h4>
<div></div>
</div>
<div>Facets give us another way to view the comparisons. In the first version, with the facets stacked vertically, it emphasizes comparisons on the X axis - that is, over time. In this way, we can see that the four locations have had smoke events that have occurred at different times.
<pre><code>station_name &lt;- c(&quot;Burns Lake Fire Centre&quot;, &quot;Prince George Plaza 400&quot;, </code></pre>
<pre><code>                  &quot;Smithers St Josephs&quot;, &quot;Victoria Topaz&quot;)<br><br><br>PM25_data %>%<br>  filter(STATION_NAME %in% station_name) %>%<br>  ggplot() +<br>  geom_line(aes(x = DATE_PST, y = REPORTED_VALUE)) +<br>  facet_grid(STATION_NAME ~ .) +<br>  labs(title = glue(&quot;Air quality: Burns Lake, Prince George, Smithers, Victoria&quot;),<br>       subtitle = &quot;one hour average, Âµg/m3 of PM2.5&quot;, <br>       caption = &quot;data: B.C. Ministry of Environment and Climate Change Strategy&quot;) +<br>  theme(axis.text.x=element_text(size=rel(0.75), angle=90),<br>        axis.title = element_blank())</code></pre>
<pre><code><br></code></pre>
<div><a rel="nofollow" href="https://i2.wp.com/3.bp.blogspot.com/-sow8rXAmMqA/W4lIrpTHX3I/AAAAAAAAJQw/2af2X1IIDkUWDh4mAJawZQVQbEw1qaJfQCEwYBhgL/s1600/facets-1.png?ssl=1" target="_blank"><img border="0" src="https://i1.wp.com/3.bp.blogspot.com/-sow8rXAmMqA/W4lIrpTHX3I/AAAAAAAAJQw/2af2X1IIDkUWDh4mAJawZQVQbEw1qaJfQCEwYBhgL/s320/facets-1.png?resize=450%2C961&ssl=1" data-original-width="450" width="450" data-original-height="961" data-recalc-dims="1"></a></div>
<pre><code><br></code></pre>
<pre><code><br></code></pre>
<p>In the second version, the facets are placed horizontally, making comparisons on the Y axis clear. The smoke events in the four locations have been of very different magnitudes. </p>
<pre><code>PM25_data %>%<br>  filter(STATION_NAME %in% station_name) %>%<br>  ggplot() +<br>  geom_line(aes(x = DATE_PST, y = REPORTED_VALUE)) +<br>  facet_grid(. ~ STATION_NAME) +<br>  labs(title = glue(&quot;Air quality: Burns Lake, Prince George, Smithers, Victoria&quot;),<br>       subtitle = &quot;one hour average, Âµg/m3 of PM2.5&quot;, <br>       caption = &quot;data: B.C. Ministry of Environment and Climate Change Strategy&quot;) +<br>  theme(axis.text.x=element_text(size=rel(0.75), angle=90),<br>        axis.title = element_blank())</code></pre>
<div><a rel="nofollow" href="https://i0.wp.com/4.bp.blogspot.com/-7ysvVQ_diQw/W4lIsBlCDkI/AAAAAAAAJQo/LGyqFgdkgR0Z9ODmRkWukjKCpSupYAL-ACEwYBhgL/s1600/facets-2.png?ssl=1" target="_blank"><img border="0" src="https://i0.wp.com/4.bp.blogspot.com/-7ysvVQ_diQw/W4lIsBlCDkI/AAAAAAAAJQo/LGyqFgdkgR0Z9ODmRkWukjKCpSupYAL-ACEwYBhgL/s320/facets-2.png?resize=450%2C961&ssl=1" data-original-width="450" width="450" data-original-height="961" data-recalc-dims="1"></a></div>
<p>These two plots show not only that Burns Lake and Prince George have had the most extreme smoke events, but that they have had sustained periods of poor air quality through the whole month. While the most extreme event in Smithers exceeds that of Victoria, there hasn't been a prolonged period of smoke in the air like the other three locations.</p>
<p>-30-</p>
<p></div>

				
		
<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://bayesball.blogspot.com/2018/08/smoke-from-distant-fire.html"> Bayes Ball</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 rel="nofollow" href="https://www.r-project.org/" title="The R Project for Statistical Computing">R</a> news and <a rel="nofollow" href="https://www.r-bloggers.com/search/tutorial" title="R tutorials">tutorials</a> on topics such as: <a rel="nofollow" href="https://www.r-bloggers.com/search/data%20science" title="Data science">Data science</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/Big%20Data" title="Big Data">Big Data, <a rel="nofollow" href="https://www.r-users.com/" title="R jobs">R jobs</a>, visualization (<a rel="nofollow" href="https://www.r-bloggers.com/search/ggplot2" title="ggplot and ggplot2 tutorials">ggplot2</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/boxplot" title="Boxplots using lattice and ggplot2 tutorials">Boxplots</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/map" title="Maps and gis">maps</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/animation" title="Animation in R">animation</a>), programming (<a rel="nofollow" href="https://www.r-bloggers.com/search/RStudio" title="RStudio IDE for R">RStudio</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/sweave" title="Sweave and literate programming">Sweave</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/LaTeX" title="LaTeX in R">LaTeX</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/SQL" title="SQL and databases">SQL</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/eclipse" title="Eclipse IDE for R">Eclipse</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/git" title="git and github, Version Control System">git</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/hadoop" title="Large data in R using Hadoop">hadoop</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/Web+Scraping" title="Web Scraping of google, facebook, yahoo, twitter and more using R">Web Scraping</a>) statistics (<a rel="nofollow" href="https://www.r-bloggers.com/search/regression" title="Regressions and ANOVA analysis tutorials">regression</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/PCA" title="principal component analysis tutorial">PCA</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/time+series" title="Time series">time series</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/trading" title="finance trading">trading</a>) and more...
<img width="1" alt="" src="http://feeds.feedburner.com/~r/RBloggers/~4/P4E95HRyj2U" height="1"></a></div></code></pre>