GEOG393 Central Park Meteorology Project
John Hocknell
Load and prepare data acquired from NOAA.
Using NOAA’s Climate Data Online tool, I ordered the Global Summary
of the Month Data for the Central Park station from January 1920 to
December 2019. The metadata documentation can be found here.
The station’s location has been the same since 1920, but I use
unique() on the station ID, name, latitude, longitude, and elevation
just to be sure that they are the same for every month. I also separate
the date column into two new columns for the year and month to make it
easier to group by month later on. Lastly, I check for missing air
temperature values.
# read in NOAA data and check its structure
data <- read.csv("raw/CP_TAVG_1920-2019.csv")
str(data)## 'data.frame': 1200 obs. of 8 variables:
## $ STATION : chr "USW00094728" "USW00094728" "USW00094728" "USW00094728" ...
## $ NAME : chr "NY CITY CENTRAL PARK, NY US" "NY CITY CENTRAL PARK, NY US" "NY CITY CENTRAL PARK, NY US" "NY CITY CENTRAL PARK, NY US" ...
## $ LATITUDE : num 40.8 40.8 40.8 40.8 40.8 ...
## $ LONGITUDE : num -74 -74 -74 -74 -74 ...
## $ ELEVATION : num 42.7 42.7 42.7 42.7 42.7 42.7 42.7 42.7 42.7 42.7 ...
## $ DATE : chr "1920-01" "1920-02" "1920-03" "1920-04" ...
## $ TAVG : num 23.4 28.5 39.9 48.1 58.7 69.2 73.8 74.2 68.5 61.1 ...
## $ TAVG_ATTRIBUTES: chr ",0" ",0" ",0" ",0" ...names(data) <- c("station", "name", "latitude", "longitude", "elevation", "date", "tavg", "tavg_attributes")
# the CP station changed location once, so just in case - check that there is one unique value for each of the following
unique(data$station)## [1] "USW00094728"unique(data$name)## [1] "NY CITY CENTRAL PARK, NY US"unique(data$latitude)## [1] 40.77898unique(data$longitude)## [1] -73.96925unique(data$elevation)## [1] 42.7# separate date column to create year and month columns
# also separate attribute column by comma according to metadata - x,y where
# x is the number of days missing that month and y is the dataset source
data <- data %>% separate(date, c("year", "month"), "-") %>% separate(tavg_attributes, c("days_missing", "dataset_source"), ",")
# check for missing air temperature, year, and month values
sum(is.na(data$tavg)) # first check for missing tavg values## [1] 0sum(is.na(data$year))## [1] 0sum(is.na(data$month))## [1] 0# investigate the days missing and dataset sources
unique(data$days_missing) # second check for missing tavg values## [1] "" "2"sum(as.numeric(data$days_missing), na.rm = TRUE) # total days missing throughout entire dataset## [1] 8unique(data$dataset_source)## [1] "0" "Z" "X" "W"Calculate monthly mean and meadian values for the entire dataset.
The results are combined into one data frame to make them easily accessible, showing the monthly mean and median for 1920 to 2019.
# calculate the monthly mean for tavg
monthly_mean <- aggregate(tavg ~ month, data, mean, na.rm = TRUE)[2]
# calculate the monthly median for tavg
monthly_median <- aggregate(tavg ~ month, data, median, na.rm = TRUE)[2]
# combine calculated stats into one data frame by month
monthly_stats <- data.frame(unique(data$month), monthly_mean, monthly_median)
stats_colnames <- c("month", "mean", "median")
names(monthly_stats) <- stats_colnames
monthly_stats## month mean median
## 1 01 32.430 32.20
## 2 02 33.989 34.00
## 3 03 41.716 41.85
## 4 04 52.215 52.30
## 5 05 62.488 62.55
## 6 06 71.445 71.45
## 7 07 76.692 76.55
## 8 08 75.210 75.30
## 9 09 68.499 68.30
## 10 10 57.712 57.55
## 11 11 47.213 47.10
## 12 12 36.913 37.25# remove variables to clean up environment
rm(list = c("monthly_mean", "monthly_median"))Calculate monthly mean and meadian values for both time periods.
As per the project guidelines, the monthly statistics are also to be analyzed for each time period. The first time period is from 1920 to 1969 and the second time period is from 1970 to 2019.
# calculate the monthly stats for period 1 (1920-1969) and put into a data frame
period1_mean <- aggregate(tavg ~ month, data[data$year < 1970, ], mean, na.rm = TRUE)[2]
period1_median <- aggregate(tavg ~ month, data[data$year < 1970, ], median, na.rm = TRUE)[2]
period1_monthlies <- data.frame(unique(data$month), period1_mean, period1_median)
names(period1_monthlies) <- stats_colnames
period1_monthlies## month mean median
## 1 01 32.318 31.95
## 2 02 32.946 33.40
## 3 03 40.826 40.25
## 4 04 51.156 50.50
## 5 05 61.898 61.60
## 6 06 71.244 71.30
## 7 07 76.202 76.20
## 8 08 74.476 74.35
## 9 09 68.214 68.25
## 10 10 57.900 57.40
## 11 11 46.692 46.40
## 12 12 35.648 35.75# calculate the monthly stats for period 2 (1970-2019) and put into a data frame
period2_mean <- aggregate(tavg ~ month, data[data$year > 1969, ], mean, na.rm = TRUE)[2]
period2_median <- aggregate(tavg ~ month, data[data$year > 1969, ], median, na.rm = TRUE)[2]
period2_monthlies <- data.frame(unique(data$month), period2_mean, period2_median)
names(period2_monthlies) <- stats_colnames
period2_monthlies## month mean median
## 1 01 32.542 33.05
## 2 02 35.032 35.50
## 3 03 42.606 42.35
## 4 04 53.274 53.45
## 5 05 63.078 63.05
## 6 06 71.646 71.70
## 7 07 77.182 77.40
## 8 08 75.944 75.75
## 9 09 68.784 68.30
## 10 10 57.524 57.75
## 11 11 47.734 47.85
## 12 12 38.178 38.40# remove variables to clean up environment
rm(list = c("period1_mean", "period1_median", "period2_mean", "period2_median"))
# test to check if the mean of both periods is roughly equal to the mean of all years (should all be true if calculated correctly)
round(monthly_stats$mean, 3) == round((period1_monthlies$mean + period2_monthlies$mean) / 2, 3)## [1] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUECalculate annual statistics.
Calculate the annual mean and median values for all three monthly tables.
# calculate annual stats for the entire period, put it in the first row
temp <- data.frame("ANN", mean(monthly_stats$mean), mean(monthly_stats$median))
names(temp) <- stats_colnames
monthly_stats <- rbind(temp, monthly_stats)
# calculate the standard deviations for every month over the entire period
monthly_stats$st_dev <- c(NA, aggregate(tavg ~ month, data, sd, na.rm = TRUE)[[2]])
# calculate annual stats for period 1 (1920-1969)
temp <- data.frame("ANN", mean(period1_monthlies$mean), mean(period1_monthlies$median))
names(temp) <- stats_colnames
period1_monthlies <- rbind(temp, period1_monthlies)
# calculate annual stats for period 1 (1970-2019)
temp <- data.frame("ANN", mean(period2_monthlies$mean), mean(period2_monthlies$median))
names(temp) <- stats_colnames
period2_monthlies <- rbind(temp, period2_monthlies)
# remove variables to clean up environment
rm(list = c("temp", "stats_colnames"))Create the three barplots outlined in the project guide (deliveralbe #3).
The plots are formatted to be saved, so they may be cut-off in RPubs.
# simple barplots to investigate the data
# months_lab <- c("ANN", "J", "F", "M", "A", "M", "J", "J", "A", "S", "O", "N", "D")
# # entire period barplot
# barplot(monthly_stats$mean, names.arg = months_lab, main = "Central Park Monthly Mean Temperature 1920-2019", xlab = "Month", ylab = "Degrees F", col = "purple", border = NA, las = 1, ylim = c(0, 80))
# # period 1 barplot
# barplot(period1_monthlies$mean, names.arg = months_lab, main = "Central Park Monthly Mean Temperature 1920-1969", xlab = "Month", ylab = "Degrees F", col = "purple", border = NA, las = 1, ylim = c(0, 80))
# # period 2 barplot
# barplot(period2_monthlies$mean, names.arg = months_lab, main = "Central Park Monthly Mean Temperature 1970-2019", xlab = "Month", ylab = "Degrees F", col = "purple", border = NA, las = 1, ylim = c(0, 80))
# x axis labels for final plots outlined in the project guide deliverables
x_labs <- c("J", "F", "M", "A", "M", "J", "J", "A", "S", "O", "N", "D", "ANN")
# Deliverable #3 Plot #1 from project guide
# monthly-mean statistics and standard deviations around the bars of the means for the entire period
plot1 <- ggplot(monthly_stats, aes(x = month, y = mean)) +
geom_bar(stat = "identity", width = .7, fill = "purple") +
geom_errorbar(aes(ymin = mean - (st_dev/2), ymax= mean + (st_dev/2)), width=.2) +
coord_cartesian(ylim = c(20, 80)) +
labs(title = "Central Park Monthly Mean Temperature and Standard Deviation 1920-2019",
caption = "Annual (ANN) and monthly mean temperatures from Central Park weather station. Lat = 40.77898, Lon = -73.96925. Bars show mean values, and vertical lines\non the bars show the standard deviations.",
tag = "Figure 1.") +
xlab("Month") +
ylab("Degrees F") +
theme(
plot.title = element_text(face = "bold"),
axis.title.x = element_text(face = "bold"),
axis.title.y = element_text(face = "bold"),
plot.caption = element_text(hjust = 0),
plot.tag = element_text(face = "bold")
) +
scale_y_continuous(minor_breaks = seq(20, 80, by = 5)) +
scale_x_discrete(labels = x_labs)
plot1
# save plot 1
ggsave(file = "plots/plot1.png", plot = plot1, width = 10, height = 7, bg = "gray")
# Deliverable #3 Plot #2 from project guide
# side-by-side bar plot comparing the mean and median air temperature of each month for the entire period
# long pivot so that each mean and median for each month is its own row to plot with ggplot
monthlies_longPivot <- monthly_stats[, 1:3] %>% pivot_longer(!month, names_to = "statistic", values_to = "value")
plot2 <- ggplot(monthlies_longPivot, aes(x = month, y = value, fill = statistic)) +
geom_bar(stat = "identity", position = "dodge", width = .7) +
coord_cartesian(ylim = c(20, 80)) +
labs(title = "Central Park Monthly Mean and Median Temperature 1920-2019",
caption = "Annual (ANN) and monthly mean and median temperatures from Central Park weather station. Lat = 40.77898, Lon = -73.96925.",
tag = "Figure 2.") +
xlab("Month") +
ylab("Degrees F") +
theme(
plot.title = element_text(face = "bold"),
axis.title.x = element_text(face = "bold"),
axis.title.y = element_text(face = "bold"),
legend.title = element_blank(),
legend.position = "bottom",
plot.caption = element_text(hjust = 0),
plot.tag = element_text(face = "bold")
) +
scale_y_continuous(minor_breaks = seq(20, 80, by = 5)) +
scale_x_discrete(labels = x_labs) +
scale_fill_manual(values = c("orange", "black"))
plot2
# save plot 2
ggsave(file = "plots/plot2.png", plot = plot2, width = 10, height = 7, bg = "gray")
# # Deliverable #3 Plot #3 from project guide
# side-by-side bar plot of monthly means comparing the earlier and later time period for each month
# combine monthly means from both periods into one table with a new column to track which period the mean is from
means_byPeriod <-
period1_monthlies %>% select(month, mean) %>% mutate(period = "1920-1969") %>%
bind_rows(
period2_monthlies %>% select(month, mean) %>% mutate(period = "1970-2019"))
plot3 = ggplot(means_byPeriod, aes(x = month, y = mean, fill = period)) +
geom_bar(stat = "identity", position = "dodge", width = .7) +
coord_cartesian(ylim = c(20, 80)) +
labs(title = "Central Park Monthly Mean Temperature",
caption = "Annual (ANN) and monthly mean temperatures from Central Park weather station. Lat = 40.77898, Lon = -73.96925. Comparison of mean from 1920-1969 versus\n1970-2019.",
tag = "Figure 3.") +
xlab("Month") +
ylab("Degrees F") +
theme(
plot.title = element_text(face = "bold"),
axis.title.x = element_text(face = "bold"),
axis.title.y = element_text(face = "bold"),
legend.title = element_blank(),
legend.position = "bottom",
plot.caption = element_text(hjust = 0),
plot.tag = element_text(face = "bold")
) +
scale_y_continuous(minor_breaks = seq(20, 80, by = 5)) +
scale_x_discrete(labels = x_labs) +
scale_fill_manual(values = c("pink", "lightblue"))
plot3
# save plot 3
ggsave(file = "plots/plot3.png", plot = plot3, width = 10, height = 7, bg = "gray")Conduct a t-test (deliverable #4).
# remove rows with the ANN (annual) statistic before conducting t-test
test_data <- means_byPeriod[-(which(means_byPeriod$month == "ANN")),]
# conduct a paired t-test
# paired because we have pairs of observations (the same months from two different time periods)
paired_ttest <- t.test(mean ~ period, data = test_data, paired = TRUE, conf.level = 0.95)
paired_ttest##
## Paired t-test
##
## data: mean by period
## t = -4.644, df = 11, p-value = 0.000712
## alternative hypothesis: true mean difference is not equal to 0
## 95 percent confidence interval:
## -1.7200934 -0.6139066
## sample estimates:
## mean difference
## -1.167