# Analiza danych hydrologicznych (Rzeka Soła, wodowskaz w Rajczy



Pobieranie 1,01 Mb.
Data17.11.2017
Rozmiar1,01 Mb.

# Analiza danych hydrologicznych (Rzeka Soła, wodowskaz w Rajczy (Beskid śląski)

> dir()


[1] "przeplywy_RAJCZA1983.csv" "przeplywy_RAJCZA1983.txt"

[3] "przeplywy_RAJCZA1983.xls" "przeplywy_SOL1983.xls"

[5] "stany_wody_RAJCZA1983.csv" "stany_wody_RAJCZA1983.txt"

[7] "stany_wody_RAJCZA1983.xls" "stany_wody_SOL1983.xls"


Przepływy Rajcza1983 – Excel – wygląd danych

Dane eksportowane do Rajcza1983.csv wyglądają następująco:


Rzeka ;SOŁA;;;Rok;1983,00;;93,00;Profil;RAJCZA;;(limnigraf);;;;;

Km;75,00;;;A = ;254,00;km2;;P.z. ;481,20;m nad Kr.;;;;;;

Dz.;XI;XII;I;II;III;IV;V;VI;VII;VIII;IX;X;;;;

1;0,84;0,84;3,04;13,00;4,00;17,80;3,42;1,40;1,48;0,99;0,26;0,63;;;;

2;0,78;0,84;3,04;7,80;3,42;19,00;2,84;1,40;1,10;0,99;0,56;0,63;;;;

3;0,78;0,84;3,04;7,80;3,42;26,00;32,10;1,13;1,10;0,99;0,56;0,56;;;;

4;0,72;0,84;3,04;4,68;3,42;16,60;8,84;1,13;1,00;0,99;0,63;0,63;;;;

5;0,72;0,84;26,00;4,00;2,84;12,00;8,84;1,00;0,90;1,15;0,70;1,31;;;;

6;0,72;0,60;12,00;4,00;2,84;10,90;6,96;4,00;0,86;1,15;0,70;0,84;;;;

7;0,72;1,04;19,00;4,00;23,20;10,90;6,96;2,32;0,82;1,31;0,63;0,70;;;;

8;0,60;1,04;17,80;3,42;23,20;10,90;4,68;1,40;0,79;1,15;0,56;0,70;;;;

9;0,60;1,04;17,80;3,13;26,00;10,90;4,00;1,26;0,76;1,15;0,56;0,70;;;;

10;0,60;1,04;6,96;3,42;26,00;23,20;7,80;1,13;0,73;1,15;0,56;4,40;;;;

11;0,60;2,48;4,68;3,42;26,00;16,60;6,12;1,40;0,73;0,99;0,56;1,31;;;;

12;0,60;3,80;12,00;2,84;19,00;23,20;4,68;1,26;0,86;0,99;0,63;2,28;;;;

13;0,60;2,92;8,84;2,84;7,80;13,00;4,00;1,13;6,06;0,99;0,70;1,64;;;;

14;0,60;2,26;6,96;2,13;7,80;8,84;2,84;1,00;1,29;0,99;0,70;1,64;;;;

15;3,80;1,60;4,00;1,30;4,68;6,12;2,84;1,00;18,60;1,15;0,70;0,84;;;;

16;3,80;1,60;4,00;1,30;4,68;7,80;2,09;2,09;5,22;0,99;0,63;0,84;;;;

17;1,60;33,60;3,52;1,30;6,96;8,84;2,09;2,55;3,72;0,91;0,63;0,77;;;;

18;1,60;19,00;24,60;1,30;12,00;7,80;2,09;2,55;1,80;0,91;0,70;16,70;;;;

19;1,60;13,00;7,80;1,30;24,60;7,80;1,40;3,13;1,48;0,91;0,99;4,44;;;;

20;6,02;7,80;6,12;1,30;26,00;8,84;1,40;48,30;1,31;0,56;0,70;2,76;;;;

21;3,80;4,00;4,00;1,30;26,00;8,84;1,40;14,20;1,31;0,56;0,70;1,64;;;;

22;1,84;4,00;6,12;1,20;33,60;10,90;1,40;6,00;1,31;0,56;0,63;1,64;;;;

23;1,60;5,36;13,00;1,20;24,60;10,90;1,40;4,16;1,48;0,56;0,56;1,31;;;;

24;1,60;4,00;7,80;1,20;12,00;8,84;1,40;3,20;1,48;0,56;0,56;0,90;;;;

25;1,60;3,52;4,68;1,00;13,00;6,96;2,32;2,66;1,31;0,56;0,56;0,99;;;;

26;1,60;2,56;4,68;1,10;29,00;6,96;1,86;2,48;1,16;0,56;0,77;0,84;;;;

27;1,60;2,56;2,00;1,19;14,20;4,68;2,32;2,30;1,15;0,56;0,77;0,99;;;;

28;1,60;3,52;17,40;3,42;7,80;4,00;2,09;3,92;0,99;0,56;0,70;0,84;;;;

29;0,94;3,04;51,00;;8,84;3,42;1,86;4,72;0,99;0,56;0,63;0,84;;;;

30;0,94;3,04;38,60;;17,80;3,42;1,40;2,05;0,99;0,56;0,63;0,77;;;;

31;;2,56;23,20;;17,80;;1,40;;1,31;0,56;;0,70;;;;


#czytamy dane (po korekcie plikow *.csv - korekta polega na wyrzuceniu czterech ; na końcu linijki i pozostawieniu jedynie danych dla poszczególnych miesięcy, możemy to zrobić edytorem ASCII, np. Notepad++)

Przetransformowane dane zapisano jako „przepływy_Rajcza1983.txt”:

1;0,84;0,84;3,04;13,00;4,00;17,80;3,42;1,40;1,48;0,99;0,26;0,63

2;0,78;0,84;3,04;7,80;3,42;19,00;2,84;1,40;1,10;0,99;0,56;0,63

3;0,78;0,84;3,04;7,80;3,42;26,00;32,10;1,13;1,10;0,99;0,56;0,56

4;0,72;0,84;3,04;4,68;3,42;16,60;8,84;1,13;1,00;0,99;0,63;0,63

5;0,72;0,84;26,00;4,00;2,84;12,00;8,84;1,00;0,90;1,15;0,70;1,31

6;0,72;0,60;12,00;4,00;2,84;10,90;6,96;4,00;0,86;1,15;0,70;0,84

7;0,72;1,04;19,00;4,00;23,20;10,90;6,96;2,32;0,82;1,31;0,63;0,70

8;0,60;1,04;17,80;3,42;23,20;10,90;4,68;1,40;0,79;1,15;0,56;0,70

9;0,60;1,04;17,80;3,13;26,00;10,90;4,00;1,26;0,76;1,15;0,56;0,70

10;0,60;1,04;6,96;3,42;26,00;23,20;7,80;1,13;0,73;1,15;0,56;4,40

11;0,60;2,48;4,68;3,42;26,00;16,60;6,12;1,40;0,73;0,99;0,56;1,31

12;0,60;3,80;12,00;2,84;19,00;23,20;4,68;1,26;0,86;0,99;0,63;2,28

13;0,60;2,92;8,84;2,84;7,80;13,00;4,00;1,13;6,06;0,99;0,70;1,64

14;0,60;2,26;6,96;2,13;7,80;8,84;2,84;1,00;1,29;0,99;0,70;1,64

15;3,80;1,60;4,00;1,30;4,68;6,12;2,84;1,00;18,60;1,15;0,70;0,84

16;3,80;1,60;4,00;1,30;4,68;7,80;2,09;2,09;5,22;0,99;0,63;0,84

17;1,60;33,60;3,52;1,30;6,96;8,84;2,09;2,55;3,72;0,91;0,63;0,77

18;1,60;19,00;24,60;1,30;12,00;7,80;2,09;2,55;1,80;0,91;0,70;16,70

19;1,60;13,00;7,80;1,30;24,60;7,80;1,40;3,13;1,48;0,91;0,99;4,44

20;6,02;7,80;6,12;1,30;26,00;8,84;1,40;48,30;1,31;0,56;0,70;2,76

21;3,80;4,00;4,00;1,30;26,00;8,84;1,40;14,20;1,31;0,56;0,70;1,64

22;1,84;4,00;6,12;1,20;33,60;10,90;1,40;6,00;1,31;0,56;0,63;1,64

23;1,60;5,36;13,00;1,20;24,60;10,90;1,40;4,16;1,48;0,56;0,56;1,31

24;1,60;4,00;7,80;1,20;12,00;8,84;1,40;3,20;1,48;0,56;0,56;0,90

25;1,60;3,52;4,68;1,00;13,00;6,96;2,32;2,66;1,31;0,56;0,56;0,99

26;1,60;2,56;4,68;1,10;29,00;6,96;1,86;2,48;1,16;0,56;0,77;0,84

27;1,60;2,56;2,00;1,19;14,20;4,68;2,32;2,30;1,15;0,56;0,77;0,99

28;1,60;3,52;17,40;3,42;7,80;4,00;2,09;3,92;0,99;0,56;0,70;0,84

29;0,94;3,04;51,00;;8,84;3,42;1,86;4,72;0,99;0,56;0,63;0,84

30;0,94;3,04;38,60;;17,80;3,42;1,40;2,05;0,99;0,56;0,63;0,77

31;;2,56;23,20;;17,80;;1,40;;1,31;0,56;;0,70
Podobnie potraktowano stany wody

Teraz możemy wczytać dane:

dane1 = read.table( "stany_wody_RAJCZA1983.txt",sep=";",dec=",")

dane2 = read.table( "przepływy_RAJCZA1983.txt",sep=";",dec=",")


# Ogladamy dane1, jak nazywaja sie kolumny?

head(dane1)

V1 V2 V3 V4 V5 V6 V7 V8 V9 V10 V11 V12 V13

1 1 194 200 206 230 210 238 208 200 202 200 194 195

2 2 193 200 206 220 208 240 206 200 200 200 194 195

3 3 193 200 206 220 208 250 258 198 200 200 194 194

4 4 192 194 206 212 208 236 222 198 199 200 195 195

5 5 192 194 250 210 206 228 222 197 198 201 196 202

6 6 193 190 228 210 206 226 218 210 197 201 196 198

# pierwsza kolumna jest zbedna


# chcemy polaczyc kolumny w jeden wektor
stany = t(dane1[1]) # wpisujemy pierwsza kolumne,używamy funkcji t(), która zamienia kolumny na wiersze
for (i in 2:13) stany = c(stany, t(dane1[i])) # do wektora stany dokładamy w petli

# kolejne kolumny z dane1


# teraz musimy zrobic porzadek
stany = stany[-(1:31)] # wyrzucamy pierwsze 31 liczb
stany = stany[!is.na(stany)] # wyrzucamy "puste" dane (NA)
# Ta sama operacja na przeplywach
przeplywy = t(dane2[1])

for (i in 2:13) przeplywy = c(przeplywy, t(dane2[i])) # do wektora przeplywy dokladamy w

# petli kolejne kolumny z dane1

# teraz musimy zrobic porządek w wektorze przepływów


przeplywy = przeplywy[-(1:31)] # wyrzucamy pierwsze 31 liczb
przeplywy = przeplywy[!is.na(przeplywy)] # wyrzucamy "puste" dane (NA)

plot(stany,przeplywy, cex=0.4) # możemy narysować zależność przepływów od stanów wody




W hydrologii standardowo modeluje się tę zależność wielomianem drugiego stopnia
Rajcza = lm(przeplywy ~ I(stany^2) + stany + 1) # budujemy model Rajcza
summary (Rajcza) # sprawdzamy jak to się udało
Call:

lm(formula = przeplywy ~ I(stany^2) + stany + 1)


Residuals:

Min 1Q Median 3Q Max

-4.5984 -0.1784 0.0127 0.2633 14.5312
Coefficients:

Estimate Std. Error t value Pr(>|t|)

(Intercept) 1.752e+02 8.652e+00 20.26 <2e-16 ***

I(stany^2) 5.417e-03 1.761e-04 30.76 <2e-16 ***

stany -1.953e+00 7.832e-02 -24.93 <2e-16 ***

---


Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 1.209 on 362 degrees of freedom

Multiple R-squared: 0.9746, Adjusted R-squared: 0.9745

F-statistic: 6946 on 2 and 362 DF, p-value: < 2.2e-16
# Chcemy narysować wyniki modelu
stany_model=data.frame(stany=seq(190,280,by=0.5)) # generujemy ramke danych z kolumną

# o nazwie „stany”

# Teraz obliczamy z modelu przepływy dla tej gęstej sieci punktów
przeplywy_model = predict( lm(przeplywy ~ I(stany^2) + stany + 1),new=stany_model)
lines(t(stany_model), przeplywy_model, col="green") # stany_model to kolumna liczb,

# musimy ją zamienić na wektor



# przy użyciu funkcji t()

Mamy model umożliwiający przewidywać przepływy przy znanych stanach wody

modele nieliniowe - świat, niestety, jest nieliniowy
nie dysponujemy jednym dobrym narzedziem do budowy takich modeli

Istnieje wiele różnych metod dopasowywania modeli nieliniowych


Jednym z takich narzędzi jest używana w R funkcja nls
nls - nonlinear least-squares (nieliniowa metoda najmniejszych kwadratów)
Użycie:
model = nls(zm_zal ~ wzór zawierajacy zm_niez i parametry,

+ start=list(odgadniete wartosci poczatkowe parametrow))


UWAGA! Funkcja nls nie może być używana do „bezbłędnych” danych

Ma to znaczenie gdyby ktoś użył sztucznie generowanych danych z użyciem

wzorów matematycznych – do takich danych koniecznie należy dodać zaburzenia

- „szum”.


x = 1:10

y = 2*x + 3 # dokladne dane, linia prosta to szczególny

# przypadek krzywej
yeps = y + rnorm(length(y), sd = 0.01) # dane z szumem

nls(yeps ~ a + b*x, start = list(a = 0.12345, b = 0.54321))

Nonlinear regression model

model: yeps ~ a + b * x

data: parent.frame()

a b


3.009 1.998

residual sum-of-squares: 0.0007632


Number of iterations to convergence: 2

Achieved convergence tolerance: 1.471e-08


## nie potrafi uzbieżnić wynikow dla dokładnych danych!

nls(y ~ a + b*x, start = list(a = 0.12345, b = 0.54321)) # bierzemy dokładne y

Błąd w nls(y ~ a + b * x, start = list(a = 0.12345, b = 0.54321)) :

liczba iteracji przekroczyła maksymalną liczbę 50

niekiedy możemy zdać sie na odgadnięcie parametrów

# początkowych przez nls
x <- -(1:100)/10 # tworzymy wektor x

y <- 100 + 10 * exp(x / 2) + rnorm(x)/10 # do wartości funkcji dokładamy niewielki szum


nlmod <- nls(y ~ Const + A * exp(B * x)) # nie mamy opcji start =

Komunikat ostrzegawczy:

In nls(y ~ Const + A * exp(B * x)) :

No starting values specified for some parameters.

Initializing ‘Const’, ‘A’, ‘B’ to '1.'.

Consider specifying 'start' or using a selfStart model

> nlmod

Nonlinear regression model



model: y ~ Const + A * exp(B * x)

data: parent.frame()

Const A B

99.9909 9.9561 0.4974

residual sum-of-squares: 0.8216
Number of iterations to convergence: 9

Achieved convergence tolerance: 2.634e-06


Jak widać nls bardzo dobrze odtworzył wartości naszej sztucznej funkcji, zaczął od ostrzeżenia, że jako wartości początkowe parametrów użyje liczby 1
Możemy narysować dane i wyniki:
plot(x,y, main = "nls(*), dane, prawdziwa funkcja i model, n=100")
curve(100 + 10 * exp(x / 2), col="red", add = TRUE) # malujemy funkcje
lines(x, predict(nlmod), col="blue") # dorysowujemy wyniki modelu

# przyklad drugi


q = seq(0,2,by=0.2) # tworzymy wektor zm. niezależnej

z = sin(3*pi*q)*exp(2.5*q) # obliczamy wartosci



plot(q,z)

z = sin(3*pi*q)*exp(2.5*q)+0.1*rnorm(length(q)) # "zasmiecamy" wyniki

# dodajac troche szumu

plot(q,z) # rysujemy z jako funkcje q



# teraz dopasowujemy model o nazwie "krzywa"

# widac, ze funkcja jest okresowa i szybko rosnie

krzywa = nls(z ~ sin(A*q + B)* exp(C*q), # <== ogolna postac modelu

start=list(A=1, B=1, C=1)) # <== wartosci poczatkowe

# parametrów

summary (krzywa) # patrzymy co wyszlo
Formula: z ~ sin(A * q + B) * exp(C * q)
Parameters:

Estimate Std. Error t value Pr(>|t|)

A 9.418e+00 4.891e-03 1925 <2e-16 ***

B -1.255e+01 9.519e-03 -1319 <2e-16 ***

C 2.501e+00 7.966e-04 3139 <2e-16 ***

---


Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 0.132 on 8 degrees of freedom
coef(krzywa) # mozemy wyswietlic znalezione parametry modelu

A B C


9.417534 -12.553562 2.500630

points(q,fitted(krzywa),pch=3,col="red") # wyswietlamy wartosci modelu

# w punktach danych

# funkcja fitted pozwala znalezc wartosci uzyskiwane z



# modelu w punktach q

# teraz zaczynamy rysowac uzyskana krzywa

# dorysujemy do rysunku poprzez uzycie komendy lines,

# ktora laczy odcinkami linii prostej podane punkty


nq = seq(0,2,by=0.05) # bierzemy gesta siatke punktow
lines(nq,predict(krzywa,newdata=list(q=nq)),col="red")

# wartosci funkcji (modelu) w punktach nq obliczamy

# stosujac funkcje predict(krzywa,newdata=list(q=nq))

# jej argumenty to nazwa modelu i lista definiujaca nowe wartosci

# zmiennej niezaleznej modelu (tutaj "q")

# Dla porownania mozemy dorysowac oryginalna, tzn. pozbawiona szumow

# funkcje
fun = function(x) sin(3*pi*x)*exp(2.5*x)

curve(fun,0,2,add=TRUE, col="green")


Linia zielona – dokładna funkcja

Linia czerwona wyniki modelu
Szukanie zer funkcji i rozwiązywanie równań nieliniowych
Do takich działań służy funkcja uniroot()
Jej argumenty to


  1. badana funkcja (nazwa lub wzór)

  2. granice, w których mieści się szukane miejsce zerowe funkcji

definiowane jako lower=dolna granica przedziału,

upper= górna granica przedziału,

lub równoważnie: interval=c(lower, upper)

3) tolerance=dokładnośc z jaką ma być wyznaczony pierwiastek


Przykład:


x = seq(-5,6,by=0.5)

y = 2.4523*x^2 + 1.7252*x - 11.3475 # y jest funkcją x

plot(x,y)

abline(h=0) # możemy ocenić gdzie znajdują się miejsca zerowe funkcji y



Nasza funkcja to wielomian, x i y to tabelka wartości tej funkcji. Do zbudowania samej funkcji możemy użyć pakietu polynom


library(polynom)

p = polynomial(coef=c(-11.3475,1.7252,2.4523)) # tworzymy wielomian o zadanych

#wspólczynnikach

p=as.function(p) # zamieniamy go na funkcję


x1 = uniroot(p,lower=-3, upper=-2) # x1 leży gdzieś pomiędzy -3 a -2

x1

$root



[1] -2.531413
$f.root

[1] -0.0002312473


$iter

[1] 4
$estim.prec

[1] 6.103516e-05
x2 = uniroot(p,lower=1, upper=2) # teraz szukamy drugiego pierwiastka

x2

x2



$root

[1] 1.827933


$f.root

[1] 1.049456e-05


$iter

[1] 4
$estim.prec

[1] 6.103516e-05
porównujemy z danymi "analitycznymi"

a=2.4523


b=1.7252

c=-11.3475

delta=b^2-4*a*c
x1a = (-b + sqrt(delta))/(2*a)

x2a = (-b - sqrt(delta))/(2*a)

x1a

[1] 1.827932



x2a

[1] -2.531434

x1a-x2$root # odejmujemy pierwiastki dodatnie

[1] -9.816735e-07

x2a-x1$root # odejmujemy pierwiastki ujemne (numeracja wyników się różni!)

[1] -2.163127e-05


# wyniki sa dokładne na 5-tym miejscu po przecinku
Metodyka pracy przy szukaniu zer funkcji:

  1. Narysować wykres funkcji,

  2. Odczytać przybliżone położenia pierwiastków

  3. Użyć funkcji uniroot do uściślenia rozwiązań w poszczególnych przedziałach

Rozwiązywanie równań nieliniowych:


Ogólna postać równania:
f(x) = stała lub funkcja
Tworzymy funkcję g(x) = f(x)-prawa strona równania i szukamy jej zer.

# równanie exp(t) = 4t^2

# przenosimy wszystkie wyrazy na lewa strone,

# zamieniamy równanie na definicje funkcji f5

f5 = function(t) exp(t)-4*t^2
curve(f5,0,1) # zawsze trzeba sporzadzic wykres

abline(h=0) # i znaleźć przybliżone położenie zera funkcji



uniroot(f5,interval=c(0.7,0.8)) # szukamy rozwiązania

$root

[1] 0.7148075


$f.root

[1] -5.727013e-06


$iter

[1] 3
$estim.prec



[1] 6.103516e-05
W rozdziale 8 podręcznika Carline Soetaert „Using R for scientific computing” można znaleźć przykłady użycia funkcji uniroot do rozwiązywania bardziej skomplikowanych zadań z chemii morza.



©operacji.org 2017
wyślij wiadomość

    Strona główna