The climate dice are loaded: A replication

^{(Source: Popovich & Pearce (2017): *[It’s Not Your Imagination. Summers Are Getting Hotter](https://www.nytimes.com/interactive/2017/07/28/climate/more-frequent-extreme-summer-heat.html)*, New York Times)}

On July 28, 2017, the New York Times published an article by Nadja Popovich and Adam Pearce with the intriguing title "It’s Not Your Imagination. Summers Are Getting Hotter". As the title suggests, the article makes the claim that summers are getting hotter and hotter. The main argument of the article is the visualization depicted above. It shows "how summer temperatures have shifted toward more extreme heat over the past several decades".

This is a very effective visualization because it communicates relatively complicated scientific information (a distribution of temperatures) in an approachable way. All you need to understand is that temperatures move to the right, which means summers get hotter. (The article also features an animated visualization, which shows how the temperature distribution has changed over the decades.). This post replicates the visualization and reflects on the underlying data preparation and data wrangling activities.

What does the visualization show?¶

In a nutshell, the visualization shows that our chances to experience hot weather during summer months are rising. The visualization that the distribution of temperature measurements on Earth's Northern hemisphere during the summer months (more specifically, June, July, and August) shifts toward hotter temperature over the decades.

The visualization displays a probability density functions for each decade since the 1950s. Probability density functions are fairly complicated statistical instruments, when compared to the typical bar charts or scatter plots in most visualizations. However, the reader does not necessarily need to understand the meaning of probability density functions. Instead, the authors use color (blue, orange, and red) to denote colder and hotter temperatures and animations to visualize the change of temperatures:

• The (grey) area represents the base distribution of temperatures. This baseline is calculated from the summer temperature measurements between 1951 to 1980. The colored areas represent the summer temperature distribution between 2005 and 2015.
• As the probability density function shifts to the right, the colored areas change accordingly. For instance, in the initial probability density function, the reader will probably not notice the small dark red area that denotes extreme temperatures. However, in the final probability density function, the dark red area is substantially larger.
• The animation allows the reader to understand the shift of colors as the visualization cycles through the decades.

I encourage you to visit the original article in the New York Times to experience the full effect of the visualization.

Where does the data come from?¶

The article is a consumer-friendly version of visualizations reported in the following research paper that has been published in the Proceedings of the National Academy of Sciences of the United States of America (PNAS):

James Hansen, Makiko Sato, Reto Ruedy (2012): Perception of Climate Change. Proceedings of the National Academy of Sciences, 109(37) E2415-E2423; DOI: 10.1073/pnas.1205276109

The paper offers the fundamental idea for the visualization. While these original visualizations convey the same content, the visualizations have a distinct "scientific" look and feel ;)

^{(Source: Hansen et al. (2012): *[Perception of Climate Change](http://www.columbia.edu/~mhs119/PerceptionsAndDice/)*, Columbia University)}

It is important to note the following characteristics of the data in Hansen et al. 2012:

• The shift in the average temperature is only for the Northern hemisphere.
• The paper uses averaged temperature data of 250km squares.

A replication¶

I replicate the Hansen et al. (2012) using the following steps:

1. A replication of the data analysis presented by Hansen et al. (2012)
2. A replication of the visualizations created by Hansen et al. (2012)
3. An attempt to get as close as possible to the look and feel of the visualization in the New York Times.

Setup¶

We set up our environment to allow multiple outputs for each code cell. Also, we use the magic command to display matplotlib plots in the Jupyter notebook.

For data wrangling, we use the Pandas library, one of the most versatile workhorses for data analytics. Additionally, we need the NumPy library for general scientific computing.

We format the display of rows and the number of decimals. (These commands setup Jupyter to my personal preferences. They are not required.)

Furthermore, we use altair to render visualizations.

Data¶

The Hansen et al. (2012) paper uses data from the Goddard Institute for Space Studies (GISS), which includes data from the Global Historical Climatology Network (GHCN). Obtaining the GHCN data is surprisingly difficult because it is located on FTP servers, which frequently throw errors (at least on my machine).

However, browsing through the description of the GISS data will lead you to the description of the station data. The station data consists of average monthly temperatures for stations around the globe.

The data for the stations is downloadable here. The data is being continuously updated, so if you download this notebook, some values may be slightly different.

The data is wrapped in a zip file that contains a csv file called step2.v3.csv. The following paper describes the dataset:

Lawrimore, J. H., M. J. Menne, B. E. Gleason, C. N. Williams, D. B. Wuertz, R. S. Vose, and J. Rennie (2011), An overview of the Global Historical Climatology Network monthly mean temperature data set, version 3, J. Geophys. Res., 116, D19121, doi: 10.1029/2011JD016187.

Let's read the data.

A quick look at the data reveals three potential issues:

1. The data has no header information.
2. The data has strange values that initially do not look like temperatures.
3. The data clearly has missing values which are denoted as -9999.

This README file explains the data.

For the fun of it, we nevertheless explore the data and confirm the structure of the data.

Data Exploration¶

Let's interpret, validate, and confirm our understanding of the data.

Column 0¶

The first column seems to be the station identifier. You can confirm this assumption by entering a few station identifiers in the search form here. They check out.

Here are the links to the stations above:

• 425004277140 in Scipio, Utah, USA.
• 425002537350 in Hebron, Nebraska, USA.
• 501947140000 in Cootamundra, New South Wales, Australia.

(If you run this notebook on your machine, these values may change because of the underlying data changes).

Column 1¶

The second column seems to contains years. Let's confirm this.

This looks very much like years ranging from 1880 to 2018.

Columns 2 to 13¶

Then the file has additional 12 columns. They represent the 12 months of the year. Again, you can confirm this using the search form here.

For instance, the average temperature in January in Santa Cruz, California is 9.8 degrees Celcius (49.6 degrees Fahrenheit). The closest station is 425000479160 in Santa Cruz.

The station in Santa Cruz reports an average temperature of 9.83 degrees Celsius in January for the period between 1893 and 2018.

The 12 columns contain the average temperature per month. The values are in hundredths of a degree Celsius. For example, 983 means an average of 9.83 degree Celsius (49.694 degrees Fahrenheit).

-9999 denotes missing values. Again, the README file explains the values (and many more).

^{(Source: J. H. Lawrimore, M. J. Menne, B. E. Gleason, C. N. Williams, D. B. Wuertz, R. S. Vose, and J. Rennie (2011), [An overview of the Global Historical Climatology Network monthly mean temperature data set, version 3](https://agupubs.onlinelibrary.wiley.com/doi/abs/10.1029/2011JD016187), J. Geophys. Res., 116, D19121, doi: [10.1029/2011JD016187](https://doi.org/10.1029/2011JD016187))}

Data Preparation¶

With this interpretation, we rename the original file to reflect our interpretation. We read the data as follows:

• Missing values (na_values) are -9999
• The first column is called station_id
• The second column is called year
• The remaining 12 columns are called jan, feb,mar,apr,may,jun,jul,aug,sep,oct,nov,dec

We have a look at the first five rows of the dataset.

The dataset contains temperature information from 5991 stations.

Data Wrangling¶

The Hansen et al. (2012) paper focusses on the summer months only. The paper argues that the shift of the probability density function is more pronounced in the summer months (See charts above).

We focus our data on the three months of summer: June, July, August.

To make our life easier, we drop rows with any missing values.

Despite removing missing values, we have data for at least one summer for each station.

The probability density functions should represent the full summer. After all, the title of the New York Times article is "Summers Are Getting Hotter". The Hansen et al. (2012) paper is using the mean of the summer months. This is what we do, too.

We create a mean for each row in the summer data frame and store it in a new column called JJA, which stands for

• June,
• July,
• August.

We do not need the columns for the individual months any longer, so we drop them.

The individual stations are all around the globe. A station near the equator will measure vastly different temperatures than a station near the south pole. However, we are only interested in changes in the mean temperature over time.

Thus, we standardize the JJA column by calculating the z-score for each temperature measurement. You calculate z-scores by subtracting the mean from each value and then divide by the standard deviation. The result of this process are values with a mean of 0 and a standard deviation of 1.

A first attempt¶

To calculate the z-score, we need the mean and the standard deviation for each station. We use Pandas' groupby and transform functions to generate the z-score.

The transform function only returns the transformed values but does not change the order, we can simply write it back as a new column.

A quick look at the summer data shows that we now have the z-score for each summer and station.

As expected the mean for each station is 0 and the standard deviation is 1.

We use the hist function to visualize the histogram (a visual representation of the probaility density functions). Since the z-score is a continuous measure, we create one hundred bins.

This looks very much like the end result we would like to have. However, currently, all the temperatures are averages across the whole period. We do not have the baseline yet. The original idea by Hansen et al. (2012) is to use the temperatures of the years between 1951 and 1980 as the baseline.

Understanding the distribution of data across the time periods¶

Before we establish the baseline, we need to have a clear understanding of what subsets of the data are relevant for the baseline. In general, we have three subsets:

• Data prior to the baseline suggested by Hansen et al. (2012)
• Data of the baseline between 1951 and 1980.
• Data after the baseline.

• 3884 stations report data from prior to the baseline.
• 5784 stations report data on the baseline period. This is an increase by 1900 stations compared to prior to the baseline.
• 5151 stations report data beyond the baseline period. This is a decrease by 633 stations compared to the baseline.

We analyze the distribution of stations across the timeperiod.

The following table shows the distribution of stations over the relevant periods (prior to baseline, during baseline, beyond baseline).

Establishing the baseline¶

Currently, the z-scores are calculated using the full dataset. We establish the baseline by calculating the mean and the standard deviation using station data from the time period of 1951 to 1980.

We create a new data set that only contains the baseline data. We copy the summer data, group the data by the station_id, and use the agg function to create the required statistical information mean and standard deviation. Additionally, we create the count to make sure that the mean and standard deviation was derived from a decent number of measurements. We reset the index to obtain a standard data frame.

We have 76 stations that offer no data on the base period. These are the 7 "resurrected" stations and the 69 stations installed after 1980 (See the table above).

Let's have a look at the result.

Now, we add the baseline count, mean, and standard deviation to the original summer data.

The data makes sense. The minimum mean temperature is -66.86 degrees Celsius in Vostok on Antarctica. The maximum mean temperature is 37.28 degrees Celcius in Sibi in Pakistan.

The number of stations also checks out. 5784 of the 5991 stations report data on the baseline. We lose 18 (5784-5766) stations because it is not possible to calculate the standard deviation. These stations have only 1 year within the baseline (either 1951 or 1980).

Calculating the z-score is now very easy. We simply take the JJA, subtract the mean and divide by std.

The min and max values show that there are some outliers.

These outliers appear in the dataset because we have calculated the mean and the standard deviation only on the baseline time period. That means a station that began operation in 1978 only has three observations even though it may have far more observations on the complete dataset. See the following example.

We ignore them for now by restricting the data to +/- 5 standard deviations.

The baseline histograms¶

Let's have a look at the histogram, again.

If we now look at the data since 2000, we notice the shift to the right. This looks promising.

Include the decades¶

We add a column to the data that contains the decade starting with 1951, 1961, 1971 and so on.

Let's have a look at the data.

Let's plot the data per decade.

That is very close to the results in Hansen et al. (2012). It is important to note that we use the complete station data (Not just the data from the Northern hemisphere as in the two articles.)

Focus on the Northern hemisphere¶

The Goddard Institute for Space Studies (GISS) also provides more details on the stations here. This is a text file(!) with very problematic formatting. Simply reading the file does not work. We need a fixed-width reader. Lucky us, Pandas provides one.

Strangely, the station ids from the temperature dataset have an additional zero at the end, which we need to remove.

We have more stations than being used in our temperature dataset. We extract the stations that are in temperature dataset.

We expect latitude in the range between -90 degrees and 90 degrees. Similarly, longitude data should be in the range between -180 degrees and 180 degrees.

Now, we can merge the station information with the weather information.

Let's have a look and drop all rows with missing values.

Next, we need to filter the information so that we only get the Northern hemisphere. The station information explains:

Lat =latitude in degrees, negative = South of Equator

Roughly 83% of all stations are in the Northern hemisphere.

Histograms by decade¶

In the decades leading up to the base period, the temperatures meander around the mean.

The results for the Northern hemisphere are almost identical to our result above. Also, the results are almost identical to the Hansen et al. (2012) paper (first chart on the left):

^{(Source: Hansen et al. (2012): *[Perception of Climate Change](http://www.columbia.edu/~mhs119/PerceptionsAndDice/)*, Columbia University)}

The Southern hemisphere¶

The results for the Southern hemisphere also show the shift to the right. However, the shift is far less pronounced. Also, we have far fewer stations for the Southern hemisphere.

Correcting for Data Biases¶

Hansen et al. (2012) note two characteristics of the results. First, as clearly visible, the temperature moves to the right. Second and more subtle, the distribution of temperatures becomes wider and flatter. Thus, Hansen et al. (2012) argue that the variance of the temperatures increases, too. However, this interpretation is not without challenges.

A response article to Hansen et al. (2012) challenges the second interpretation, the increase in variance, of the data:

Andrew Rhines, Peter Huybers (2013): Frequent summer temperature extremes reflect changes in the mean, not the variance, Proceedings of the National Academy of Sciences, 110(7) E546; DOI: 10.1073/pnas.1218748110

The authors argue that the data is biased towards higher variance for the following reasons:

1. Normalizing time series with the standard deviation of the baseline only introduces a bias because "the mean temperature of individual time series generally differs from zero during the later period". This difference in the means across the time series increases the variance.
2. Trends in the time series also increase the variance.
3. We have fewer stations in the dataset for later periods. Recall that 5784 stations report data on the baseline period but only 5151 report data beyond the baseline.

The original variance ratio between baseline and the time period 1981-2010 is 1.88. Rhines and Huybers (2013) suggest to debias the data in the following way:

1. Removing the sample means independently will reduce the ratio to 1.49.
2. Removing differences in the spatial variance between the baseline and the time period 1981-2010 will reduce the ratio to 1.24
3. Accounting for the increase in errors due to fewer stations will reduce the ratio to 0.98.

Before we replicate these approaches to debiasing, let's have a closer look at the data from the perspective of Rhines and Huybers (2013).

Comparing the variances¶

First, we look at the variances in the various time periods. We use the data from the Northern hemisphere.

For the baseline, the mean is 0 (as expected) and the variance is close to 1 (it is not exactly 1 because we removed stations with standard deviations with an absolute value of larger than 5. These outliers will have a significant effect on the variance.

For the time period 1981 to 2010, the mean shifts to the right and the variance increases, too. The variance is noticeably different from Rhines and Huybers (2013) because of the removal of outliers.

If we extend the later period to all available data beyond the baseline, the mean continues to shifts to the right and the variance also increases.

We can also have a look at how the mean and the variance shifts over the decades.

We see that the deviation of from the baseline temperature (left chart) changes from negative to positive values. There are two noticeable spikes in the data. One around 1890 and the other around 1930. Browsing Wikipedia reveals that at least for the United States, these spikes correspond to documented heat waves in 1896 and 1936. Not that these explain the spikes but at least the data is consistent with the events. (Remember the data shows the summer months only.)

The variance does not have the same visual trend as the mean temperature.

The black line shows the variance for each decade as a ratio of the baseline. The orange line shows the number of stations reporting temperature measurements for each decade. The more stations report data, the lower is the variance. Intuitively, this makes sense and is in line with the third argument, the measurement errors, by Rhines and Huybers (2013).

Now that we have a relatively good understanding of the biases in the data, let's correct the data for the biases.

Removing the sample means independently¶

We generate the additional information for the first correction suggested by Rhines and Huybers (2013), which means that we calculate the standard deviation for the comparison periods (1981 to 2010 and 1981 to today) separately. We then calculate the z-scores using the baseline mean but the standard deviations from the comparison periods.

Calculating the z-scores with the period-specific standard deviations reduced the variance significantly. For the comparison period 1981-2010 we have a variance of 1.23. For the comparison period 1981-today we have a variance of 1.24.|

The shapes of the distributions of temperatures for the three periods look very similar now.

Correcting the data for spatial variance¶

Spatial variance may have been introduced because the temperatures were measured at different spatial locations.

This is the easiest corrections because it is also accounted for. The Hansen et al. (2012) paper uses temperature data in a 250km grid, which means that the measurements of all the stations within a given grid have been averaged.

Here, we use the actual station data. Thus, our analysis does not suffer from this bias. Hence, our corrected variance of 1.24 is the same as reported by Rhines and Huybers (2013) after they have corrected for the spatial variance.

Correcting for the reduced number of stations¶

A decrease in stations that measure temperatures means that the density of stations in each of the 250km grids decreases and in turn, the likelihood of measurement errors increases. Rhines and Huybers (2013) argue that the measurement error is 1 degree Celsius in each grid. They correct the data by multiplying the error with the ratio of stations available in the baseline and comparison periods.

Since we use the fine-granular station data, we can create the subset of stations that report temperatures for both the baseline and the comparison periods. Of course, we also need to apply the correction for the normalization.

We use the data of 3,195 stations. These stations report across the whole timeframe.

The difference in the variance is further reduced but not to the extent that Rhines and Huyber (2013) suggest.

We test whether there is a significant difference in the variance. We use Brown and Forsythe's (174) test for homogeneity of variance.

The difference in variance is significant for both comparison periods. Despite the reduction in the variance due to the bias corrections, the remaining variance in the comparison periods is still significantly different from the baseline period. In conclusion, it seems that Rhines and Huyber (2013) are correct when pointing out the biased variance but Hansen et al. (2012) are still correct in their assessment that both the mean and the variance change.

Interestingly, the question of whether the variance is changing appears to be an open question in the literature. The following article offers a quick summary of the various positions on the issue:

Alexander and Perkins (2013): Debate heating up over changes in climate variability. Environmental Research Letters (8:2013). DOI: 10.1088/1748-9326/8/4/041001

Data Wrangling for further processing¶

The dataset is very large. Let's condense it for further visualization work:

• To display the histograms, we need the z_scores and the decade.
• We only need the data since 1951.

For better visualization, we are only interested in the left value of the interval of each bin.

Final visualization¶

We have removed all the axis labels to let the data speak for itself. Obviously, this is just a quick attempt to replicate the visualization itself. The New Times version has comprehensive annotations that help to interpret the visualization. This requires a visualization framework that is more powerful in developing custom visualizations (e.g., D3.js). I will address this in a future post.

However, when comparing our version with the version in the New York Times, we notice that the shift in the mean is similar but the difference in the general shape of the distributions is far less pronounced.

Final Thoughts¶

Replicating the original visualizations of the Hansen et al. (2012) paper was relatively easy. Correcting the underlying data for the biases as recommended by Rhines and Huyber (2013) was far more difficult because we had to interpret the recommendations. This highlights the importance of publishing code and data (e.g., as Jupyter files). Then, everyone can understand the process from raw data to results to visualization.

The effect of the recommended corrections to the data attenuates the implications of the data. While the shift toward higher temperatures seems to be clear, the debate is still open on whether we also have more extreme weather (temperatures in the tails of the distribution). This highlights one of the biggest challenges in data analytics: without a deep understanding of the data and its limitations, you cannot derive robust implications (The analysis above is just a very superficial attempt!).

Our analysis also highlights that data preparation, data wrangling, and robustness checks take up most of the time and effort. The final result does not communicate the amount of work that went into. Thus, it is very important to be clear about the objectives of a data analytics product. If you just want to explore a dataset, quick and dirty is enough, if you want to inform your paying audience, you need considerable more resources (as evidenced by the New York Times article). If you want to inform political decisions and investments into the fight against climate change, you need the combined earnest work of countless researchers.

Literature¶

• Alexander and Perkins (2013): Debate heating up over changes in climate variability. Environmental Research Letters (8:2013). DOI: 10.1088/1748-9326/8/4/041001
• Hansen, Sato, Ruedy (2012): Perception of Climate Change. Proceedings of the National Academy of Sciences, 109(37) E2415-E2423; DOI: 10.1073/pnas.1205276109
• Rhines and Huybers (2013): Frequent summer temperature extremes reflect changes in the mean, not the variance, Proceedings of the National Academy of Sciences, 110(7) E546; DOI: 10.1073/pnas.1218748110

This project (including the data as of October 10, 2018) is available at https://github.com/mschermann/loaded_climate_dice.