Message Boards Message Boards

[GIF] Animating climate change through an annual window

GROUPS:

Yesterday I posted a little piece of code to visualise the changing mean temperature of the earth. The animation tried to point out "how the temperature spirals out of control". In order to account for positive and negative deviations from the "standard" temperature, we needed to use a transformation from temperature difference to radius. This was non-linear and one might argue that lower temperatures are underrepresented.

Here is an alternative attempt.

enter image description here

I do not post it in the same thread as the last post, because the gif is quite large and adding more the original thread might make it very slow to load. In spite of that I only have uploaded a lower quality video with every second frame. You can download the full quality file from here. Here's the implementation of the new idea, which I have seen in a similar fashion somewhere online:

tempdata = Transpose[{Interpreter["Date"][#[[All, 1]]], #[[All, 2]]}] &@ 
Import["http://www.metoffice.gov.uk/hadobs/hadcrut4/data/current/time_series/HadCRUT.4.4.0.0.monthly_ns_avg.txt", "Data"];
frames2 = 
Table[ListLinePlot[#, PlotStyle -> Table[Opacity[0.7/(1 + Min[k, Length[Partition[tempdata[[All, 2]], UpTo[12]]]] - j)^2 + 0.3], {j, 1, Min[k, Length[Partition[tempdata[[All, 2]], UpTo[12]]]]}], PlotRange -> {{1, 12}, Evaluate[MinMax[tempdata[[All, 2]]] + {-0.1, 0.1}]}, 
Frame -> True, Background -> Black, FrameStyle -> Directive[White, Bold, 16], LabelStyle -> Directive[White, Bold, 16], ImageSize -> Large, 
ColorFunction -> (ColorData["Temperature"][2^(#2) - 0.5] &), ColorFunctionScaling -> False, FrameTicks -> {{True, None}, {Transpose@{Range[12], {"Jan", "Feb", "Mar", "Apr", "May", "Jun", "Jul", "Aug", "Sep", "Oct", "Nov", "Dec"}}, None}}, FrameLabel -> {"Month", "\[CapitalDelta]T"}, 
ImagePadding -> 75, AspectRatio -> 1, Epilog -> {Text[Style[1849 + Min[k, Length[Partition[tempdata[[All, 2]], UpTo[12]]]], Red, 21], {2.5, 1.}]}] & @
Partition[tempdata[[All, 2]], UpTo[12]][[1 ;; Min[k, Length[Partition[tempdata[[All, 2]], UpTo[12]]]]]], {k, 1, Length[Partition[tempdata[[All, 2]], UpTo[12]]] + 10}];

You an list animate it:

ListAnimate[frames2]

I exported the frames:

Monitor[Do[Export["~/Desktop/ClimateGraph/frame" <> ToString[1000 + k] <> ".jpg", frames2[[k]], ImageResolution -> 100], {k, 1, Length[frames2], 1}], k]

and use the terminal command:

convert -delay 15 -loop 0 frame*.jpg animatedfull.gif

to create the animation above.

If you prefer the curves a bit smoother you can use:

frames3 = 
Table[ListLinePlot[#, PlotStyle -> Table[Opacity[0.7/(1 + Min[k, Length[Partition[tempdata[[All, 2]], UpTo[12]]]] - j)^2 + 0.3], {j, 1, Min[k, Length[Partition[tempdata[[All, 2]], UpTo[12]]]]}], PlotRange -> {{1, 12}, Evaluate[MinMax[tempdata[[All, 2]]] + {-0.1, 0.1}]}, InterpolationOrder->2,
Frame -> True, Background -> Black, FrameStyle -> Directive[White, Bold, 16], LabelStyle -> Directive[White, Bold, 16], ImageSize -> Large, 
ColorFunction -> (ColorData["Temperature"][2^(#2) - 0.5] &), ColorFunctionScaling -> False, FrameTicks -> {{True, None}, {Transpose@{Range[12], {"Jan", "Feb", "Mar", "Apr", "May", "Jun", "Jul", "Aug", "Sep", "Oct", "Nov", "Dec"}}, None}}, FrameLabel -> {"Month", "\[CapitalDelta]T"}, 
ImagePadding -> 75, AspectRatio -> 1, Epilog -> {Text[Style[1849 + Min[k, Length[Partition[tempdata[[All, 2]], UpTo[12]]]], Red, 21], {2.5, 1.}]}] & @
Partition[tempdata[[All, 2]], UpTo[12]][[1 ;; Min[k, Length[Partition[tempdata[[All, 2]], UpTo[12]]]]]], {k, 1, Length[Partition[tempdata[[All, 2]], UpTo[12]]] + 10}];

And proceed as before. The last frame looks like this:

enter image description here

Cheers,

Marco

POSTED BY: Marco Thiel
Answer
1 year ago

The simplicity of this is fascinatingly impressive and brutal. Reminds me the films where characters are forced to grasp for a thin layer of air at the ceiling in a tank flooded with water. Nicely done.

POSTED BY: Sam Carrettie
Answer
1 year ago

Like in The Drowning Pool (the scene was also in the book, I believe). I wish few more of Ross MacDonald's Lew Archer books were made into movies.

POSTED BY: Daniel Lichtblau
Answer
1 year ago

enter image description here - another post of yours has been selected for the Staff Picks group, congratulations !

We are happy to see you at the tops of the "Featured Contributor" board. Thank you for your wonderful contributions, and please keep them coming!

POSTED BY: Moderation Team
Answer
1 year ago

On my system, the kernel quits during the import and conversion. I broke the code up into smaller segments and used a local copy and got it to work. I think the problem is that the kernel times out using Interpreter[] so much. If I do the dates first, using Interpreter[], then the temperatures, it works.

I tried this with Mathematica 10.4.1 on a 13 inch and 15 inch MacBook Pro, running tOS X 10.11.5

I agree that the spiral interpretation, although neat looking, distorts the data.

POSTED BY: George Woodrow III
Answer
1 year ago

Dear George,

you are quite right, sometimes there are issues with Interpreter when applied to many entries. In fact I have another version where I use a bit more of input modification and then the function DateObject instead. This is much faster and never appears to cause problems.

Interpreter is more concise though and takes away lot of the fiddling with the data. I don't think that this is an issue of the OS or laptop performance. I used simple MacBook to generate the plot, but also checked on a couple of MacBook Pros, an iMac, a windows machine and a MacPro. The internet speed seems to be important. In another post people suggested that also the Wolfram Server you connect to is quite important for the performance.

Best wishes,

Marco

POSTED BY: Marco Thiel
Answer
1 year ago

Great data visualization, it would be nice to see this on a billboard in some major city, rather than a Cola sign or Ikea furniture advert...

Vince

POSTED BY: Vincent Meade
Answer
1 year ago

Thanks for this great data visualization. I just joined Wolfram pro in order to try to understand the temperature spike starting in October, 2015. Below is the quarterly data from NASA (http://data.giss.nasa.gov/gistemp/tabledata_v3/GLB.Ts+dSST.txt) from 1977-2016 I plotted these with Excel and came up with the following graph.

I would appreciate any thoughts/analysis of this. I am still learning Wolfram and don't know how to enter a data file this large for Wolfram statistical analysis.

Quarterly NASA temperature anomalies--1977-2016

15 27 23 9 8 14 -5 8 6 12 10 27 40 32 24 23 39 35 33 18 21 7 9 14 45 36 22 28 21 24 12 13 3 15 11 12 28 27 13 9 32 23 36 32 49 46 40 30 28 29 27 29 39 28 37 39 44 42 49 37 40 34 16 5 32 30 21 13 18 31 31 40 55 40 47 42 35 33 37 29 36 42 44 62 69 65 72 47 57 33 38 42 44 52 42 35 40 56 54 59 68 71 56 59 56 58 56 65 67 56 38 63 60 67 65 77 65 53 63 69 81 71 60 61 37 59 51 66 56 59 69 71 72 85 64 71 50 61 69 59 49 67 61 76 58 60 63 75 63 80 68 81 82 80 76 96 119

Note that the quarterly data for Apr, May, and June will come out in another week or two. So far the average of April and May is 120.

POSTED BY: Jack Fleck
Answer
1 year ago

Your graph includes sea surface temperatures. US SST's are suspect because they have recently been corrected to account for buoy vs. ship measurements - that's the excuse. Of course the official US position is that global warming is the biggest threat.to our safety.

POSTED BY: Douglas Kubler
Answer
1 year ago

Please consider this graph of global temperatures from 1977 - 2016. Does it seem more likely to you that the trend is exponential or linear or something else? Data source: http://data.giss.nasa.gov/gistemp/tabledata_v3/GLB.Ts+dSST.txt global temperature increases 1977-2016

NASA temperature anomalies--1977-2016 (2016 projected using June, 2016 for remainder of year) 18 7 17 28 33 13 30 15 12 19 34 40 29 44 42 23 24 32 46 34 48 63 42 42 54 63 61 54 69 63 66 54 64 72 60 63 65 74 87 94

POSTED BY: Jack Fleck
Answer
1 year ago

Interesting but the data is a fantasy. HADCRUT4 is a blend of CRUTEM4 land-surface air temperature dataset and the HadSST3 sea-surface temperature (SST) dataset. I can't believe the entire world, including places man had never explored or sailed, had accurate thermometers recording every year. The graph should be labeled "for entertainment only". The 2016 bump is an El Nino phenomena that doesn't upset the political climate so it will not be corrected and homogenized and gridded out of existence.

POSTED BY: Douglas Kubler
Answer
1 year ago

Dear Jack,

thank you very much for posting another data set with the global temperature anomalies.

Please consider this graph of global temperatures from 1977 - 2016. Does it seem more likely to you that the trend is exponential or linear or something else? Data source: http://data.giss.nasa.gov/gistemp/tabledata_v3/GLB.Ts+dSST.txt

Well, that is a very complicated question. Probably "something else" is the best answer here. There is a huge scientific machinery behind the numbers/measurements and the predictions resulting from them. In fact, there is a beautiful documentary on BBC 4 on Climate Change by Numbers where Dr Hannah Fry, Prof Norman Fenton and Prof David Spiegelhalter explain some of the numbers and methods behind the predictions. All three presenters are absolutely brilliant and it's certainly worth watching. There is of course a lot of scientific literature on all of that, but the presenters manage to at least hint at the complexities involved in all of this.

Something similar applies to your question. There is a range of models of climate change. Some are very sophisticated and certainly too complicated for a post on the community. They often have "mechanistically components" that condense knowledge about processes into equations and models. Not all models are equally suited of course and the following is certainly much to simplistic. Mathematica offered stochastic models that we can use to make predictions. These models allow us to check different assumptions we might want to make. Let's use features of Mathematica out of the box. Here's your data:

data = Import["http://data.giss.nasa.gov/gistemp/tabledata_v3/GLB.Ts+dSST.txt", "TSV"];
dataclean = (Select[StringSplit[#, WhitespaceCharacter ..] & /@ 
       Flatten[data[[8 ;;]], 1], Length[#] > 1 && #[[1]] != "Year" &][[1 ;; -5]] //. {a___, c_, b___} /; StringContainsQ[c, "*"] -> {a, "NA", b});

We can plot it - just as you did:

ListLinePlot[
 Transpose@{ToExpression /@ dataclean[[1 ;;, 1]], (1/100. Mean /@ (DeleteCases[#, 
 NA] & /@ (ToExpression /@ (dataclean[[1 ;;, 2 ;; 13]]))))}, Mesh -> All, MeshStyle -> {Red, Thick}, ImageSize -> Large , 
 PlotTheme -> "Marketing", FrameLabel -> {{"Delta Temparature", None}, {"Year", None}}, 
 LabelStyle -> {FontFamily -> "Times New Roman", 18, GrayLevel[0], Bold}] 

enter image description here

We can then make an "averaged time series (over years)".

averagedTimeSeries = Transpose@{ToExpression /@ dataclean[[1 ;;, 1]], (1/100. Mean /@ (DeleteCases[#, NA] & /@ (ToExpression /@ (dataclean[[1 ;;, 2 ;; 13]]))))} tsm = TimeSeriesModelFit[TimeSeries[averagedTimeSeries]]

The last line makes Mathematica look for a good model in some underlying range of stochastic models. These are not perfect for this type of prediction, but might illustrate how important a good model is for the predictions. These lines use the model to forecast. They also compute -based on the model- estimates of the errors of our prediction.

forecast = TimeSeriesForecast[tsm, {60}]
\[Alpha] = .95;
q = Quantile[NormalDistribution[], 1 - (1 - \[Alpha])/2];
err = forecast["MeanSquaredErrors"];
bands = TimeSeriesThread[{{1, -q}.#, {1, q}.#} &, {forecast, 
   Sqrt[err]}]

This can now be plotted:

ListLinePlot[{averagedTimeSeries, forecast, bands["PathComponent", 1],
   bands["PathComponent", 2]}, Filling -> {3 -> {4}}, 
 AspectRatio -> 1/3, PlotStyle -> {Automatic, Automatic, Red, Red}, 
 Frame -> True, Mesh -> All, MeshStyle -> {Red, Thick}, 
 ImageSize -> Full , PlotTheme -> "Marketing", 
 FrameLabel -> {{"Delta Temparature", None}, {"Year", None}}, 
 LabelStyle -> {FontFamily -> "Times New Roman", 18, GrayLevel[0], 
   Bold}]

enter image description here

The green line shows the "prediction" and the red ones are the error bands. The thing is that the model restricts the possible types of behaviour to some extent. Take for example the following "models":

modelexp = NonlinearModelFit[TimeSeries[averagedTimeSeries], (a + c^(l t )), {a, c, l}, t]
modelconst = NonlinearModelFit[TimeSeries[averagedTimeSeries], a, {a}, t]
modellin = NonlinearModelFit[TimeSeries[averagedTimeSeries], a + b t, {a, b}, t]
modelsquare = NonlinearModelFit[TimeSeries[averagedTimeSeries], a + b t + c t^2, {a, b, c}, t]
modelcubic = NonlinearModelFit[TimeSeries[averagedTimeSeries], a + b t + c t^2 + d t^3, {a, b, c, d}, t]
modelsin = NonlinearModelFit[TimeSeries[averagedTimeSeries], amp Sin[o t] + (a + b t + c t^2), {a, b, c, d, o, amp}, t]

Here are the plots:

Grid[{{Show[ListPlot[averagedTimeSeries, PlotTheme -> "Scientific"], 
    Plot[modelconst[t], {t, 1889, 2016}, PlotStyle -> Green, 
     PlotTheme -> "Scientific"], PlotLabel -> "Constant"], 
   Show[ListPlot[averagedTimeSeries, PlotTheme -> "Scientific"], 
    Plot[modellin[t], {t, 1889, 2016}, PlotStyle -> Green, 
     PlotTheme -> "Scientific"], PlotLabel -> "Linear"]}, {Show[
    ListPlot[averagedTimeSeries, PlotTheme -> "Scientific"], 
    Plot[modelsquare[t], {t, 1889, 2016}, PlotStyle -> Green, 
     PlotTheme -> "Scientific"], PlotLabel -> "Quadratic"], 
   Show[ListPlot[averagedTimeSeries, PlotTheme -> "Scientific"], 
    Plot[modelcubic[t], {t, 1889, 2016}, PlotStyle -> Green, 
     PlotTheme -> "Scientific"], PlotLabel -> "Cubic"]}, {Show[
    ListPlot[averagedTimeSeries, PlotTheme -> "Scientific"], 
    Plot[modelexp[t], {t, 1889, 2016}, PlotStyle -> Green, 
     PlotTheme -> "Scientific"], PlotLabel -> "Exponential"], 
   Show[ListPlot[averagedTimeSeries, PlotTheme -> "Scientific"], 
    Plot[modelsin[t], {t, 1889, 2016}, PlotStyle -> Green, 
     PlotTheme -> "Scientific"], PlotLabel -> "Sine-Quadratic"]}}]

enter image description here

None of these models appear to be particularly meaningful and their predictions will be doubtful at best. I think that the data is only meaningful and interpretable if we use suitable model assumptions and the models that are meaningful here are rather complicated; thousands of scientists work on them for example in the IPCC. So there is no easy answer - it is hard work. Of course you might want to argue that whatever the functional form of the trend is, under relatively weak assumptions, you can probably use a low order polynomial to approximate the trend at least for short prediction horizons.

I would compare it to a situation were we have measurements of a falling stone. From basic physics we know that the model that dominates that situation is $$m\ddot x=-m g.$$, with suitable initial conditions. We can solve this:

sol= DSolve[{x''[t] == -g, x[0] == x0, x'[0] == v0}, x[t], t]

which gives:

x[t] -> 1/2 (-g t^2 + 2 t v0 + 2 x0)

where g is the earths acceleration, v0 the initial speed of the stone, and x0 the initial hight. Let's suppose that we drop the stone form a height of 10 metres and with zero velocity:

sol /. {x0 -> 10, v0 -> 0}

which gives:

x[t] -> 1/2 (20 - g t^2)

Ok. So we have a quadratic function with the parameter g:

Plot[Table[x[t] /. sol[[1]] /. {x0 -> 10, v0 -> 0}, {g, 5, 12, 1}], {t, 0, 1}]

enter image description here

So if we had data,

measureddata = 
  Table[{t, x[t] + RandomVariate[NormalDistribution[0, 1]] /. sol[[1]] /. {x0 -> 10, v0 -> 0, g -> 9.81}}, {t, 0, 1.5, 0.002}];

we could obviously fit our model to that to obtain:

FindFit[measureddata, 1/2 (20 - g t^2), g, t]
(*{g -> 9.78587}*)

Here is a plot:

Show[ListPlot[measureddata], Plot[Evaluate@Normal@NonlinearModelFit[measureddata, 1/2 (20 - g t^2), g, t], {t,0, 1.5}, PlotStyle -> Red]]

enter image description here

So we obtain an interpretable result within our model. If our model was inappropriate, we would still get a fit, but would not really be able to interpret it. If we had no idea of the model and would try to fit something very different:

Show[ListPlot[measureddata], Plot[Evaluate@Normal[NonlinearModelFit[measureddata, k (1 - Sin[g t + o]), {k, g, o}, t]], {t, 0, 1.5}, PlotStyle -> Red]]

the results might look good but might not be interpretable.

Show[ListPlot[measureddata], 
 Plot[Evaluate@Normal[NonlinearModelFit[measureddata, k (1 - Sin[g t + o]), {k, g, o}, t]], {t, 0, 1.5}, PlotStyle -> Red]]

enter image description here

Note that the two models perform similarly:

mod1 = NonlinearModelFit[measureddata, k (1 - Sin[g t + o]), {k, g, o}, t];
mod2 = NonlinearModelFit[measureddata, 1/2 (20 - g t^2), g, t];
mod1["RSquared"]
mod2["RSquared"]
(*{0.975745,0.978336}*)

So goodness of fit is not the only point and the predictions might be (very) different for different model assumptions. In the IPCC reports they study different models and different emission scenarios come up with their predictions. When you ask whether the trend is exponential or otherwise, that is basically asking for a suitable model.

Mathematica has lots of tools to estimate model parameters if you have data. It can suggest plausible (often stochastic) models - but you will need some serious modelling and the answer to your question will be quite complicated, I guess.

Cheers,

M.

POSTED BY: Marco Thiel
Answer
1 year ago

Group Abstract Group Abstract