# Reproducing the XKCD “Consensus New Year” plot

This xkcd strip for December 31, 2018 had the following plot: How may I use CountryData with DateListPlot (or other related functions) to reproduce this plot?

I know I am able to retrieve the population for each country from a call to CountryData, but I had been stuck trying to figure out how much of a country's population is in a given time zone. Is there any function I am missing to be able to make this plot? I'll appreciate any pointers!

Building off the other answer, CountryData has the data for the time zones for each country, as well as their population. So we can split each country proportionally into its individual timezones:

countries = CountryData[];
populations = CountryData[#, "Population"] & /@ countries;

timeZones = CountryData[#, "TimeZones"] & /@ countries;

data = Flatten[SortBy[Transpose[{timeZones, populations}] /. {a_, b_} :>
Thread[{a, b/Length[a]}] /; ListQ[a], First], 1];

gathered = Reverse@SortBy[
Map[{#[[1, 1]], Total[#[[All, 2]]]} &, GatherBy[data, First]], First];

houroffsets = (Quantity[#, "Hours"] & /@ (-gathered[[All, 1]])) +
Quantity[14 - 5, "Hours"];
ESTtimes = DateString[DatePlus[DateObject[{2018, 12, 31, 10, 0, 0},
TimeZone -> -5], #],
{"Hour12", ":", "Minute", "AMPM", " EST"}] & /@ houroffsets;
xticks = Transpose[{houroffsets, ESTtimes}][[{10, 18, 26, 33, 36, 40}]]

plot = ListLinePlot[Transpose[{houroffsets, Accumulate[gathered[[All, 2]]]/Total[gathered[[All, 2]]]}],
FrameTicks -> {xticks, {{0, "0%"}, {0.2, "20%"}, {0.4, "40%"}, {0.6,
"60%"}, {0.8, "80%"}, {1.0, "100%"}}}, Joined -> True,
PlotLabel -> "PERCENTAGE OF THE WORLD'S POPULATION LIVING IN 2019",
Frame -> {True, True, False, False}, InterpolationOrder -> 1,
ImageSize -> 800, GridLines -> {xticks[[All, 1]], {0.5}}, Axes -> False] Note that this graph differs from the xkcd one, as the spaces on the xkcd plot are not spaced according to the actual time, just based on the individual timezones.

Additionally we can use the xkcdconvert code defined in this answer to change the style (although the labels start to overlap at this size).

xkcdConvert[plot] • Hallo, this looks very good! (Thank you also to David for first response.) I am still not comfortable with just associating a country's population and time zone, however. As an instance, Russia spans a lot of time zones, and the population is different in e.g. St. Petersburg or Siberia. Still, this is appreciated! – 河上小太郎 Jan 3 '19 at 10:22
• Sure, that is fair. You might not want to accept my answer then, and see if anyone else can come up with something better. I guess how much you want to reproduce the xkcd plot exactly and how much you want an actual answer to the same question is something you should clarify. – SPPearce Jan 3 '19 at 10:23

While it takes a lot longer to get the data, one can go a bit more granular by using "AdministrativeDivision" rather than "Country". This could likely be done more efficiently, but I break the steps down quite a bit here:

countries = EntityList["Country"];
divisions =
Map[EntityValue[
"Entities"] &, countries];
flat = divisions // Flatten;


Get the associated populations:

tzp = EntityValue[flat, {"TimeZones", "Population"}];


Now a lot of stuff happens in the following: first I get rid of missing data, then I use the timezone information to convert to a time and date (using TimeZoneConvert), finally I group the results and take the Total. (You might have noticed that I make an approximation here, if an AdministrativeDivision contains multiple time zones, I just take the first one. I don't think this is too terrible to do).

results =
KeySort@GroupBy[MapAt[
TimeZoneConvert[
DateObject[{2019, 1, 1, 0, 0}, TimeZone -> First[#]], -5] &,
DeleteMissing[tzp, 1, 1], {All, 1}], First -> Last, Total];

Short[results] Now to visualize it. I use Block to set the time zone used by DateListPlot to be UTC-5.

Block[{\$TimeZone = -5},
DateListPlot[
Accumulate@Normalize[QuantityMagnitude@Values@results, Total]}],
InterpolationOrder -> 0, Frame -> True, Axes -> False,
, DateTicksFormat -> {"Hour12Short", ":", "Minute", "AMPM", " EST",
"\n", "MonthNameShort", " ", "DayShort", "st"},
FrameTicks -> {{{{0, "0%"}, {1/2, "50%"}, {1, "100%"}},
Automatic}, {{DateObject[{2018, 12, 31, 5}],
DateObject[{2018, 12, 31, 11}],
DateObject[{2018, 12, 31, 13, 30}],
DateObject[{2018, 12, 31, 19, 0}],
DateObject[{2019, 1, 1, 0, 0}], DateObject[{2019, 1, 1, 3, 0}],
DateObject[{2019, 1, 1, 6, 0}]}, None}}, ImageSize -> 800,
FrameTicksStyle -> {Directive["Label", 12],
Directive["Label", 10]}, AspectRatio -> 1/2,
GridLines -> {Automatic, {1/2}},
PlotLabel -> "Percentage of the World's Population Living in 2019",
LabelStyle -> 18]] • Thank you so much, I had not known of AdministrativeDivision before! (This program is surely vast.) Your simplifying assumption should hopefully not introduce too much error in the plot. (I now actually wonder which divisions span the most time zones.) – 河上小太郎 Jan 3 '19 at 18:57

Here is the outline of the solution, illustrated for just the UN countries. (Replace with CountryData[] to get all countries.)

List of UN countries:

theCountries = CountryData["UN"];


The longitudes and populations of these countries:

theLongs = Longitude[#] & /@ theCountries;
thePops = CountryData[#, "Population"] & /@ theCountries;


The data sorted by longitude (going around the world):

allData = SortBy[Transpose[{theLongs, thePops}], First]


Rounding each country to be in a time zone (15 degree separation in longitude):

roundedData = {Round[#[], 15], #[]} & /@ allData


Accumulating (and normalizing) the populations as we go "around the world":

finalData =
Transpose[{roundedData[[All, 1]],
Accumulate[roundedData[[All, 2]]]/
Last[Accumulate[roundedData[[All, 2]]]]}]


Plot it:

        ListPlot[finalData,
Joined -> True,
Ticks -> {Transpose[{Range[-165, 165, 15],
Table[Rotate["GMT+" <> ToString[i], π/2], {i, 23}]}],
{{0.4,  "40%"}, {0.6, "60%"}, {0.8, "80%"}, {1.0, "100%"}}},
PlotLabel ->
"CONSENSUS NEW YEAR:  AS OF 1:30PM EASTERN TIME (6:30PM UTC)\n A MAJORITY OF THE WORLD'S POPULATION WILL BE LIVING IN 2019",
ImageSize -> 500] (Of course, this makes the simplification that all a country's population is in the same time zone... of course not quite correct for many large countries.)

If you want to be really ambitious, plot the data in the XKCD style.

• If you reverse the order it resembles the plot very well, it seems they did the same assumptions. – Kuba Jan 3 '19 at 8:43

OK this is an interesting and funny Q. I figured I should make my answer as modular as possible. By 'modular' I refer to the year the "...consensus" is sought for.

It should be noted that CountryData does not report "Population" figures for recent years. That was my original motivation.

As far as the specifics of the implementation are concerned, I am implementing the accepted answer's approach for countries with multiple time zones, namely I allot each timezone an equal fraction of the total population. No other corrections are made.

## Discrepancies with the original drawing

At this point I should note that this answer displays some discrepancies with the original drawing, that I will address in the following.

There are only four vertical gridlines and three horizontal ones. The later are the same as the ones in the original drawing (drawn at 0, 50 and 100% respectively).

As far as the former ones are concerned, I decided to plot-in addition to the lines corresponding to the starting and ending timestamps-the vertical line at midnight before the new year and the gridline corresponding to the consensus time instant (the time when 50% or more of the world population are living in the new year).

The label of the plot contains a reference to the time zone (eg UTC-5) and the format for the date ticks includes the year (not present in the original drawing).

Lastly, there is no frame around the plot because Framed was not cooperating nice with xkcdConvert and the label below the image was not reproduced because it didn't look good (to me) in the absence of a frame.

## Examples

The package that contains the function plotPopYYYY can be found here (which you can load directly by evaluating Get["https://pastebin.com/raw/a5CYf4ka"]) and also at the bottom of this answer. All the displayed examples were XKCD-ed by the code here.

The procedure to obtain the finished plot is fairly straightforward: evaluate plotPopYYYY[year, zone] for the year you are interested in and the respective time zone; after you obtain the graphic, resize it manually to the dimensions most pleasing to your eyes and then copy that graphic and use it as an input to xkcdConvert. That's all there is to it! Enjoy!

Here are some example output drawings:

1. The last year for which there are available "Population" data (year=2015) and time zone=-5 (the image has been downsized): 2. year = 2009 and zone=2: ## code:

entry[name_, zone_, pop_] := <|
"name" -> name,
"zone" -> zone,
(* population density per zone *)
"density" -> pop
|>

makeEntry[entity_, year_] := Module[{zones = CountryData[entity,"TimeZones"], name = CountryData[entity, "Name"], pop},
(* obvious oversimplification: population uniformly distributed over time zones in each country *)
pop = If[
(* is population missing for requested year? *)
MissingQ[#],
(* then return Missing *)
#,
(* else calc population density per zone *)
#/Quantity[Length[zones], IndependentUnit["zone"]] // N
] &@CountryData[entity, {"Population", year}];

Select[entry[name, #, pop] & /@ zones, FreeQ[#, Missing] &]
]

entriesByZone[year_] := Module[{zf = #["zone"] &},
(* sort (east to west) and gather up entries by zone *)
GatherBy[
SortBy[
Flatten[makeEntry[#, year] & /@ CountryData[]], (-zf[#] &)], zf]
]

(**)

summaryEntry[countries_, count_, zone_, density_] := <|
"countries" -> countries,
"n" -> count,
"zone" -> zone,
"density" -> density
|>

(* input: list of entries by zone *)
makeSummary[entries_List] := Module[{countries, count, density, iCountries, zone},
(* list of countries in the zone *)
countries = Query[All, "name"][entries];
count = Length[countries];
(* population density by zone *)
density = Query[Total, "density"][entries] Quantity[1, IndependentUnit["zone"]];

iCountries = Iconize[countries, "countries"];
(* all entries are from the same zone... *)
zone = Query[1, "zone"][entries];

(* output *)
summaryEntry[iCountries, count, zone, density]
]

(**)

accumulatedEntry[count_, zone_, density_, time_, y_] := <|
"n" -> count,
"zone" -> zone,
"density" -> density,
"time" -> time,
"y" -> y
|>

accumulateByZone[summary_, year_, zone_] := Module[{entry, density, time, y, total},

(* total population *)
total = Query[Total, "density"][summary];

density = Quantity[0, "People"];
time = {year, 12, 31, 23, 0, 0};
y = Quantity[0, "Percent"];

entry = accumulatedEntry[0, zone, density, time, y];

FoldList[
accumulatedEntry[
(* number of countries *)
#1["n"] + #2["n"],
(* zone *)
#2["zone"],
(* population density per zone - in persons *)
#1["density"] + #2["density"],
(* time *)
DatePlus[#1["time"], {#1["zone"] - #2["zone"], "Hours"}],
(* y *)
Quantity[100 (#1["density"] + #2["density"])/total , "Percent"]
] &, entry, summary] // Rest
]

plotPopYYYY[YEAR_: 2009, ZONE_: 2] := Module[{s, offset, a, d, t, y, tfst, tlst, t50, tmdn, xticks, yticks, label, tspec},

s = makeSummary /@ entriesByZone[YEAR - 1];

(* check if there are available data *)
If[
s =!= {},

offset = Query[1, "zone"][s] + 1;

a = accumulateByZone[s, YEAR - 1, offset];
d = Values@Query[All, {"time", "y"}][a];
{t, y} = Transpose[d];

t = DatePlus[#, {ZONE - (offset - 1), "Hours"}] & /@ t;

d = Transpose[{t, y}];

{tfst, tlst} = Through[{#[] &, #[[-1]] &}[t]];

t50 = First[
Query[Select[#[[-1]] >= Quantity[50, "Percent"] &], First][d]];

tmdn = {YEAR - 1, 12, 31, 24, 0, 0};
xticks = {tfst, t50, tmdn, tlst};

yticks = Quantity[#, "Percent"] & /@ {0, 50, 100};

label = StringTemplate[
"PERCENTAGE OF WORLD'S POPULATION LIVING IN Y (UTCZ)"];

tspec = {"Hour12", ":", "Minute", "AMPM", "\n", "MonthNameShort", " ", "Day", "\n", "Year"};

DateListPlot[

d,

PlotStyle -> Thick,
DateTicksFormat -> tspec,
FrameTicksStyle -> Directive[Large],
FrameTicks -> {{yticks, None}, {xticks, None}},
GridLines -> {Automatic, yticks},
FrameLabel -> {{None, None}, {None, label[<|"Y" -> YEAR, "Z" -> ZONE|>]}},
FrameStyle -> Directive[Large],
PlotRange -> Full,
ImageSize -> Large
]
],

Missing["NotAvailable"]

]

]

• This looks good, but why is the code at pastebin instead of in the post? Better to include it in the answer. How far back does the data go? – SPPearce Jan 4 '19 at 10:27
• @KraZug happy new year; I enjoyed your answer; as far as I know the data stop at 2015 (I tested for the range 2001-2015, eg CountryData["France",{"Population",2015}] evaluates to Missing); I posted the code on pastebin for brevity; if there's interest I could include the package here too, I guess – user42582 Jan 4 '19 at 13:20
• Happy new year to you too. I think the general stackexchange policy is for answers to be self-contained if possible, and given that it is only ~150 lines it isn't too much to include. The Population data appears to go back significantly further than 2001 for a number of countries, although it is often only provided at single years as you get further back. Would be very interesting to see this for something like 1900 – SPPearce Jan 4 '19 at 13:28
• You could extract the data that is available for all the years and generate an interpolation function in order to be able to estimate population. And then produce an animation of how it changes over time... – SPPearce Jan 4 '19 at 13:35
• sure, that sounds interesting; I actually thought about doing population projections when I realized there were no current data available but then again I figured that would probably take too long before something worth reporting started taking shape so I focused on producing something modular, for the time being... – user42582 Jan 4 '19 at 21:05