Wczytanie danych

Wczytanie danych przez funkcję read.csv.

dane <- read.csv("../dane/USGS 05427718 YAHARA RIVER AT WINDSOR, WI.csv",sep=";", header=T)

Wyświetlenie pierwszych 10 wierszy:

kable(dane[1:10,1:22])%>%
  kable_styling() %>%
  scroll_box(width = "100%", height = "400px")

---

Dane zostały przygotowane w excelu. Przygotowanie danych polegało na zamianie wartości ‘A’ na ‘,’ ‘P’ na ‘,’, zamianie numeróW danych na nazwę danych np. nr “158128_00010” oznacza “Temperature, water, degrees Celsius”. Następnie zmodyfikowanie nazw kolumn na nazwy “przyjazne” dla R. Nazwy puste zamieniłem na wartości “-999” aby później łatwo zmienić je w R na wartości NA.

Przygotowanie danych

Aby wartości NA były prawidłowo odczytane przez R zamienię wartości -999 na NA

dane[dane == -999] <- NA

Zobaczmy czy zadziałało.

kable(dane[1:10,1:22])%>%
  kable_styling() %>%
  scroll_box(width = "100%", height = "400px")

---

Działa! To teraz przydałoby się rozdzielić datę na rok, miesiąc cyfrowo, miesiąc słownie, dzień, dzień roku.

dateN <- dane$date
dateN <- as.Date(dateN, format("%Y-%m-%d"))

dateN <- data.frame(dateN = dateN,
                 year = as.numeric(format(dateN, format = "%Y")),
                 month = as.numeric(format(dateN, format = "%m")),
                 month.name = format(dateN, format = "%b"),
                 day = as.numeric(format(dateN, format = "%d")),
                 day.of.year = as.numeric(format(dateN, format = "%j"))
                )
head(dateN)

##        dateN year month month.name day day.of.year
## 1 2011-01-01 2011     1        sty   1           1
## 2 2011-01-02 2011     1        sty   2           2
## 3 2011-01-03 2011     1        sty   3           3
## 4 2011-01-04 2011     1        sty   4           4
## 5 2011-01-05 2011     1        sty   5           5
## 6 2011-01-06 2011     1        sty   6           6

---

Aby połączyć te wartości z danymi trzeba użyć polecenia bind_cols() z pakietu dyplr.

library(dplyr)

dane = bind_cols(dane,dateN)

Wyświetlmy nagłówki kolumn z przygotowanej już tabeli do pracy z danymi.

names(dane)

##  [1] "agency_cd"                                                                         
##  [2] "site_number"                                                                       
##  [3] "date"                                                                              
##  [4] "Discharge.cubic.feet.per.second.mean"                                              
##  [5] "Gage.height.feet.mean"                                                             
##  [6] "Stream.water.level.elevation.mean"                                                 
##  [7] "Temperature.water.max"                                                             
##  [8] "Temperature.water.min"                                                             
##  [9] "Temperature.water.mean"                                                            
## [10] "Precipitation.inch."                                                               
## [11] "Phosphorus.water.mean"                                                             
## [12] "Phosphorus.water.unfiltered.pounds.per.day.of.phosphorus.mean"                     
## [13] "Orthophosphate.milligrams.per.liter.as.phosphorus.mean"                            
## [14] "Orthophosphate.water.dissolved.pounds.per.day.mean"                                
## [15] "Ammonia..NH3.NH4...water..dissolved.per.day.nitrogen.mean"                         
## [16] "Ammonia..NH3.NH4...water.filtered.milligrams.per.liter.as.nitrogen.mean"           
## [17] "Ammonia.plus.organic.nitrogen.water.unfiltered.pounds.of.nitrogen.per.day.mean"    
## [18] "Ammonia.plus.organic.nitrogen.water.unfiltered.milligrams.as.nitrogen.per.day.mean"
## [19] "Nitrate.plus.nitrite.water.dissolved.pounds.of.nitrogen.per.day.mean"              
## [20] "Suspended.sediment.short.tons.per.day.mean"                                        
## [21] "Suspended.sediment.concentration.milligrams.per.liter.mean"                        
## [22] "Nitrate.plus.nitrite.water.filtered.milligrams.per.liter.as.nitrogen.Mean"         
## [23] "dateN"                                                                             
## [24] "year"                                                                              
## [25] "month"                                                                             
## [26] "month.name"                                                                        
## [27] "day"                                                                               
## [28] "day.of.year"

---

Analiza 1

Na początek zróbmy wykres średniej, minimalnej i maksymalnej temperatury wody w dla poszczególnych miesięcy w 2016 roku

Najpierw trzeba pogrupwać dane po miesiącach, i wyciągnąć minimalną, średnią i maksymalną temperaturę.

d1 = dane[dane$year==2017,]

d1 %>%
  group_by(month) %>%
    summarise(
      water_temp_min = min(Temperature.water.min, na.rm = TRUE),
      water_temp_max = max(Temperature.water.max, na.rm = TRUE),
      water_temp_mean = mean(Temperature.water.mean, na.rm = TRUE)
    ) -> wynik1

Teraz wczytajmy ggplot2

library(ggplot2)

i stwórzmy wykres!

p <- ggplot(wynik1, aes(x = month))
p <- p + geom_line(aes(y = water_temp_min, colour = "min"))
p <- p + geom_line(aes(y = water_temp_max, colour = "max"))
p <- p + geom_line(aes(y = water_temp_mean, colour = "mean"))
p <- p + scale_colour_manual(values = c("blue", "red","green"))
p <- p + labs(y = "water temperature [°C]",
              x = "month number",
              colour = "Parameter")
p <- p + theme(legend.position = c(0.9, 0.9))

wykres temperatur wody rok 2017

Analiza 2

Po takiej rozgrzewce myślę, że można zrobić coś w miarę prostego w plotly.
Zróbmy wykres przepływu wody w \(m^3/s\) dla roku 2018 i dodajmy jako drugi wykres średnią przepływu wody
od roku 2011 do 2018.

Wybierzmy dane dla roku 2018 i obliczmy średnią dla każdego dnia.

d2 = dane[dane$year==2018,]

d2 %>% 
  group_by(day.of.year) %>%
  summarise(
    Streamflow.mean = mean(Discharge.cubic.feet.per.second.mean, na.rm = TRUE)
  ) -> SF2018

W sumie to nie trzeba było tego robić z summarise ale przyda nam się to dla policzenia
średniej przepływu dla każdego dnia w roku z danych od roku 2011 do 2018

dane %>% 
  group_by(day.of.year) %>%
  summarise(
    Streamflow.mean_2011to2018 = mean(Discharge.cubic.feet.per.second.mean, na.rm = TRUE)
  ) -> SF
SF = SF[-366,]
SF = SF[,2]

Ok mamy już dane pogrupowane teraz chcemy aby przepływ był wyrażony w \(m^3/s\) a nie \(ft^3/s\).

Napiszmy funkję, która zmieni nam te wartości. Przyjmijmy, że \(1 ft^3\) to \(0.02831685 m^3\).

ft3_to_m3 <- function(x) {
  m3 <- x*0.02831685
}

Teraz wywołajmy tę funkcję.

SF2018$Streamflow.mean.m3 = ft3_to_m3(x = SF2018$Streamflow.mean)
SF$Streamflow.mean.m3.2011to2018 = ft3_to_m3(x = SF$Streamflow.mean_2011to2018)

Następnie połączmy przygotowane dane.

dane_plot = bind_cols(SF2018,SF)

Na koniec przygotujmy dane do zrobienia wykresu.

data = seq(as.Date("2018-01-01"),length = 365, by = "day")
lista1 = c(dane_plot$Streamflow.mean.m3)
lista2 = c(dane_plot$Streamflow.mean.m3.2011to2018)
df = data.frame(data,lista1,lista2)

Wczytajmy bibliotekę plotly, która zrobi nam wykres.

library(plotly)

Zróbmy wykres!

p <- plot_ly(
  type = "scatter",
  x = df$data, 
  y = df$lista1,
  name = 'stream_flow_2018',
  mode = "lines",
  line = list(
    color = '#00994C'
  )) %>%
  add_trace(
    type = "scatter",
    x = df$data, 
    y = df$lista2,
    name = 'Stream_flow_mean',
    mode = "lines",
    line = list(color = '#FFB266', width = 4, dash = 'dot')
    ) %>%
  layout(
    title = "Streamflow, m^3/s",
    xaxis = list(
      type = 'date',
      tickformat = "%Y-%m-%d"
    ))
p

Analiza 3

Zajmijmy się może opadami. Zróbmy boxplot dla opadów, który będzie pogrupowany po miesiącach.

Przygotujmy więc dane.

d3 = dane[1005:3040,] #dane od początku pomiaru opadów

d3 %>%
  group_by(year, month, month.name) %>%
  summarise(
    Precipitation.inch. = mean(Precipitation.inch., na.rm = TRUE)
  ) -> opady_SR

Teraz zróbmy porządek z miesiącami, ponieważ wiemy, że na wykresach układają się w losowej kolejnosci.

nazwy = opady_SR$month.name[4:15]
ms = as.character(nazwy)

opady_SR$month.name = factor(opady_SR$month.name , levels = ms)

ggplot został już wczytany wcześniej, więc możemy zrobić wykres.

ggplot(opady_SR, aes(y= Precipitation.inch. ,x=month.name)) +
  geom_boxplot()

Analiza 4

Teraz może coś ciekawszego, bardziej reprezentatywnego, zróbmy wykres słupkowy w plotly, który
będzie pokazywał średnie, najmniejsze i największe opady dla każdego miesiąca w 2017 roku.

No to znowu wczytajmy interesujący nas zakres danych.

d7 = dane[dane$year==2017,]

Amerykanie używają zawsze innych jednostek niż my. W Polsce opady podaje się zazwyczaj w milimetrach,
więc napiszmy funkcję zmieniającą cale na milimetry.

1 cal to 25.4 milimetra.

inch_to_mm <- function(x) {
  Precipitation_mm <- x*25.4
}

Wywołajmy ją.

d7$Precipitation_mm = inch_to_mm(x = d7$Precipitation.inch.)

Policzmy same opady, tzn. nie liczmy dni w którch nie padało w ogóle.
W takim razie zamieńmy warość 0 na NA.

d7[d7 == 0] <- NA

Mamy już prawie gotowe dane. Teraz trzeba je pogrupować po miesiącach i pominąć wartości NA.

d7 %>%
  group_by(month, month.name) %>%
  summarise(
    min = min(Precipitation_mm, na.rm = TRUE),
    mean = mean(Precipitation_mm, na.rm = TRUE),
    max = max(Precipitation_mm, na.rm = TRUE)
  ) -> opady_2017

To samo co wcześniej, robimy porządek z miesiącami.

nazwy = opady_2017$month.name[1:12]
ms = as.character(nazwy)

opady_2017$month.name = factor(opady_2017$month.name , levels = ms)

i tworzymy wykres!

p <- plot_ly(opady_2017, x = ~month.name, y = ~min, type = 'bar', name = 'Precipitation (mm) - min', marker = list(color = '#7fa9ee')) %>%
  add_trace(y = ~mean, name = 'Precipitation (mm) - mean', marker = list(color = '#00112b')) %>%
  add_trace(y = ~max, name = 'Precipitation (mm) - max', marker = list(color = '#207c88')) %>%
  layout(xaxis = list(title = "", tickangle = -45),
         yaxis = list(title = ""),
         margin = list(b = 100),
         barmode = 'group')
p

Analiza 5

No dobra ale, żeby nie było, że dni w których nie padało nie są liczone to zróbmy też wykresy dla
tych wartości, gdzie opad był równy “0”.

Użyjemy do tego funkcji facet wrap, gdzie zrobimy wykresy słupkowe dla każdego dnia grupowanego po miesiącach w 2017 roku.

Wczytajmy zakres danych.

d7 = dane[dane$year==2017,]

zamieńmy cale na milimetry.

inch_to_mm <- function(x) {
  Precipitation_mm <- x*25.4
}

d7$Precipitation_mm = inch_to_mm(x = d7$Precipitation.inch.)

W tym kroku trzeba troszkę oszukać system i stworzyć wektor danych z dniami dla każdego miesiąca. Oszukać, bo przecież po co nam daty tylko dla stycznia? Lepiej by było tylko dni podać a nie całą datę.

day_date = c(seq(as.Date("2017-01-01"),length = 31, by = "day"),
             seq(as.Date("2017-01-01"),length = 28, by = "day"),
             seq(as.Date("2017-01-01"),length = 31, by = "day"),
             seq(as.Date("2017-01-01"),length = 30, by = "day"),
             seq(as.Date("2017-01-01"),length = 31, by = "day"),
             seq(as.Date("2017-01-01"),length = 30, by = "day"),
             seq(as.Date("2017-01-01"),length = 31, by = "day"),
             seq(as.Date("2017-01-01"),length = 31, by = "day"),
             seq(as.Date("2017-01-01"),length = 30, by = "day"),
             seq(as.Date("2017-01-01"),length = 31, by = "day"),
             seq(as.Date("2017-01-01"),length = 30, by = "day"),
             seq(as.Date("2017-01-01"),length = 31, by = "day")
             )

Kolejny raz porządek z miesiącami.

nazwy = c("sty","lut","mar","kwi","maj","cze","lip","sie","wrz","paź","lis","gru")
ms = as.character(nazwy)

d7$month.name = factor(d7$month.name , levels = ms)

No i robimy wykresy!

  d7 %>%
  ggplot(aes(x = day_date, y = Precipitation_mm)) +
  geom_bar(stat = "identity", fill = "darkorchid4") +
  facet_wrap(~month.name, ncol = 3) +
  labs(title = "Dzienne opady rok 2017",
       subtitle = "Dane wykreślone po miesiącu",
       y = "Precipitation (mm)",
       x = "Day") + theme_bw(base_size = 15) +
  scale_x_date(date_labels = "%d")

Analiza 6

Na koniec zróbmy wykresy dla wszystkich danych w naszej tabeli. Użyjemy oczywiście facet wrap.

Wczytajmy tylko te dane w których są jakieś właściwości, nie róbmy wykresu np. z miesięcy, bo to bez sensu.

library(tidyr)

d_all = dane[,4:22] #dane z danymi
d_all$year = dane$year

Teraz użyjemy funkcji gather, która Zbiera wiele kolumn i zwija je w pary klucz-wartość, duplikując wszystkie inne kolumny według potrzeb.

d6=gather(d_all, "zmienna", "value", -year)

Dane, które pobrałem, są odczytywane od różenego roku, np. temperatura wody jest odczytywana dopiero
od 3 kwietnia 2016 roku. Dane gdzie są wszystkie wartości to rok 2017 i 2018 dlatego zróbmy wykresy wszystkich wartości tylko dla roku 2017.

rok2017=d6[d6$year==2017,]

Stwórzmy wykresy.

ggplot(rok2017, aes(x=as.numeric(value))) +
  geom_histogram()+ facet_wrap(~as.factor(zmienna),scales = "free")

## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

## Warning: Removed 555 rows containing non-finite values (stat_bin).

---

2019 copyright by Piotr Felczak

---