<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/scada-spikes-in-water-treatment-data/"></div></div><div></div></div><div></div></div>
<div>
(This article was first published on <strong><a href="http://r.prevos.net/scada-spikes/"> The Devil is in the Data</a></strong>, and kindly contributed to <a rel="nofollow" href="https://www.r-bloggers.com/">R-bloggers)</a> </div>
<p><img sizes="(max-width: 300px) 100vw, 300px" data-image-description="" data-orig-file="https://i2.wp.com/r.prevos.net/wp-content/uploads/sites/11/2017/01/dirty-water.jpg?fit=400%2C300" data-orig-size="400,300" data-image-title="dirty-water" class="wp-image-274" data-comments-opened="1" data-recalc-dims="1" data-attachment-id="274" data-permalink="http://r.prevos.net/scada-spikes/dirty-water/" alt="Turbid water" data-medium-file="https://i2.wp.com/r.prevos.net/wp-content/uploads/sites/11/2017/01/dirty-water.jpg?fit=300%2C225" data-large-file="https://i2.wp.com/r.prevos.net/wp-content/uploads/sites/11/2017/01/dirty-water.jpg?fit=400%2C300" data-image-meta="{&quot;aperture&quot;:&quot;0&quot;,&quot;credit&quot;:&quot;&quot;,&quot;camera&quot;:&quot;&quot;,&quot;caption&quot;:&quot;&quot;,&quot;created_timestamp&quot;:&quot;0&quot;,&quot;copyright&quot;:&quot;&quot;,&quot;focal_length&quot;:&quot;0&quot;,&quot;iso&quot;:&quot;0&quot;,&quot;shutter_speed&quot;:&quot;0&quot;,&quot;title&quot;:&quot;&quot;,&quot;orientation&quot;:&quot;0&quot;}" src="https://i2.wp.com/r.prevos.net/wp-content/uploads/sites/11/2017/01/dirty-water.jpg?resize=300%2C228">SCADA spikes are events in the data stream of water treatment plants or similar installations. These SCADA spikes can indicate problems with the process and could result in an increased risk to public health.</p>
<p>The <a target="_blank" rel="nofollow" href="https://www.wsaa.asn.au/publication/health-based-targets-manual">WSAA Health Based Targets Manual </a>specifies a series of decision rules to assess the performance of filtration processes. For example, this rule assesses the performance of conventional filtration:</p>
<blockquote>
<p>“Individual filter turbidity ≤ 0.2 NTU for 95% of month and not > 0.5 NTU for ≥ 15 consecutive minutes.”</p>
</blockquote>
<p><a target="_blank" rel="nofollow" href="https://en.wikipedia.org/wiki/Turbidity">Turbidity</a> is a measure for the cloudiness of a fluid because of large numbers of individual particles otherwise invisible to the naked eye. Turbidity is an important parameter in water treatment because a high level of cloudiness strongly correlates with the presence of microbes. This article shows how to implement this specific decision rule using the R language.</p>
<h2>Simulation</h2>
<p>To create a minimum working example, I first create a simulated SCADA feed for turbidity. The <em>turbidity</em> data frame contains 24 hours of data. The <a target="_blank" rel="nofollow" href="https://stat.ethz.ch/R-manual/R-devel/library/base/html/seq.POSIXt.html"><em>seq.POSIXt</em></a> function creates 24 hours of timestamps at a one-minute spacing. In addition, the <a target="_blank" rel="nofollow" href="https://stat.ethz.ch/R-manual/R-devel/library/stats/html/Normal.html"><em>rnorm</em></a> function creates 1440 turbidity readings with an average of 0.1 NTU and a standard deviation of 0.01 NTU. The image below visualises the simulated data. The next step is to assess this data in accordance with the decision rule.</p>
<pre>
# Simulate data
set.seed(1234)
turbidity &lt;- data.frame(DateTime = seq.POSIXt(as.POSIXct(&quot;2017-01-01 00:00:00&quot;), by = &quot;min&quot;, length.out=24*60),
                        Turbidity = rnorm(n = 24*60, mean = 0.1, sd = 0.01)
 )
</pre>
<p>The second section simulates five spikes in the data. The first line picks a random start time for the spike. The second line in the for-loop picks a duration between 10 and 30 minutes. In addition, the third line simulates the value of the spike. The mean value of the spike is determined by the <em><a target="_blank" rel="nofollow" href="http://stat.ethz.ch/R-manual/R-devel/library/stats/html/Binomial.html">rbinom</a></em> function to create either a low or a high spike. The remainder of the spike simulation inserts the new data into the turbidity data frame.</p>
<pre>
# Simulate spikes
for (i in 1:5) {
   time &lt;- sample(turbidity$DateTime, 1)
   duration &lt;- sample(10:30, 1)
   value &lt;- rnorm(1, 0.5 * rbinom(1, 1, 0.5) + 0.3, 0.05)
   start &lt;- which(turbidity$DateTime == time)
   turbidity$Turbidity[start:(start+duration - 1)] &lt;- rnorm(duration, value, value/10)
}
</pre>
<p>The image below visualises the simulated data using the mighty <a target="_blank" rel="nofollow" href="http://ggplot2.org/">ggplot</a>. Only four spikes are visible because two of them overlap. The next step is to assess this data in accordance with the decision rule.</p>
<pre>library(ggplot2)
ggplot(turbidity, aes(x = DateTime, y = Turbidity)) + geom_line(size = 0.2) + 
   geom_hline(yintercept = 0.5, col = &quot;red&quot;) + ylim(0,max(turbidity$Turbidity)) + 
   ggtitle(&quot;Simulated SCADA data&quot;)
</pre>
<p><img sizes="(max-width: 584px) 100vw, 584px" data-image-description="&lt;p&gt;Simulated SCADA data with spikes&lt;/p&gt;
" data-orig-file="https://i0.wp.com/r.prevos.net/wp-content/uploads/sites/11/2017/02/spikes.png?fit=1718%2C962" data-orig-size="1718,962" data-image-title="scada spikes" data-comments-opened="1" class="wp-image-262" data-recalc-dims="1" data-attachment-id="262" title="Simulated SCADA data with spikes" data-permalink="http://r.prevos.net/scada-spikes/spikes/" alt="Simulated SCADA data with spikes" data-medium-file="https://i0.wp.com/r.prevos.net/wp-content/uploads/sites/11/2017/02/spikes.png?fit=300%2C168" data-large-file="https://i0.wp.com/r.prevos.net/wp-content/uploads/sites/11/2017/02/spikes.png?fit=584%2C327" data-image-meta="{&quot;aperture&quot;:&quot;0&quot;,&quot;credit&quot;:&quot;&quot;,&quot;camera&quot;:&quot;&quot;,&quot;caption&quot;:&quot;&quot;,&quot;created_timestamp&quot;:&quot;0&quot;,&quot;copyright&quot;:&quot;&quot;,&quot;focal_length&quot;:&quot;0&quot;,&quot;iso&quot;:&quot;0&quot;,&quot;shutter_speed&quot;:&quot;0&quot;,&quot;title&quot;:&quot;&quot;,&quot;orientation&quot;:&quot;0&quot;}" src="https://i0.wp.com/r.prevos.net/wp-content/uploads/sites/11/2017/02/spikes.png?resize=584%2C327"></p>
<h2>SCADA Spikes Detection</h2>
<p>The following code searches for all spikes over 0.50 NTU using the <a target="_blank" rel="nofollow" href="https://stat.ethz.ch/R-manual/R-devel/library/base/html/rle.html">run length function</a>. This function transforms a vector into a vector of values and lengths. For example, the run length of the vector c(1, 1, 2, 2, 2, 3, 3, 3, 3, 5, 5, 6) is:</p>
<ul>
<li>lengths: int [1:5] 2 3 4 2 1</li>
<li>values : num [1:5] 1 2 3 5 6</li>
</ul>
<p>The value 1 has a length of 1, the value 2 has a length of 3 and so on. The spike detection code creates the run length for turbidity levels greater than 0.5, which results in a boolean vector. The <a target="_blank" rel="nofollow" href="http://stat.ethz.ch/R-manual/R-devel/library/base/html/cumsum.html">cumsum</a> function calculates the starting point of each spike which allows us to calculate their duration.</p>
<p>The code results in a data frame with all spikes higher than 0.50 NTU and longer than 15 minutes. The spike that occurred at 11:29 was higher than 0.50 NTU and lasted for 24 minutes. The other three spikes are either lower than 0.50 NTU. The first high spike lasted less than 15 minutes.</p>
<pre>
# Spike Detection
spike.detect &lt;- function(DateTime, Value, Height, Duration) {
 runlength &lt;- rle(Value > Height)
 spikes &lt;- data.frame(Spike = runlength$values,
 times &lt;- cumsum(runlength$lengths))
 spikes$Times &lt;- DateTime[spikes$times]
 spikes$Event &lt;- c(0,spikes$Times[-1] - spikes$Times[-nrow(spikes)])
 spikes &lt;- subset(spikes, Spike == TRUE &amp; Event > Duration)
 return(spikes)
}
spike.detect(turbidity$DateTime, turbidity$Turbidity, 0.5, 15)
</pre>
<p>This approach was used to prototype a <a target="_blank" rel="nofollow" href="http://digitaledition.awa.asn.au/?iid=131096&startpage=page0000064#folio=64">software package to assess water treatment plant data in accordance with the Health-Based Targets Manual</a>. The finished product has been written in SQL and is available under an Open Source sharing license.</p>
<p>The post <a target="_blank" rel="nofollow" href="http://r.prevos.net/scada-spikes/">SCADA spikes in Water Treatment Data</a> appeared first on <a target="_blank" rel="nofollow" href="http://r.prevos.net/">The Devil is in the Data</a>.</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://r.prevos.net/scada-spikes/"> The Devil is in the 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/9-FK2AB8Kdc"></a></div>