# File-Name:       ufo_sightings.R           
# Date:            2012-02-10                                
# Author:          Drew Conway (drew.conway@nyu.edu)
# Purpose:         Kod do rozdziału 1.  Będziemy tu używać podstawowych funkcji R i ćwiczyć
#                   stosowanie paradygmatów kodowania wykorzystywanych w dalszych ćwiczeniach.
#                   Przećwiczymy wczytywanie, podglądanie i oczyszczanie surowych danych oraz
#                   podstawowe wizualizacje danych. To konkretne studium przypadku opiera się
#                   na danych  o obserwacjach UFO; sprawdzimy, czy w danych widać (i jakie)
#                   trendy okresowe.
# Data Used:       http://www.infochimps.com/datasets/60000-documented-ufo-sightings-with-text-descriptions-and-metada
# Packages Used:   ggplot2, plyr, scales

# All source code is copyright (c) 2012, under the Simplified BSD License.  
# For more information on FreeBSD see: http://www.opensource.org/licenses/bsd-license.php

# All images and materials produced by this code are licensed under the Creative Commons 
# Attribution-Share Alike 3.0 United States License: http://creativecommons.org/licenses/by-sa/3.0/us/

# All rights reserved.

# UWAGA: w przypadku uruchamiania z konsoli R (w trybie interaktywnym) należy wykonać polecenie 'setwd'
# w celu zmiany katalogu roboczego na katalog zawierający plik skryptu.
# Inaczej może dojsć do błędów wczytywania danych i zapisu obrazków!

# Wczytanie bibliotek i danych
library(ggplot2)    # do wszystkich wizualizacji
library(plyr)       # do manipulowania danymi
library(scales)     # do poprawiania formatu danych pod wykresy

# Plik zawiera dane rozdzielane spacjami, więc używamy funkcji 'read.delim' z separatorem
# ustawionym na znak tabulacji. Zmieniamy też dwa ustawienia domyślne: 1) ciągi znaków nie
# mają być konwertowane na typy kategoryzowane; 2) dane nie posiadają wiersza nagłówkowego,
# więc pierwszy wiersz interpretujemy jako wiersz danych.
ufo <- read.delim(file.path("data", "ufo", "ufo_awesome.tsv"),
                  sep = "\t",
                  stringsAsFactors = FALSE,
                  header = FALSE, 
                  na.strings = "")
# Plik jest dość duży, więc wczytywanie może chwilę potrwać

# Podejrzymy ramkę danych
summary(ufo)
head(ufo)

# Na podstawie pliku opisu danych ustawiamy odpowiednio nazwy kolumn (za pomocą funkcji 'names')
names(ufo) <- c("DateOccurred", "DateReported",
                "Location", "ShortDescription",
                "Duration", "LongDescription")

# Aby przetwarzać daty zamienimy je z postaci ciągów znaków w formacie RRRRMMDD na obiekty dat
# języka R za pomocą funkcji 'strptime'

# Ale coś w danych jest nie tak; na razie zignorujemy błędy, usuwając niepoprawnie przetworzone wiersze.
# Wiemy, że ciągi dat mają zawsze 8 znaków, więc odstępstwo od tej reguły to sygnał do zignorowania wiersza.
# Wektor oznaczeń błędnych wierszy zbudujemy za pomocą funkcji 'ifelse'
good.rows <- ifelse(nchar(ufo$DateOccurred) != 8 |
                    nchar(ufo$DateReported) != 8,
                    FALSE,
                    TRUE)
length(which(!good.rows))      # 731 niepoprawnie zapisanych wierszy to niby dużo, ale łącznie jest ich ponad 60K
ufo <- ufo[good.rows, ]        # więc pozbywamy się zaledwie 0.6% łącznej liczby rekordó∑

# Teraz możemy już skonwertować ciągi znaków dat na obiekty dat języka R
ufo$DateOccurred <- as.Date(ufo$DateOccurred, format = "%Y%m%d")
ufo$DateReported <- as.Date(ufo$DateReported, format = "%Y%m%d")

# Przydałoby się utworzyć z kolumny Location osobne kolumny dla miasta i stanu.
# Użyjemy do tego funkcji 'strsplit' z wyrażeniem regularnym.
# Uwaga: nie wszystkie wpisy w kolumnie Location mają postać "Miasto, stan"; 
# użyjemy funkcji 'tryCatch', która w takich przypadkach zwróci parę [NA, NA].
# Potem jeszcze usuniemy spację sprzed nazw miast i stanów za pomocą funkcji 'gsub'.
get.location <- function(l)
{
  split.location <- tryCatch(strsplit(l, ",")[[1]],
                             error = function(e) return(c(NA, NA)))
  clean.location <- gsub("^ ","",split.location)
  if (length(clean.location) > 2)
  {
    return(c(NA,NA))
  }
  else
  {
    return(clean.location)
  }
}

# Funkcja 'lapply' zwraca listę wektorów [City, State]
city.state <- lapply(ufo$Location, get.location)

# Za pomocą funkcji 'do.call' izamienimy listę na macierz N-by-2
location.matrix <- do.call(rbind, city.state)

# Dodajemy dane o mieście i stanie do ramki danych o UFO. Możemy do tego użyć funkcji 'transform'
ufo <- transform(ufo,
                 USCity = location.matrix[, 1],
                 USState = location.matrix[, 2],
                 stringsAsFactors = FALSE)

# W następnym kroku wytniemy obserwacje spoza Stanów Zjednoczonych

# Dla obserwacji w miastach spoza USA wstawiamy NA
ufo$USState <- state.abb[match(ufo$USState, state.abb)]

# Na koniec za pomocą funkcji 'subset' wybierzemy tylko obserwacje z USA i zamienimy nazwy stanów na
# wartości zmiennej skategoryzowanej.
ufo.us <- subset(ufo, !is.na(USState))

# Teraz jesteśmy gotowi do analiz! Zobaczmy, jak wyglądają wstępnie przetworzone dane
summary(ufo.us)
head(ufo.us)

# Powyższe funkcje podsumowujące pokazują, że dane pochodzą z naprawdę dawnych obserwacji (do roku 1440!).
# Sprawdźmy więc, w jakim okresie występuje większość obserwacji. Możemy w tym celu utworzyć histogram
# częstości obserwacji UFO w czasie
quick.hist <- ggplot(ufo.us, aes(x = DateOccurred)) +
  geom_histogram() + 
  scale_x_date(breaks = "50 years", labels = date_format("%Y")) +
  labs(x = "Data obserwacji", y = "Liczba obserwacji")
  
ggsave(plot = quick.hist,
       filename = file.path("images", "quick_hist.pdf"),
       height = 6,
       width = 8)

# Jak widać, w danych jest sporo bardzo starych obserwacji; dla naszych celów ograniczymy się więc do
# incydentów z okresu od 1990 roku do obecnego momentu.
ufo.us <- subset(ufo.us, DateOccurred >= as.Date("1990-01-01"))

# Zobaczmy, jak histogram wygląda teraz
new.hist <- ggplot(ufo.us, aes(x = DateOccurred)) +
  geom_histogram(aes(fill='white', color='red')) +
  scale_fill_manual(values=c('white'='white'), guide="none") +
  scale_color_manual(values=c('red'='red'), guide="none") +
  scale_x_date(breaks = "5 years", labels = date_format("%Y")) +
  labs(x = "Data obserwacji", y = "Liczba obserwacji")

ggsave(plot = new.hist,
       filename = file.path("images", "new_hist.pdf"),
       height = 6,
       width = 8)

# Skoro dane wyglądają tak jak chcemy, zacznijmy je agregować. Użyjemy funkcji 'ddply' z pakietu plyr.
# Ale najpierw utworzymy nową kolumnę podającą rok i miesiąc obserwacji.
ufo.us$YearMonth <- strftime(ufo.us$DateOccurred, format = "%Y-%m")

# Poniższe wyrażenie zwróci liczbę obserwacji UFO z podziałem na rok-miesiąc i stan
sightings.counts <- ddply(ufo.us, .(USState,YearMonth), nrow)

# Zgodnie z oczekiwaniem, istnieje kilka wyróżniających się kombinacji rok-miesiąc-stan, w których w ogóle nie
# obserwowano UFO. Policzymy je jako zero.
# Najpierw utworzymy nowy wektor, posiadający wszystkie pary rok-miesiąc z naszego przedziału czasowego (1990-2010)
date.range <- seq.Date(from = as.Date(min(ufo.us$DateOccurred)),
                       to = as.Date(max(ufo.us$DateOccurred)),
                       by = "month")
date.strings <- strftime(date.range, "%Y-%m")

# Aby uzupełnić brakujące dane z ramki 'sightings.counts' utworzymy osobną ramkę danych z kolumną rok-miesiąc.
states.dates <- lapply(state.abb, function(s) cbind(s, date.strings))
states.dates <- data.frame(do.call(rbind, states.dates),
                           stringsAsFactors = FALSE)

# Za pomocą funkcji 'merge' złączymy posiadane dane o obserwacjach z zerami wpisywanymi dla brakujących pozycji.
# Zauważmy, że wszystkie kolumny uczestniczące w złączeniu musimy podać w wywołaniu, i ustawić parametr 'all'
# na TRUE, co spowoduje dobranie danych uzupełniających (NA) również dla tych pozycji, które nie mają obserwacji w
# ramce źródłowej.
all.sightings <- merge(states.dates,
                       sightings.counts,
                       by.x = c("s", "date.strings"),
                       by.y = c("USState", "YearMonth"),
                       all = TRUE)

# Wystarczy już tylko trochę oczyścić złączoną ramkę danych.
# Ustawiamy nazwy kolumn na znaczące etykiety
names(all.sightings) <- c("State", "YearMonth", "Sightings")

# Zamieniamy NA na zera
all.sightings$Sightings[is.na(all.sightings$Sightings)] <- 0

# Zamieniamy z powrotem pary rok-miesiąc na obiekty dat
all.sightings$YearMonth <- as.Date(rep(date.range, length(state.abb)))

# Zamieniamy oznaczenie stanu na wielkie litery i konwertujemy na wartości zmiennej skategoryzowanej
all.sightings$State <- as.factor(all.sightings$State)

# Sezonowość moglibyśmy wykrywać na wiele różnych sposobów, ale podstawowa metoda sprowadza się do
# wizualnej analizy trendu. Skonstruujemy więc wykres, który pokaże ewentualne trendy osobno dla
# wszystkich 50 stanów USA.

# Na początek tworzymy obiekt ggplot2 i tworzymy w nim warstwę geom, która tutaj jest linią. Warto zaznaczyć, że:
# (1) funkcja facet_wrap() utworzy osobne wykresy dla poszczególnych stanów, rozłożone na siatce 10x5
# (2) funkcja theme_bw() zmieni domyślny styl wykresu ggplot2 z szarego na biały (kwestia gustu)
# (3) funkcja scale_color_manual() ustawi kolor linii wykresu na ciemnoniebieski
# (4) funkcja scale_x_date() przeskaluje oś x jako daty, z wyznaczeniem linii siatki pomocniczej co 5 lat
# (5) funkcje xlab() i ylab() ustawią etykiety osi wykresu
# (6) funkcja opts() ustawi tytuł wykresu
state.plot <- ggplot(all.sightings, aes(x = YearMonth,y = Sightings)) +
  geom_line(aes(color = "darkblue")) +
  facet_wrap(~State, nrow = 10, ncol = 5) + 
  theme_bw() + 
  scale_color_manual(values = c("darkblue" = "darkblue"), guide = "none") +
  scale_x_date(breaks = "5 years", labels = date_format('%Y')) +
  xlab("Lata") +
  ylab("Liczba obserwacji") +
  ggtitle("Liczba obserwacji UFO w kolejnych miesiącach, z podziałem na stany (1990-2010)")

# Zapiszmy wykres w pliku PDF
ggsave(plot = state.plot,
       filename = file.path("images", "ufo_sightings.pdf"),
       width = 14,
       height = 8.5)


# Utworzymy nowy wykres, na którym liczba obserwacji będzie znormalizowana względem populacji danego stanu
state.pop <- read.csv(file.path('data/census.csv'), stringsAsFactors=FALSE)

state.pop$abbs <- sapply(state.pop$State, function(x) state.abb[grep(paste('^', x, sep=''), state.name)])
all.sightings$Sightings.Norm <- sapply(1:nrow(all.sightings), 
    function(i) all.sightings$Sightings[i] / state.pop$X2000[which(state.pop$abbs== all.sightings$State[i])])
    
    
state.plot.norm <- ggplot(all.sightings, aes(x = YearMonth,y = Sightings.Norm)) +
  geom_line(aes(color = "darkblue")) +
  facet_wrap(~State, nrow = 10, ncol = 5) + 
  theme_bw() + 
  scale_color_manual(values = c("darkblue" = "darkblue"), guide = "none") +
  scale_x_date(breaks = "5 years", labels = date_format('%Y')) +
  xlab("Lata") +
  ylab("Liczba obserwacji na mieszkańca (wg spisu z 2000 r.)") +
  ggtitle("Liczba obserwacji UFO w kolejnych miesiącach, z podziałem na stany (1990-2010)")
  
  
# Zapiszmy wykres jako plik PDF
ggsave(plot = state.plot.norm,
     filename = file.path("images", "ufo_sightings_norm.pdf"),
     width = 14,
     height = 8.5)
