Take the 2-minute tour ×
Mathematica Stack Exchange is a question and answer site for users of Mathematica. It's 100% free, no registration required.

At work, we were discussing when is it the best time to change to winter tires for bikes and/or cars.

Using WeatherData[] and DateListPlot[], it was fairly straightforward for me to create the diagram below:

Mean temperature per day in Stockholm

Fig 1 Mean temperature per day in Stockholm, red is negative and a risk without winter tires.

The code for this is

cityTemp = WeatherData["Stockholm", "MeanTemperature", {{1977, 1, 1}, {2011, 12, 31}, "Day"}];
iceRiskDays = Select[cityTemp, Last[#] < 0 &];

yearStrip[ dataItem_] := {{0, Part[First[dataItem], 2], Last[First[dataItem]]}, Last[dataItem]}

DateListPlot[{yearStrip[#] & /@ cityTemp, yearStrip[#] & /@ iceRiskDays} ]

My question is: How do I calculate for each day the proportion of values for that day that are below 0 Celsius? (e.g. for a date at the end of November, the proportion is likely to be bigger than 0.5)?

My attempts to do this ended with trying to create separate lists for each date, but I felt that this was a less elegant way and also creates less "fit" with DateListPlot.

share|improve this question
    
You know about Count[]? –  J. M. Oct 14 '12 at 12:29
    
Yes, but where I get stuck is how to "map" Count to the list and to "count" the number of values for a date and also count the number of values that - in this case - are negative. –  FredrikD Oct 14 '12 at 12:33
3  
Something like Count[cityTemp, {d_List, t_?Negative}]? BTW, you might also be interested in the functions Split[]/SplitBy[]. –  J. M. Oct 14 '12 at 12:39
1  
+1, But your question contains a wrong assumption: As soon as the temperature drops below +7 C the rubber of summer tires is too hard and winter tires are recommended ;-) –  halirutan Oct 14 '12 at 15:17
2  
Another misleading assumption is that the only days with risk of ice are those with mean temperatures below $0$. A better indicator would be days with a low near $0$. (Ice can be common even with lows slightly above $0$, but a low below $0$ is a pretty clear indicator that moisture will freeze. :-) –  whuber Oct 14 '12 at 15:26

4 Answers 4

up vote 7 down vote accepted

Nice question. First I wrote your code in this way. Where we don't need the function yearStrip

{date,year,month,day,temp}={1,1,2,3,2}
cityTemp=WeatherData["Stockholm","MeanTemperature",{{1977,1,1},{2011,12,31},"Day"}];
cityTemp[[All,date,year]]=0
iceRiskDays=Select[cityTemp,#[[temp]]<0&];
DateListPlot[{cityTemp,iceRiskDays}]

Where we get the same result:

enter image description here

Now we can group the day-month information using GatherBy in this way:

dataGather=GatherBy[cityTemp,#[[date,{month,day}]]&];
getNegativeProportion[data_]:=N@Total[UnitStep[-data]]/Length[data]
negativeProportionList=Sort[{#[[1,1]],getNegativeProportion[#[[All,temp]]]}&/@dataGather];
movAvg=MovingAverage[{N@AbsoluteTime@#[[date]],#[[temp]]}&/@negativeProportionList,14];
DateListPlot[{negativeProportionList,movAvg},Joined-> True,PlotStyle-> {Blue,Red}]

Where we get:

enter image description here

With the red line as the moving average for 14 days!

share|improve this answer
    
I am new to MMA, so I have to ask for a pointer to the documentation of first assignment {date, ...} = {1,1..} –  FredrikD Oct 14 '12 at 16:10
1  
Welcome to Mathematica. This is just a form of make the rest of the code simpler to read when using Parts. It's equivalent to do date=1; year=1 and so on. –  Murta Oct 14 '12 at 18:01

Following @J.M. if you are after the likelihood of icy days as a function of month

ff[x_] := Select[cityTemp, #[[1, 2]] == x &] // (Count[#, {d_List, t_?Negative}]/Length[#]) &

{Range[12] - 1/2, ff /@ Range[12]} // Transpose // 
ListLinePlot[#, InterpolationOrder -> 0,AxesLabel -> {"Month", "likelihood"}] &

which suggests in January its about 15% (from where I live) and 60 % in Stockholm (Brrr...)

Mathematica graphics

Comparing Paris to Stockholm (using the nice package autoLegend)

ff[x_, temp_] := Select[temp, #[[1, 2]] == x &] // (Count[#, {d_List, t_?Negative}]/
 Length[#]) &

{{Range[12] - 1/2, ff[#, cityTempS] & /@ Range[12]} // Transpose,
{Range[12] - 1/2, ff[#, cityTempP] & /@ Range[12]} // Transpose} // 
ListLinePlot[#, InterpolationOrder -> 0,
AxesLabel -> {"Month", "likelihood"}] &//autoLegend[#, {"Stockholm","Paris"}] &

Mathematica graphics

Stealing yet again another suggestion of @J.M. we can look at global warming ;-) This time we look at the faction of days with ice as a function of year:

 dat = Map[Count[#, {d_List, t_?Negative}]/Length[#] &, 
  SplitBy[cityTempP, {#[[1, 1]] &, #[[1, 2]] &}], {2}] // Rest;
 fP = LinearModelFit[Mean /@ dat, {1, x}, x];

 dat = Map[Count[#, {d_List, t_?Negative}]/Length[#] &, 
 SplitBy[cityTempS, {#[[1, 1]] &, #[[1, 2]] &}], {2}] // Rest;
 fS = LinearModelFit[Mean /@ dat, {1, x}, x];

 Show[Plot[{fS[x], fP[x]}, {x, 1, 31}, AxesLabel -> {"year", "Likelihood"}],
 ListPlot[{fS["Data"], fP["Data"]}]]// autoLegend[#, {"Stockholm", "Paris"}] &

Mathematica graphics

And finally for the sake of contributing something original (the code for FilledErrorPlot), let us launch the climate-skeptic versus climato-believer debate while illustrating the famous Lies, damned lies, and statistics methodology on two different ways of plotting the same data; guess which is which ;-)

Show[dat // Transpose // FilledErrorPlot[#] &,
Plot[{fS[x]}, {x, 1, 33}, 
PlotStyle -> Directive[AbsoluteThickness[2], ColorData[10][1]], 
AxesLabel -> {"year", "Likelihood"}], AspectRatio -> 1.5, 
PlotRange -> {0.035, 0.45}, AxesOrigin -> {0., 0.035}, 
AxesLabel -> {"year", "Fraction of \n icy days"}]

Mathematica graphics

versus

Show[dat // Transpose // FilledErrorPlot[#, ErrorOnMean -> False] &,
Plot[{fS[x]}, {x, 1, 33}, 
PlotStyle -> Directive[AbsoluteThickness[2], ColorData[10][1]], 
AxesLabel -> {"year", "Likelihood"}], AspectRatio -> 0.5, PlotRange -> {-0.5, 1}, 
AxesLabel -> {"year", "Fraction of \n icy days"}]

Mathematica graphics

share|improve this answer
2  
To compute the proportions for all the months of those years: Map[Count[#, {d_List, t_?Negative}]/Length[#] &, SplitBy[cityTemp, {#[[1, 1]] &, #[[1, 2]] &}], {2}] –  J. M. Oct 14 '12 at 12:55
    
Yes, it is a likelyhood function that I am after but for a days (there is about 12000 / 400 = 30 values per days which is ok for me). If I understand your counting to months is done in the {Range[12]... code, will experiment with this. (Stockholm is considered to be in the warm, comfortable, south part of Sweden..) –  FredrikD Oct 14 '12 at 13:07
    
will try autoLegend. Check the order of "Paris" and "Stockholm" in the AxesLabel part, I think you put in the wrong order. –  FredrikD Oct 14 '12 at 13:16
    
@FredrikD oops... well spotted. –  chris Oct 14 '12 at 13:22

Here is a slightly different approach using GatherBy.

Reformat dates to be just {month, day}:

shortDates = {DateString[First@#, { "Month","Day",}], Last@#} & /@  cityTemp;

Group the values by day and month, count the days on which the temperature is not positive:

probs = Count[#[[All, 2]], t_?NonPositive]/Length@# & /@    
        GatherBy[Sort@shortDates, #[[1]] &];

Plot the likelihoods:

   ListLinePlot[probs, AxesLabel -> { "Day Number","Likelihood"}]

Mathematica graphics

You can calculate different likelihoods by adjusting the short date conversion string, for a month value just switch the date reformatting to months:

shortDates = {DateString[First@#, {"Month"}], Last@#} & /@ cityTemp;
share|improve this answer

If my understanding of the question is correct, this does the job:

cityTemp = WeatherData["Stockholm", "MeanTemperature",
                       {{1977, 7, 1}, {2011, 12, 31}, "Day"}];

(* proportions; I ignore February 29 for this *)
props = Array[(Count[#, {_List, _?Negative}]/Length[#]) &[Cases[cityTemp,
              {Prepend[Rest[DatePlus[{2011, 1, 1}, #]], yr_], temp_}, Infinity]] &, 365, 0];

Take[N[props], 5] (* sample *)
   {0.705882, 0.588235, 0.705882, 0.676471, 0.705882}

To plot these proportions, you can do this:

DateListPlot[props, {{2011, 1, 1}, {2011, 12, 31}, "Day"},
             FrameLabel -> {None, "Proportion"}, Joined -> True]

plot of proportion of cold days

Pick the day of the year with the most number of "cold" instances:

Pick[Array[Rest[DatePlus[{2011, 1, 1}, #]] &, 365, 0], props, Max[props]]
   {{2, 18}}

So, at least for Stockholm, February 18 is often a cold day.

share|improve this answer

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

Not the answer you're looking for? Browse other questions tagged or ask your own question.