Climate intro

Computational Topics

• Reading in data from the web into the R / RStudio environment
• Become familiar with variations in CSV / tabular data formats and how to handle them
• Encountering missing data
• Working with dates and date-time objects
• Plotting timeseries data
• Subsetting, reshaping data
• apply functions

Statistical Topics

• Interpret data visualizations
• Explore noise vs seasonality vs trends
• Understand the use of windowed averages

Possible Final Project topics from this module

• Applying regressions, multivariate regression, and forecasting to these data
• Extending analysis to other dimensions of climate
• exploring regional / spatial trends in climate data
• Synthesizing analysis across multiple data sets

Evidence for Global Climate Change

In this module, we will explore several of the most significant data sources on global climate change. An introduction to these data sources can be found at NASA’s Climate Vital Signs website, http://climate.nasa.gov/vital-signs

We will begin by examining the carbon dioxide record from the Mauna Loa Observatory (pictured in title slide).

Why C02?

Carbon dioxide (CO2) is an important heat-trapping (greenhouse) gas, which is released through human activities such as deforestation and burning fossil fuels, as well as natural processes such as respiration and volcanic eruptions.

Parsing tabular data

One of the most common formats we will interact with is tabular data. Tabular data is often presented in plain text, which is not as simple as it sounds, (as we shall see in a moment). NASA points us to a raw data file maintained by NOAA on one of it’s FTP servers: ftp://aftp.cmdl.noaa.gov/products/trends/co2/co2_mm_mlo.txt .

Parsing tabular data

So, where does this data come from? How does one measure atmospheric CO2 levels anyway?

Data Provenance

Knowing where our data come from and how values are measured is essential to proper interpretation of the results. Data scientists usually speak of raw data and derived data, but the intepretation of these terms is always relative. Typically raw simply means “the data I started with” and derived “the data I produced.” Thus our “raw” data is almost always someone else’s “derived” data, and understanding how they got to it can provide important insights to our analysis. One of the first questions we should ask of any data is “where does it come from?”

Data Provenance

In particular, we usually want to make note of three things:

1. What is the uncertainty in the data?

2. What is the resolution of the data?

3. What do missing values mean?

1. What is the uncertainty in the data?

Almost all measurements come with some degree of uncertainty, or measurement error. Often we will not be able to know this precisely. Rather, we seek a a qualtiative understanding of the measurement process to give us some idea of the relative importance of errors in the measurement process influencing the value. We may later be able to infer a more precise description of measurement error from the data itself, but this will always require assumptions about both the data-generating process itself.

2. What is the resolution of the data?

Derived data often summarize raw data in some way. For instance, global climate data is frequently reported as monthly or even annual averages, even though the raw data may be collected day by day, or even minute by minute. Data may be averaged over space as well as time, such as weather measurements made in at separate stations. Weighted averages and more complex techniques are often used as well.

3. What do missing values mean?

Real world data almost always has missing values. Here, it is important we try to understand why values are missing so we know how to handle them appropriately. If there is a systematic reason behind why data are missing (say, days where snowfall or storms made the weather station inaccessible) they could bias our analysis (underestimating extreme cold days, say). If data are missing for an unrelated reason (the scientist is sick, or the instrument fails) then we may be more justified in simply ommitting the data. Often we cannot know the exact reason certain data are missing and this is just something we must keep in mind as a caveat to our infererence. Frequently our results will be independent of missing data, and sometimes missing data can be accurately inferred from the data that is available.

Measuring C02 levels

So how are atmospheric CO2 levels measured?

Researchers shine an infrared light source of a precise intensity through dry air in a container of precisely controlled volume & pressure, ensuring a consistent number of atoms in the chamber. CO2 obsorbs some of this radiation as it passes through the chamber, and then a sensor on the opposite end measures the radiation it recieves, allowing researchers to calculate the amount obsorbed and infer the CO2 concentration. The data are reported in parts per million (ppm), a count of the number of CO2 molecules per million molecules of dry air. These calculations are calibrated by comparing against chambers that are prepared using known concentrations of CO2. For more information, see NOAA documentation.

Measurement uncertainty:

Importantly, the measurement error introduced here is rather small, roughly 0.2 ppm. As we shall see, many other factors, such as local weather and seasonal variation also influence the measurement, but the measurement process itself is reasonably precise. As we move to other sources of data these measurment errors can become much more significant.

Resolution:

What is the resolution of the CO2 data? Already we see our data are not the actual “raw” measurements the researchers at Mauna Loa read off their instruments each day, but have been reported as monthly averages.

Missing values:

The last column of the data set tells us for how many days that month researchers collected data. We see that they only started keeping track of this information in 1974, but have since been pretty diligent – collecting data almost every day of the month (no breaks for weekends here! What do you think accounts for the gaps? How might you test your hypothesis? Would these introduce bias to the monthly averages? Would that bias influence your conclusion about rising CO2 levels?)

Spatially our Mauna Loa data has no aggregation – the data is collected at only one location. How might the data differ if it were aggregated from stations all over the globe?

Importing Data

In Data Science, 80% of time spent prepare data, 20% of time spent complain about need for prepare data.

— Big Data Borat (@BigDataBorat) February 27, 2013

Importing Data

Our first task is to read this data into our R environment. To this, we will use the read.csv function. Reading in a data file is called parsing, which sounds much more sophisticated. For good reason too – parsing different data files and formats is a cornerstone of all pratical data science research, and can often be the hardest step.

Importing Data

So what do we need to know about this file in order to read it into R?

library("tidyverse")

## Let's try:
co2 <- read_table("ftp://aftp.cmdl.noaa.gov/products/trends/co2/co2_mm_mlo.txt")
co2

# A tibble: 784 x 1
   `# --------------------------------------------------------------------`
                                                                      <chr>
 1                                                  # USE OF NOAA ESRL DATA
 2                                                                        #
 3             # These data are made freely available to the public and the
 4       # scientific community in the belief that their wide dissemination
 5        # will lead to greater understanding and new scientific insights.
 6         # The availability of these data does not constitute publication
 7   # of the data.  NOAA relies on the ethics and integrity of the user to
 8     # insure that ESRL receives fair credit for their work.  If the data
 9       # are obtained for potential use in a publication or presentation,
10      # ESRL should be informed at the outset of the nature of this work.
# ... with 774 more rows

Importing Data

hmm… no luck. Let’s try defining the comment symbol:

co2 <- read_tsv("ftp://aftp.cmdl.noaa.gov/products/trends/co2/co2_mm_mlo.txt",
                comment = "#")
co2

# A tibble: 712 x 1
   `1958   3    1958.208      315.71      315.71      314.62     -1`
                                                               <chr>
 1   1958   4    1958.292      317.45      317.45      315.29     -1
 2   1958   5    1958.375      317.50      317.50      314.71     -1
 3   1958   6    1958.458      -99.99      317.10      314.85     -1
 4   1958   7    1958.542      315.86      315.86      314.98     -1
 5   1958   8    1958.625      314.93      314.93      315.94     -1
 6   1958   9    1958.708      313.20      313.20      315.91     -1
 7   1958  10    1958.792      -99.99      312.66      315.61     -1
 8   1958  11    1958.875      313.33      313.33      315.31     -1
 9   1958  12    1958.958      314.67      314.67      315.61     -1
10   1959   1    1959.042      315.62      315.62      315.70     -1
# ... with 702 more rows

Importing Data

Getting there, but not quite done. Our first row is being interpreted as column names. The documentation also notes that certain values are used to indicate missing data, which we would be better off converting to explicitly missing so we don’t get confused.

co2 <- read.table("ftp://aftp.cmdl.noaa.gov/products/trends/co2/co2_mm_mlo.txt",
                  sep = "", comment = "#",
                  col.names = c("year", "month", "decimal_date", 
                                "average", "interpolated", 
                                "trend", "days"),
                  na.strings = c("-1", "-99.99"))
co2 %>% head()

Importing Data

Alternately, with readr::read_table from tidyverse

#devtools::install_github("tidyverse/readr")
co2 <- read_table("ftp://aftp.cmdl.noaa.gov/products/trends/co2/co2_mm_mlo.txt", 
                  comment = "#",
                  col_names = c("year", "month", "decimal_date", 
                                "average", "interpolated", "trend", "days"),
                  col_types = c("iiddddi"),
                  na = c("-1", "-99.99"))

No luck, see: https://github.com/tidyverse/readr/pull/563

Importing Data

Seems that comment arg is not fully implemented in CRAN version of readr so we must rely on skip to avoid the comment block:

co2 <- read_table("ftp://aftp.cmdl.noaa.gov/products/trends/co2/co2_mm_mlo.txt",
                  col_names = c("year", "month", "decimal_date", 
                                "average", "interpolated", "trend", "days"),
                  col_types = c("iiddddi"),
                  na = c("-1", "-99.99"),
                  skip = 72)
co2 

Success! We have read in the data. Now we’re ready to rock and roll.

Plotting Data with ggplot

Effective visualizations are an integral part of data science, poorly organized or poorly labelled figures can be as much a source of peril as understanding. Nevertheless, the ability to generate plots quickly with minimal tinkering is an essential skill. As standards for visualizations have increased, too often visualization is seen as an ends rather than a means of data analysis. See Fox & Hendler (2011) for more discussion of this.

Plotting Data with ggplot

ggplot(co2, aes(decimal_date, average)) + geom_line()

Plotting multiple series

We often would like to plot several data values together for comparison, for example the average, interpolated and trend co2 data. We can do this in three steps:

1. subsetting the dataset to the columns desired for plotting

co2_sub <- co2 %>%
    select(decimal_date, average, interpolated, trend)
co2_sub %>% head()

Plotting multiple series

2. rearranging the data into a “long” data table where the data values are stacked together in one column and there is a separate column that keeps track of the whether the data came from the average, interpolated, or trend column. Notice by using the same name, we overwrite the original co2_sub

co2_sub <- co2_sub %>%
    gather(series, ppmv, -decimal_date)
co2_sub %>% head()

Plotting multiple series

3. plotting

co2_sub %>%
 ggplot(aes(decimal_date, ppmv, col = series)) + 
  geom_line()

Plotting multiple series

Or we can take advantage of dplyr’s nifty pipping abilities and accomplish all of these steps in one block of code. Beyond being more succinct, this has the added benefit of avoiding creating a new object for the subsetted data.

co2 %>%
  select(decimal_date, average, interpolated, trend) %>%
  gather(series, ppmv, -decimal_date) %>%
  ggplot(aes(decimal_date, ppmv, col = series)) +  geom_line()

What do we see?

Our “Figure 1” shows three broad patterns:

• A trend of steadily increasing CO2 concentration from 1950 to 2015
• Small, regular seasonal oscillation is visible in the data
• Increase appears to be accelerating (convex curve)

Exploring Seasonal Oscillations

• Demonstrate that the periodic behavior truly is seasonal
• What month is the maximum? What is the minimum?
• Is that consistent throughout the data? Can you plot the month the minumum and maximum occurs in each year? (exercise for the reader)
• What do you think could explain the seasonal cycle observed here?

Exploring Seasonal Oscillations

## Seasonal oscillation
co2 %>% 
  group_by(year) %>% 
  ## Note we need month[which.max] since not all years have all 12 months
  summarise(max_month = month[which.max(average)], min_month = month[which.min(average)]) %>%
  gather(id, value, -year) %>%
  ggplot(aes(value, fill=id)) + stat_count()

How could we get months by name instead?

Further exploration:

• has the oscillation amplitude increased over time?

Trend, cycle, or noise?

“Climate is what you expect, weather is what you get”

Present-day climate data is often sampled at both finer temporal and spatial scales than we might be interested in when exploring long-term trends. More frequent sampling can reveal higher-frequency trends, such as the seasonal pattern we observe in the CO2 record. It can also reveal somewhat greater variability, picking up more random (stochastic) sources of variation such as weather patterns.

Trend, cycle, or noise?

To reveal long term trends it is frequently valuable to average out this high-frequency variation. We could spend the whole course discussing ways such averaging or smoothing can be done, but instead we’ll focus on the most common methods you will see already present in the climate data we examine. The monthly record data we analyze here already shows some averaging. How was this performed?

Moving averages

# install.packages("RcppRoll")
library(RcppRoll)
co2 %>%
 mutate(annual = RcppRoll::roll_mean(average, 
                                     n=12L, 
                                     align = "left", 
                                     fill = NA, 
                                     na.rm=TRUE, 
                                     normalize=FALSE)) ->
  co2
head(co2)

Moving averages

co2 %>% ggplot(aes(decimal_date)) + 
  geom_line(aes(y=average), col="blue") + 
  geom_line(aes(y=annual), col="red")

Warning: Removed 11 rows containing missing values (geom_path).