Przetwarzanie i wizualizacja danych w R

Opisana niżej analiza zawiera faktyczne elementy analizy danych pomiarowych, będzie stanowiła tło do wykonania projektu zaliczeniowego.

1 Analiza danych z masztu do pomiaru wiatru 1-dniowych

1.1 Zadanie - Import danych dane_w1.csv

Zaczniemy od pobrania ze strony przedmiotu i zaimportowania pliku danych csv o nazwie dane_w1.csv. Uwaga na ustawienie nagłówków jako nazw zmiennych oraz odpowiedniego separatora kolumn i kropki dziesiętnej.

1.2 Opis instalacji

Zaimportowany zbiór danych zawiera rzeczywiste dane pomiarowe z jednej doby pomiaru dla 140 metrowego masztu pomiarowego. Dla porównania, można zajrzeć na listę najwyższych budynków w Krakowie.

Przykładowy masz pomiarowy, źródło: Wikipedia

Pod linkiem Anemometr wygląd przykładowych urządzeń pomiarowych (strona producenta)

Podstawowe informacje techniczne, czyli co mamy w tych danych:

• Rejestrator danych zapisuje dane na karcie pamięci co 10 minut. Pojedynczy rekord (wiersz, fragment na screenie niżej) w pliku danych opisuje 5 parametrów dla każdego urządzenia: Avg - średnia z 10 minut pomiaru, Min i Max - wartości minimalne i maksymalne w ciągu 10 minut, StdDev - odchylenie standardowe, Count - liczba odczytów = 600 (czyli technicznie odczyt następuje co sekundę).

Taki pojedynczy rekord dla urządzenia “Anemometer1” ma postać:

• w przypadku Count różnego od 600 możemy mieć do czynienia z usterką urządzenia pomiarowego, problemów z zasilaniem lub innych problemów technicznych

• w przypadku niskiego St.Dev możemy mieć do czynienia z ciszą lub zamarznięciem urządzenia (np. w okresie zimowym).

• Urządzenia opisane jako “anemometer” i/lub “wind speed” - służą do pomiaru prędkości wiatru [m/s].

• Urządzenia opisane jako “wind vane” i/lub “wind direction” - służą do pomiaru kierunku wiatru [w stopniach].

• czynniki atmosferyczne jak temperatura [w st. Celcjusza], ciśnienie (barometric air pressure) [chyba w hPa], wilgotność względna powietrza [w %]

• parametry rejestratora typu napięcie [w V]

• inne parametry (mniej istotne).

Istotne informacje:

• wysokości zamontowanych anemometrów
  Anemometr1, 140m
  Anemometr2, 140m
  Ultrasonic, 137m
  Anemometr3, 112m
  Anemometr4, 84m

• wysokości zamontowanych vane’ów:
  Ultrasonic, 137m
  Wind Vane 1, 137m
  Wind Vane 2, 81m

1.3 Zadanie - profil pionowy wiatru

Utworzyć ramkę danych zawierającą 2 kolumny: wysokość anemometru (dane wyżej) oraz średnią dobową (średnia kolumn Avg - policzyć trzeba) - 5 wierszy, bo 5 urządzeń mamy. Narysować wykres zależności średniej prędkości wiatru od wysokości (tzw profil pionowy wiatru) - polecenie typu plot(wysokosc, srednia) powinno wystarczyć choć można użyć innego pakietu. Dodać opisy osi, tytuł wykresu.

Ramka i wykres powinny przypominać te poniżej (uwaga: dane wysokości mamy w informacji o maszcie pomiarowym wyżej, a dane prędkości należy samodzielnie wyliczyć, a nie przepisywać cyferki).

##   wysokosc predkosc
## 1      140 3.611029
## 2      140 3.612678
## 3      137 3.570795
## 4      112 3.569185
## 5       84 3.457780

1.4 Zadanie - róża wiatrów

Zainstalować pakiet o nazwie openair. Narysować różę wiatrów, poniżej gotowe polecenie.

library("openair") # ładowanie biblioteki
#?windRose #pomoc
windRose(dane_w1, ws="Ultrasonic.wind_speed.Avg", wd="Ultrasonic.wind_direction.Avg")# właściwe polecenie

2 Dane z pomiaru miesięcznego

2.1 Zadanie - Import danych dane_w2.xlsx

Ponownie zaczniemy od importu danych pobranych ze strony przedmiotu, tym razem plik dane_w2.xlsx. Warto zwrócić uwage na nazwy zmiennych w nagłówku, pierwszą kolumnę można ustawić jako typ danych - date.

2.2 Opis instalacji

Powyższe dane pochodzą z rzeczywistego pomiaru, tym razem dane z całego miesiąca. Jest to inny zestaw pomiarowy niż wcześniej, ale struktura jest podobna: średnia, odchylenie standardowe, minimum i maximum w próbce 10-minutowej.

Urządzenia (plik ze szczegółami jest “tajny”, CH - Channel, kanał):

1. Anemometry:
   CH1, 81.00m; (backup - zapasowy)
   CH4, 80.60m;
   CH5, 60.90m;
   CH6, 40.90m

2. Vane’y:
   CH7, 78.30m;
   CH8, 38.60m

3. Inne:
   CH9 wilgotność
   CH10 ciśnienie
   CH11 temperatura
   CH12 napięcie w rejestratorze

Końcówki w nazwach kolumn odnoszą się do:

• Avg - avarage, średnia z 10 minut pomiaru. Traktujemy ją jako wynik pomiaru.

• SD - odchylenie standardowe

• Min - minimalna wartość w pomiarze 10 minutowym

• Max - maksymalna wartość w pomiarze 10 minutowym

2.3 Zadanie - wykres prędkości wiatru

Narysować wykres (ggplot2 lub plotly) zależności prędkości wiatru dla wszystkich anemometrów (zmienne CH1Avg, CH4Avg, CH5Avg, CH6Avg) od czasu na jednym wykresie. Nadać wykresowi tytuł, opisy osi (w języku polskim).

Wykres powinien przypominać:

2.4 Zadanie - ggplotly

W przypadku użycia ggplot2 w poprzednim zadaniu przetransformować wykres na wersję interaktywną za pomocą ggplotly z pakietu plotly.

Przykładowy wynik:

2.5 Interpretacja wykresów

Wykresy prędkości wiatru w czasie pozwalają na wykrycie nieprawidłowości w działaniu urządzeń pomiarowych:

• Anemometry CH1 i CH4 zostały zamontowane na tej samej wysokości, w związku z czym ich wskazania powinny być zbliżone (niemal identyczne)

• Wraz ze wzrostem wysokości (na ogół) wzrasta prędkość wiatru, zatem jeżeli wysokości zamontowanych urządzeń to:
  CH1, 81.00m; (backup - zapasowy)
  CH4, 80.60m;
  CH5, 60.90m;
  CH6, 40.90m
  to wykresy powinny układać się “monotonicznie”. Co prawda, mogą zdarzać się inwersje np. przy niskich wartościach, jednak może być to wskazówka o awarii urządzenia

• W miesiącach zimowych wykres jak poniżej pozwala na szybkie zidentyfikowanie sytuacji, gdy urządzenia pomiarowe zamarzły - płaska linia na poziomie wartości 0 - tutaj widoczne problemy w końcówce listopada:
  

3 Profil pionowy wiatru

3.1 Wzór/wykładnik Hellmana-Suttona.

Prędkość (średnia) wiatru nie zależy od wysokości w sposób liniowy, jednak podlega wzorowi Hellmana-Suttona:

\[ V(h) = V_1 \cdot \left( \frac{h}{h_1} \right)^w \quad (1) \]

gdzie:
\(h_1\) - bazowa wysokość - znamy wysokości zamontowanych urządzeń, wybieramy jedno z nich
\(V_1\) - średnia prędkość zmierzona na bazowej wysokości - jesteśmy w stanie ją wyliczyć z danych pomiarowych
\(w\) - “nieznany”, ale możliwy do wyznaczenia wykładnik. Jest on charakterystyczny dla danego miejsca pomiaru.

Znając powyższe 3 wartości (po wyliczeniu są to stałe) dostajemy funkcję jednej zmiennej opisującą zależność prędkości wiatru od wysokości.

3.2 Zadanie - wzór Hellmana-Suttona

Przyjmijmy, że dane \(h_1\) i \(V_1\) dotyczą anemometru CH6 oraz znając wartości \(h_2\) (z danych o instalacji) i \(V_2\) (trzeba wyliczyć tą średnią) dla anemometru CH4 powyższe równanie sprowadza się do zależności

\[ V_2=V(h_2) = V_1 \cdot \left( \frac{h_2}{h_1} \right)^w \quad (2) \]

i zawiera jedną niewiadomą \(w\), która jesteśmy wyznaczyć (jakiś logarytm o jakiejś podstawie się przyda, sprawdzić działanie polecenia log w pomocy).

1. Wyprowadzić wzór na współczynnik \(w\) korzystając z równania (2) - “na papierze”

2. Na podstawie danych pomiarowych (średnich) CH4 i CH6 wyznaczyć wartość współczynnika \(w\), (powinno wyjść z grubsza wartość z zakresu 0.2-0.3, tak tylko orientacyjnie).

3. Zbudować funkcję jak we wzorze (1), argumentem powinna być tylko wysokość \(h\), a \(w\), \(h_1\) i \(V_1\) ustawić jako stałe.

4. Wyznaczyć za pomocą funkcji teoretyczną średnią prędkość wiatru na wysokości odpowiadającej CH5 i porównać z danymi pomiarowymi CH5

5. Jakiej, teoretycznie, wartości można się spodziewać na wysokości 130m (jest to często spotykana wysokość montażu gondoli turbiny wiatrowej - tam gdzie zaczepione są śmigła wiatraka).

6. Narysować wykres tej funkcji dla \(h \in [10,130]\)

7. Nanieść na wykres wartości z pomiaru (dla CH4 CH5 i CH6) - wskazówka: może się przydać konstrukcja wykresu z zadania 6.1 z Laboratorium 6

Wykres (wraz z częścią z następnego zadania) powinien przypominać

3.3 Zadanie - prawo logarytmiczne profilu wiatru

Inny sposób opisania zależności prędkości wiatru od wysokości to prawo logarytmiczne:

\[ V(h) = V_1 \cdot \frac{\ln(h/z_0)}{\ln(h_1/z_0)} \]

gdzie:

\(h_1\) - bazowa wysokość
\(V_1\) - średnia prędkość zmierzona na bazowej wysokości
\(z_0\) - współczynnik szorstkości terenu (roughness length). Jest on charakterystyczny dla danego miejsca pomiaru, podobnie jak wcześniej można go wyliczyć mając wartości pomiaru dla dwóch wysokości. Ma wartości zazwyczaj w przedziale 0.0002 dla obszarów wodnych (jeziora, morze) do 1.6 przy gęstej zabudowie.

Zadanie: ponowić obliczenia jak w poprzednim zadaniu, tym razem przy wykorzystaniu prawa logarytmicznego. \(z_0 \approx 1.3\) - tak mniej więcej powinno wyjść. Wykres prawa logarytmicznego nanieść na wykres z prawa Hellmanna.

Więcej informacji:

https://www.intechopen.com/chapters/17121

3.4 Wzorki, jakby ktoś nie umiał wyprowadzić:

Mając znane dwie wysokości \(h_1, h_2\) i policzone średnie prędkości \(V_1, V_2\) możemy wyliczyć

\[ w = \log_{\frac{h_2}{h_1}}{\frac{V_2}{V_1}} \]

albo równoważnie korzystając z własności logarytmu:

\[ w= \frac{\ln V_2-\ln V_1}{\ln h_2-\ln h_1} \]

oraz

\[ z_0=\left(\frac{h_1^{V_2/V_1}}{h_2}\right)^{\frac{1}{V_2/V_1-1}} \]

albo

\[ z_0=\exp\left(\frac{h_1^w \cdot \ln h_2-h_2^w \cdot \ln h_1}{h_1^w-h_2^w}\right) \]

4 Profil poziomy wiatru

4.1 Zadanie - różnica wskazań vane’ów

Instalacje nie są idealne, każde “skrzywienie” może wpłynąć na pomiar kierunku wiatru.

• Narysować wykres wskazań vane’ów (wiatrowskazów) CH7 i CH8 w zależności od czasu. Powinna być widoczna różnica we wskazaniach.

• Narysować wykres wartości bezwzględnej różnicy wskazań vane’ów CH7 i CH8. Wskazówka: można obliczyć tą różnicę i zapisać jako dodatkowy element ramki jak w Laboratorium 7, pakiet dyplr, zadanie z powierzchnią płatka

• Ocenić średnią różnicę między wskazaniami vane’ów, biorąc pod uwagę że mamy do czynienia z kątami z zakresu od 0 do 360 stopni. Przykładowo, kierunki mające wartość 359 i 1st. są właściwie takie sam, choć ich różnica jest bardzo duża. Dodatkowe utrudnienie: przy niskich prędkościach wiatru kierunek może być bardzo zmienny w czasie i między wysokościami (kto puszczał latawce, ten wie), więc wypadałoby wykluczyć z analizy prędkości poniżej 3 m/s dla CH6. Ponownie tutaj można wykorzystać pakiet dyplr, gdzie funkcja mutate w połączeniu case_when w miejsce niezdefiniowanego wyniku wstawia NA

• Narysować ponownie wykres wskazań vane’ów (wiatrowskazów) CH7 i CH8 w zależności od czasu, tym razem z uwzględnieniem wyliczonej różnicy (np. rozdzielając ją po połowie i dodając do CH8Avg i odejmując od CH7Avg tą część różnicy).

4.1.1 Przykładowy wynik:

## [1] "Obliczona poprawka to (przykladowo)"

## [1] 19.12772

4.2 Dopasowanie rozkładu do danych

Dla przypomnienia, z dopasowaniem rozkładu i rozkładami mieliśmy do czynienia na I stopniu:

• Laboratorium 7: rozkłady danych

Istotnym elementem analizy statystycznej danych jest dopasowanie rozkładu do danych empirycznych. Ograniczymy się tutaj do gotowej procedury z pakietu fitdistrplus.

library(fitdistrplus)

4.3 Zadanie - histogram i rozkład prawdopodobieństwa

Do modelowania rozkładu prędkości wiatru najczęściej stosuje się rozkład Weibulla - informacje o rozkładzie Weibulla, który ma dwa parametry: kształt shape i skala scale.

Zajmiemy się analizą prędkości wiatru dla anemometru CH4.

1. Sprawdzić podstawowe informacje dotyczące CH4 poprzez polecenie summary

2. Narysować histogram prędkości wiatru z podziałem co 0.5 m/s. Może się przydać parametr binwidth. Będziemy dopasowywać gęstość i będziemy chcieli nanieść ja na histogram, zatem może przydać się zamiana liczności na częstotliwość występowania poprzez użycie y = after_stat(count / sum(count)) w funkcji aes histogramu.

3. Za pomocą pakietu fitdistrplus wyestymować parametry rozkładu weibull. Sprawdzić “jakość” dopasowania za pomocą gofstat() - patrz laboratorium z I stopnia.

library(fitdistrplus)
summary(dane_w2$CH4Avg)
dane_rozklad = dane_w2$CH4Avg
parametry=fitdist(dane_rozklad, "weibull")

4. Nanieść wykres gęstości z otrzymanymi parametrami na wykres histogramu. Można wykorzystać geometrię: geom_function(fun = function(x) (dweibull(x, shape = ksztalt, scale = skala)), linewidth = 1, colour = "black") gdzie w miejsce ksztalt i skala należy podać wyestymowane wartości

Wykres powinien przypominać:

UWAGA: możliwe jest też przeskalowanie samej funkcji rozkładu Weibulla poprzez przemnożenie dweibull(x, shape = ksztalt, scale = skala) przez liczbę obserwacji razy szerokość okna histogramu:

UWAGA: druga skala z procentami została uzyskana przez scale_y_continuous(sec.axis = sec_axis(~./(4464), name = "Gęstość rozkładu",labels = scales::percent), name = "Liczba wystąpień")

4.4 Zadanie - parametry rozkładu Weibulla - cd

6. Załadować bibliotekę library("mixdist") . Może być potrzebne jej jednorazowe zainstalowanie.

7. Powyższą bibliotekę wykorzystamy do wyznaczenia średniej (i odchylenia standardowego)
   weibullparinv(shape=???, scale=???, loc = 0),
   więcej na
   https://search.r-project.org/CRAN/refmans/mixdist/html/weibullparinv.html

8. Czy uzyskana jak wyżej średnia zgadza się ze średnią empiryczną? (Średnia empiryczna - zwykła średnia dla CH4)

9. Cel tej analizy to oszacowanie, ile energii można wyprodukować w danym miejscu.
   Pod linkiem
   
   Przystępny kalkulator
   
   dostępny jest kalkulator.
   

   W oknie kalkulatora należy ustawić parametr \(A\) - nasz parametr skali i \(k\) - nasz parametr kształtu.

   
   Poeksperymentować z różnymi modelami turbin.
   

4.4.1 Gęstość powietrza

Gęstość powietrza można pozostawić bez zmian lub wyliczyć na podstawie danych ciśnienia - CH10 w hPa oraz temperatury CH11 w stopniach Celcjusza, korzystając ze wzoru: \[ \rho=\frac{p \cdot M}{R \cdot T} \]

gdzie:

p - ciśnienie w Pa (trzeba przeliczyć, 1hPa=100Pa)
M - efektywna masa molowa, 0.0289 kg/mol dla powietrza
R - uniwersalna stała gazowa równa 8.31446261815324 J/(mol*K)
T - temperatura w skali bezwzględnej (Kelvinach), trzeba przeliczyć
Powinno wyjść około 1.2 kg/m3