Reproducing the XKCD "Consensus New Year" plot

Question

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!

Answer 1

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]

Answer 2

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[
    Entity["AdministrativeDivision", {EntityProperty[
        "AdministrativeDivision", "ParentRegion"] -> #}], 
    "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[
  Thread[{Keys[results], 
    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]]

Answer 3

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[#[[1]], 15], #[[2]]} & /@ 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.

Answer 4

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]] &, #[[-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` (UTC`Z`)"];

    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"]

 ]

]