<p><div><div><div><div data-width="100" data-show-faces="false" data-send="true" data-layout="button_count" data-height="21" data-href="http://www.r-bloggers.com/german-temperature-data/"></div></div><div><a data-url="http://www.r-bloggers.com/german-temperature-data/" data-counturl="http://www.r-bloggers.com/german-temperature-data/" data-count="horizontal" data-text="German Temperature Data" class="twitter-share-button" data-via="rbloggers" href="https://twitter.com/share"></a></div><div></div></div><div></div></div>
<div>
(This article was first published on <strong><a href="http://opiateforthemass.es/articles/dwd-daten/"> Opiate for the masses</a></strong>, and kindly contributed to <a rel="nofollow" href="http://www.r-bloggers.com/">R-bloggers)</a> </div>
<h2>Temperatures</h2>
<p>Personally, I thrive in cooler weather. In fact, any temperature above 20°C is summer for me, and anything above 25°C is hot. Most people don’t understand me, but I prefer cold over warm. Fact: if it’s cold, you can put on another piece of clothing until you’re feeling warm again. If, however, it’s warm… Well, there’s only a limited amount of clothes you can take off until you’re accused of crimes of moral turpitude. Plus, if it’s above 25°C, even going full monty won’t get you to cool off.</p>
<p>I recently had an discussion if this year’s March was warmer than the average. Generally speaking, with global average temperatures <a target="_blank" href="http://climate.nasa.gov/vital-signs/global-temperature/">steadily increasing since the 1970s</a>, I think it’s safe to say that for a given random month, you’re better off stating that it was warmer than the average compared to colder. Which is exactly what I did. Turns out, I was wrong… But I learnt how to change the colour of a <code>ggplot()</code> line mid-way between two points, based on an arbitrary cutoff value!</p>
<h2>Data</h2>
<p>As a data scientist, I not only blurt out statements, but I am also in the (judging the current political discussions, unique) opportunity to gather facts to check my ejaculations (<em>teehee</em>, he said <a target="_blank" href="http://qi.com/infocloud/sherlock-holmes">“ejaculate”</a>). The official source for German weather data is the <a target="_blank" href="https://en.wikipedia.org/wiki/Deutscher_Wetterdienst">Deutscher Wetterdienst</a>, the DWD (translated as the German Meteorological Service). Since it runs the <a target="_blank" href="http://www.dwd.de/DE/klimaumwelt/cdc/cdc_node.html">national climate archive</a>, it also provides data on temperatures, amongst many other data like precipitation and historical <a target="_blank">vine ripeness data</a> (now, there’s another blogpost just waiting to be written, I’d say!).</p>
<p>We’re only interested in the <a target="_blank">monthly regional air temperature averages</a>. Let’s grab the monthly data for all German regions from the DWD’s FTP:</p>
<figure>
<pre><code data-lang="r"><span>f.get.dwd.data</span> <span>&lt;-</span> <span>function</span><span>()</span> <span>{</span>
  <span># load libraries needed in function
</span>  <span>library</span><span>(</span><span>readr</span><span>)</span>
  <span>library</span><span>(</span><span>tidyr</span><span>)</span>
  
  <span># generate empty list to fill with data in the loop
</span>  <span>raw.dta</span> <span>&lt;-</span> <span>vector</span><span>(</span><span>&quot;list&quot;</span><span>,</span> <span>length</span> <span>=</span> <span>12</span><span>)</span>
  
  <span>for</span><span>(</span><span>i</span> <span>in</span> <span>seq_len</span><span>(</span><span>12</span><span>)){</span>
    <span># grab all 12 months of data from DWD CDC
</span>    <span>url</span> <span>&lt;-</span> <span>paste0</span><span>(</span>
      <span>&quot;ftp://ftp-cdc.dwd.de/pub/CDC/&quot;</span><span>,</span>
      <span>&quot;regional_averages_DE/monthly/air_temperature_mean/&quot;</span><span>,</span>
      <span>&quot;regional_averages_tm_&quot;</span><span>,</span> <span>sprintf</span><span>(</span><span>&quot;%02i&quot;</span><span>,</span> <span>i</span><span>),</span> <span>&quot;.txt&quot;</span>
    <span>)</span>
    
    <span># check max length once
</span>    <span>if</span><span>(</span><span>i</span> <span>==</span> <span>1</span><span>){</span>
      <span>lngth</span> <span>&lt;-</span> <span>length</span><span>(</span><span>read_lines</span><span>(</span><span>url</span><span>))</span>
    <span>}</span>

    <span>f.get.length</span> <span>&lt;-</span> <span>function</span><span>(</span><span>lngth</span><span>)</span> <span>{</span>
      <span># if i &lt; current month, return Januar length
</span>      <span># else, length is one month less than Januar legnth
</span>      <span>if</span><span>(</span><span>i</span> <span>&lt;</span> <span>as.integer</span><span>(</span><span>format</span><span>(</span><span>Sys.Date</span><span>(),</span> <span>format</span> <span>=</span> <span>&quot;%m&quot;</span><span>))){</span>
        <span># -1 for skip, -1 for last line removed
</span>        <span>return</span><span>(</span><span>lngth</span><span>-2</span><span>)</span>
      <span>}</span> <span>else</span> <span>{</span>
        <span>return</span><span>(</span><span>lngth</span><span>-3</span><span>)</span>
      <span>}</span>
    <span>}</span>
    
    <span># read data with readr, since this automatically parses the urls
</span>    <span>txt.dta</span> <span>&lt;-</span> <span>read_lines</span><span>(</span><span>url</span><span>,</span> <span>skip</span> <span>=</span> <span>1</span><span>,</span> <span>n_max</span> <span>=</span> <span>f.get.length</span><span>(</span><span>lngth</span><span>))</span>
    
    <span># add data to list
</span>    <span># do NOT use readr, since this somehow completely messes up the data
</span>    <span>raw.dta</span><span>[[</span><span>i</span><span>]]</span> <span>&lt;-</span> <span>read.csv2</span><span>(</span><span>text</span> <span>=</span> <span>txt.dta</span><span>,</span> <span>stringsAsFactors</span> <span>=</span> <span>FALSE</span><span>,</span> <span>dec</span> <span>=</span> <span>&quot;.&quot;</span><span>)</span> <span>%>%</span> 
      <span># last row is fubar'ed, remove
</span>      <span>select</span><span>(</span><span>-</span><span>X</span><span>)</span>
  <span>}</span>
  
  <span># bind all into one data frame
</span>  <span>dta</span> <span>&lt;-</span> <span>bind_rows</span><span>(</span><span>raw.dta</span><span>)</span> <span>%>%</span> 
    <span># clean data set
</span>    <span>gather</span><span>(</span><span>region</span><span>,</span> <span>temperature</span><span>,</span> <span>-</span><span>Jahr</span><span>,</span> <span>-</span><span>Monat</span><span>)</span> <span>%>%</span> 
    <span># rename year and month to English
</span>    <span>select</span><span>(</span><span>year</span> <span>=</span> <span>Jahr</span><span>,</span> <span>month</span> <span>=</span> <span>Monat</span><span>,</span> <span>everything</span><span>())</span> <span>%>%</span> 
    <span>arrange</span><span>(</span><span>region</span><span>,</span> <span>year</span><span>,</span> <span>month</span><span>)</span>
  
  <span>return</span><span>(</span><span>dta</span><span>)</span>
<span>}</span></code></pre>
</figure>
<p>Pulling the 12 files (one for each month) is pretty straightforward with a <code>for</code> loop. I found <code>read_csv()</code> from <code>library(readr)</code> to choke on non-locale decimals and negative numbers. A lesson learnt: I spent half a day trying to figure out what went wrong with my download script. To “fix” the problem, I read the file as a text file with <code>read_lines()</code> (since <code>readr</code> functions nicely default to downloading files if it’s a web resource), and then feeding this text into standard R’s <code>read.csv()</code> function, specifying the non-standard decimal. I skip the first line (since it’s not data, but an explanation), and the last row (since it’s garbage). Skipping the last row needs some checking of the length of the file, done in <code>f.get.length()</code>.</p>
<p>Once you’ve set up the function, it’s easy to just grab the data and save it for later use:</p>
<figure>
<pre><code data-lang="r"><span>if</span><span>(</span><span>file.exists</span><span>(</span><span>&quot;dwd.Rdata&quot;</span><span>)){</span>
  <span>load</span><span>(</span><span>&quot;dwd.Rdata&quot;</span><span>)</span>
<span>}</span> <span>else</span> <span>{</span>
  <span>dta</span> <span>&lt;-</span> <span>f.get.dwd.data</span><span>()</span>
  <span>save</span><span>(</span><span>dta</span><span>,</span> <span>file</span> <span>=</span> <span>&quot;dwd.Rdata&quot;</span><span>)</span>
<span>}</span></code></pre>
</figure>
<p>All in all, relatively easy to get the data (modulo the <code>readr</code> problem).</p>
<h2>Graphing</h2>
<p>Right, let’s have a look at the data. I live in North Rhine-Westphalia currently, so let’s just look at that data, for March only and generate a rolling 10-year average.</p>
<figure>
<pre><code data-lang="r"><span>nrw</span> <span>&lt;-</span> <span>dta</span> <span>%>%</span>
  <span>filter</span><span>(</span><span>region</span> <span>==</span> <span>&quot;Nordrhein.Westfalen&quot;</span><span>)</span> <span>%>%</span> 
  <span>filter</span><span>(</span><span>month</span> <span>==</span> <span>3</span><span>)</span> <span>%>%</span> 
  <span>mutate</span><span>(</span>
    <span>avg.10.years</span> <span>=</span> <span>rollmeanr</span><span>(</span><span>x</span> <span>=</span> <span>temperature</span><span>,</span>
                             <span>k</span> <span>=</span> <span>10</span><span>,</span>
                             <span>fill</span> <span>=</span> <span>NA</span><span>),</span>
    <span>check</span> <span>=</span> <span>temperature</span> <span>>=</span> <span>avg.10.years</span>
  <span>)</span>

<span>nrw</span> <span>%>%</span> <span>ggplot</span><span>(</span><span>aes</span><span>(</span><span>x</span><span>=</span><span>year</span><span>,</span> <span>y</span><span>=</span><span>temperature</span><span>))</span><span>+</span>
  <span>geom_line</span><span>(</span><span>aes</span><span>(</span><span>group</span> <span>=</span> <span>region</span><span>,</span> <span>color</span> <span>=</span> <span>temperature</span> <span>></span> <span>avg.10.years</span><span>))</span><span>+</span>
  <span>geom_line</span><span>(</span><span>aes</span><span>(</span><span>x</span><span>=</span><span>year</span><span>,</span> <span>y</span><span>=</span><span>avg.10.years</span><span>))</span></code></pre>
</figure>
<p><img alt="DWD temperature data, messy" src="http://i1.wp.com/opiateforthemass.es/images/safferli/temp-mess.png?w=456" data-recalc-dims="1"></p>
<p>The usual Brownian motion around an average. The average which has been increasing since around the 1970s, matching the world-wise <a target="_blank" href="http://climate.nasa.gov/vital-signs/global-temperature/">NASA data</a>.</p>
<p>First problem, and the most painful one: I was wrong. March 2016 was in fact slightly cooler than the 10-year average.</p>
<p>Second problem: grouping the data by a cutoff point (the 10-year average temperature) works, but since it’s point-wise data, ggplot colours the line continuously from one dot to the next, even when it crosses the average temperature line and should change colour exactly at that point.</p>
<p>Is there a way for ggplot to handle this? With some quick googleing, I found the answer on <a target="_blank" href="http://stackoverflow.com/questions/9543402/how-to-obtain-multiple-colours-for-geom-line-conditional-on-a-specific-value">stackoverflow</a>: not directly, but you can generate a new “high resolution” dataset which would allow ggplot to set the cutoff point exactly at the intersection between the two lines.</p>
<p>I can claim no ownership on the solution, the credit goes fully to the stackoverflow user <a target="_blank" href="http://stackoverflow.com/users/489704/jbaums">jbaums</a>. I took his solution and “updated” it to newer R (<code>dplyr</code> and pipes, specifically; the answer was written four years ago, at which time my solution setup wasn’t available).</p>
<figure>
<pre><code data-lang="r"><span>#new data for testing
</span><span>tt</span> <span>&lt;-</span> <span>nrw</span>
<span># resolution scale (0.01 = factor 100)
</span><span>yres</span> <span>&lt;-</span> <span>0.01</span>

<span>tt</span> <span>%&lt;>%</span> 
  <span>mutate</span><span>(</span>
    <span># rename to use in function later
</span>    <span>x</span> <span>=</span> <span>year</span><span>,</span>
    <span>y</span> <span>=</span> <span>temperature</span><span>,</span>
    <span>z</span> <span>=</span> <span>avg.10.years</span><span>,</span> 
    <span># generate point coordinates for line end
</span>    <span>x2</span> <span>=</span> <span>lead</span><span>(</span><span>x</span><span>),</span>
    <span>y2</span> <span>=</span> <span>lead</span><span>(</span><span>y</span><span>),</span>
    <span>z2</span> <span>=</span> <span>lead</span><span>(</span><span>z</span><span>)</span>
  <span>)</span> <span>%>%</span> 
  <span># remove last row from group
</span>  <span>filter</span><span>(</span><span>!</span><span>row_number</span><span>()</span> <span>==</span> <span>n</span><span>())</span> <span>%>%</span> 
  <span># remove years in which z is NA
</span>  <span># (chokes the summarising later)
</span>  <span>filter</span><span>(</span><span>!</span><span>x</span> <span>%</span><span>in</span><span>%</span> <span>.</span><span>$</span><span>x</span><span>[</span><span>is.na</span><span>(</span><span>z</span><span>)])</span> <span>%>%</span> 
  <span># get high-resolution versions of the inter-point data
</span>  <span>rowwise</span><span>()</span> <span>%>%</span> 
  <span>do</span><span>(</span>
    <span># for each row, generate dataframe with the low-res data,
</span>    <span># which gets automatically filled up to the high-res data rows
</span>    <span>data.frame</span><span>(</span>
      <span>year</span> <span>=</span> <span>.</span><span>$</span><span>year</span><span>,</span>
      <span>month</span> <span>=</span> <span>.</span><span>$</span><span>month</span><span>,</span>
      <span>region</span> <span>=</span> <span>.</span><span>$</span><span>region</span><span>,</span>
      <span>temperature</span> <span>=</span> <span>.</span><span>$</span><span>temperature</span><span>,</span>
      <span>x</span> <span>=</span> <span>.</span><span>$</span><span>x</span><span>,</span>
      <span>x2</span> <span>=</span> <span>.</span><span>$</span><span>x2</span><span>,</span>
      <span># make a sequence (linear line) of the interpoint data
</span>      <span># we will use this data length as &quot;basis&quot; for x.hr and z.hr, below
</span>      <span>y.hr</span> <span>=</span> <span>seq</span><span>(</span><span>.</span><span>$</span><span>y</span><span>,</span> <span>.</span><span>$</span><span>y2</span><span>,</span> <span>yres</span><span>*</span><span>sign</span><span>(</span><span>.</span><span>$</span><span>y2</span><span>-</span><span>.</span><span>$</span><span>y</span><span>)),</span>
      <span>avg.10.years</span> <span>=</span> <span>.</span><span>$</span><span>z</span><span>,</span>
      <span>z2</span> <span>=</span> <span>.</span><span>$</span><span>z2</span>
    <span>)</span> <span>%>%</span> 
      <span>mutate</span><span>(</span>
        <span># make a sequence (linear line) of the remaining interpoint data
</span>        <span># take the rownumber of data needed from y.hr, above
</span>        <span>x.hr</span> <span>=</span> <span>seq</span><span>(</span><span>min</span><span>(</span><span>x</span><span>),</span> <span>min</span><span>(</span><span>x2</span><span>),</span> <span>length.out</span> <span>=</span> <span>nrow</span><span>(</span><span>.</span><span>)),</span>
        <span>z.hr</span> <span>=</span> <span>seq</span><span>(</span><span>min</span><span>(</span><span>z</span><span>),</span> <span>min</span><span>(</span><span>z2</span><span>),</span> <span>length.out</span> <span>=</span> <span>nrow</span><span>(</span><span>.</span><span>))</span>
      <span>)</span>
  <span>)</span></code></pre>
</figure>
<p>Let’s go through the code step by step (which is made easy since it’s piped code).</p>
<ol>
<li>start a new dataset for testing, and set the high-resolution scale factor (0.01, or factor 100 in this case)</li>
<li>generate point coordinates for the line end point by using <code>lead()</code></li>
<li>remove last row from each group</li>
<li>remove any <code>x</code> values for which <code>z</code> have <code>NA</code>s. These result from the rolling average, and will choke the <code>seq()</code> calls later</li>
<li>for each row, we now have the start point and end point of each “section”: <code>(x1, x2)</code>, <code>(x2, x3)</code>, …
<ul>
<li>we now want to generate a “high resolution” version of these point, essentially filling 100 data points linearly between <code>(x1, x2)</code>, and <code>(x2, x3)</code>, and so on</li>
<li>to do this, we use <code>rowwise() %>% do()</code></li>
<li>each row pulls the start points and the end points of each data section, and then using these as start and end points for a sequence of length <code>yres</code>, the high-resolution scale factor</li>
<li>all of this is done as a <code>data.frame()</code>, and by the magic of <code>rowwise-do()</code> the result is a single, large dataframe with everything inside!</li>
</ul>
</li>
</ol>
<p>With the high-res data now set up, we can graph the chart in a much nicer way:</p>
<figure>
<pre><code data-lang="r"><span>ggplot</span><span>()</span><span>+</span>
  <span># high-res version of the temperature data
</span>  <span>geom_line</span><span>(</span><span>data</span> <span>=</span> <span>tt</span><span>,</span> <span>aes</span><span>(</span><span>x</span><span>=</span><span>x.hr</span><span>,</span> <span>y</span><span>=</span><span>y.hr</span><span>,</span> <span>group</span> <span>=</span> <span>region</span><span>,</span> <span>color</span> <span>=</span> <span>y.hr</span> <span>></span> <span>z.hr</span><span>))</span><span>+</span>
  <span># add rolling average as line
</span>  <span>geom_line</span><span>(</span><span>data</span> <span>=</span> <span>tt</span><span>,</span> <span>aes</span><span>(</span><span>x</span><span>=</span><span>year</span><span>,</span> <span>y</span><span>=</span><span>avg.10.years</span><span>,</span> <span>color</span> <span>=</span> <span>&quot;10y rolling avg&quot;</span><span>))</span><span>+</span>
  <span>guides</span><span>(</span>
    <span>colour</span> <span>=</span> <span>guide_legend</span><span>(</span><span>&quot;temperature is larger than n10-year rolling avg&quot;</span><span>)</span>
  <span>)</span><span>+</span>
  <span># change line colours -- black for the line, red if above, green if below
</span>  <span>scale_color_manual</span><span>(</span><span>values</span><span>=</span><span>c</span><span>(</span><span>&quot;#000000&quot;</span><span>,</span> <span>&quot;#A5CF35&quot;</span><span>,</span> <span>&quot;#C44B4B&quot;</span><span>))</span><span>+</span>
  <span>labs</span><span>(</span><span>title</span> <span>=</span> <span>&quot;March temperatures in North Rhine-Westphalia&quot;</span><span>,</span>
       <span>x</span> <span>=</span> <span>&quot;&quot;</span><span>,</span> <span>y</span> <span>=</span> <span>&quot;temperature (March monthly average)&quot;</span><span>)</span>
<span># save everything, using golden ratio
</span><span>ggsave</span><span>(</span><span>file</span><span>=</span><span>&quot;temperature.png&quot;</span><span>,</span> <span>width</span> <span>=</span> <span>30</span><span>,</span> <span>height</span> <span>=</span> <span>30</span><span>/</span><span>((</span><span>1</span><span>+</span><span>sqrt</span><span>(</span><span>5</span><span>))</span><span>/</span><span>2</span><span>),</span> <span>units</span> <span>=</span> <span>&quot;cm&quot;</span><span>)</span></code></pre>
</figure>
<p><img alt="DWD temperature data, nice colour shifts" src="http://i0.wp.com/opiateforthemass.es/images/safferli/temperature.png?w=456" data-recalc-dims="1"></p>
<p>Ooh, it’s so pretty! ☺</p>
<h2>Conclusion</h2>
<p>One can clearly see the underlying trend of temperature increases. I believe I’ll have to find a job in <a target="_blank" href="https://en.wikipedia.org/wiki/Hammerfest">Hammerfest</a> in the coming years…</p>
<p>What could have been answered with a calculator and 2 minutes on the website turned out to be two day project to pull and clean data, fix and report bugs, and learn even more about <code>ggplot()</code>. Time well spent, I’d say!</p>
<p>Code and data for this analysis is availabe on <a target="_blank" href="https://github.com/safferli/dwd">github</a>, as always.</p>
<p><a target="_blank" href="http://opiateforthemass.es/articles/dwd-daten/">German Temperature Data</a> was originally published by Kirill Pomogajko at <a target="_blank" href="http://opiateforthemass.es/articles/dwd-daten/">Opiate for the masses</a> on May 12, 2016.</p>
<p><div><div><div><div data-width="55" data-show-faces="false" data-send="true" data-layout="box_count" data-height="62" data-href="http://www.r-bloggers.com/german-temperature-data/"></div></div><div><a data-url="http://www.r-bloggers.com/german-temperature-data/" data-counturl="http://www.r-bloggers.com/german-temperature-data/" data-count="vertical" data-text="German Temperature Data" class="twitter-share-button" data-via="rbloggers" href="https://twitter.com/share"></a></div><div></div></div><div></div></div>
<div>
<div>To <strong>leave a comment</strong> for the author, please follow the link and comment on their blog: <strong><a href="http://opiateforthemass.es/articles/dwd-daten/"> Opiate for the masses</a></strong>.</div>
<hr>
<a rel="nofollow" href="http://www.r-bloggers.com/">R-bloggers.com</a> offers <strong><a rel="nofollow" href="http://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="http://www.r-project.org/">R</a> news and <a title="R tutorials" rel="nofollow" href="http://www.r-bloggers.com/?s=tutorial">tutorials</a> on topics such as: <a title="Data science" rel="nofollow" href="http://www.r-bloggers.com/?s=data%20science">Data science</a>, <a title="Big Data" rel="nofollow" href="http://www.r-bloggers.com/?s=Big%20Data">Big Data, <a title="R jobs" rel="nofollow" href="http://www.r-users.com/">R jobs</a>, visualization (<a title="ggplot and ggplot2 tutorials" rel="nofollow" href="http://www.r-bloggers.com/?s=ggplot2">ggplot2</a>, <a title="Boxplots using lattice and ggplot2 tutorials" rel="nofollow" href="http://www.r-bloggers.com/?s=boxplot">Boxplots</a>, <a title="Maps and gis" rel="nofollow" href="http://www.r-bloggers.com/?s=map">maps</a>, <a title="Animation in R" rel="nofollow" href="http://www.r-bloggers.com/?s=animation">animation</a>), programming (<a title="RStudio IDE for R" rel="nofollow" href="http://www.r-bloggers.com/?s=RStudio">RStudio</a>, <a title="Sweave and literate programming" rel="nofollow" href="http://www.r-bloggers.com/?s=sweave">Sweave</a>, <a title="LaTeX in R" rel="nofollow" href="http://www.r-bloggers.com/?s=LaTeX">LaTeX</a>, <a title="SQL and databases" rel="nofollow" href="http://www.r-bloggers.com/?s=SQL">SQL</a>, <a title="Eclipse IDE for R" rel="nofollow" href="http://www.r-bloggers.com/?s=eclipse">Eclipse</a>, <a title="git and github, Version Control System" rel="nofollow" href="http://www.r-bloggers.com/?s=git">git</a>, <a title="Large data in R using Hadoop" rel="nofollow" href="http://www.r-bloggers.com/?s=hadoop">hadoop</a>, <a title="Web Scraping of google, facebook, yahoo, twitter and more using R" rel="nofollow" href="http://www.r-bloggers.com/?s=Web+Scraping">Web Scraping</a>) statistics (<a title="Regressions and ANOVA analysis tutorials" rel="nofollow" href="http://www.r-bloggers.com/?s=regression">regression</a>, <a title="principal component analysis tutorial" rel="nofollow" href="http://www.r-bloggers.com/?s=PCA">PCA</a>, <a title="Time series" rel="nofollow" href="http://www.r-bloggers.com/?s=time+series">time series</a>, <a title="finance trading" rel="nofollow" href="http://www.r-bloggers.com/?s=trading">trading</a>) and more...
<img height="1" alt="" width="1" src="http://feeds.feedburner.com/~r/RBloggers/~4/GziniWeXA4Q"></a></div>