R-ohjelmointi.org

Tilastotieteellistä ohjelmointia R-kielellä

Harppauskerroksen määrittäminen

Harppauskerros on vesistössä eri vesimassojen välinen vaihettumisvyöhyke. Esimerkiksi lämpötilan suhteen järveen tai mereen voi muodostua termokliini. Itämeressä on myös tyypillisesti halokliini eli kahden eri suolapitoisuudeltaan erilaisen vesimassa välinen vyöhyke.

Termokliini on helppo määrittää mittaamalla vesistön tietyn pisteen lämpötila pinnasta riittävän syvälle saakka metrin välein. Homman voi tehdä lämpömittarilla tai sopivalla anturilla, joita saa esimerkiksi kaikuluotaimiin. Tämän jälkeen pitää enää piirtää mittaustuloksista käyrä, ja etsiä käyrältä se piste, jossa veden lämpötila muuttuu ensimmäisen kerran pinnalta pohjaa kohti edettäessä vähintään yhden Celsius-asteen metrillä. Kyseinen kohta on termokliinin yläpinta.

Seuraava R-koodi auttaa harppauskerroksen määrittämisessä. Muodostetaan ensin demodata:

# Demodata
# Sarake depth = syvyys
# Sarake temp  = lämpötila
d<-structure(list(depth = 0:25, temp = c(18.1, 18.1, 18, 18, 18, 
17.8, 17.8, 17.6, 17.6, 17.5, 17, 17, 16.9, 16.5, 16.3, 12, 11.8, 
11.7, 10.4, 8.7, 8.5, 8.5, 8.4, 8.4, 8.3, 8.3), yr = c(0, 0, 
-0.100000000000001, 0, 0, -0.199999999999999, 0, -0.199999999999999, 
0, -0.100000000000001, -0.5, 0, -0.100000000000001, -0.399999999999999, 
-0.199999999999999, -4.3, -0.199999999999999, -0.100000000000001, 
-1.3, -1.7, -0.199999999999999, 0, -0.0999999999999996, 0, -0.0999999999999996, 
0), ar = c(0, 0.100000000000001, 0, 0, 0.199999999999999, 0, 
0.199999999999999, 0, 0.100000000000001, 0.5, 0, 0.100000000000001, 
0.399999999999999, 0.199999999999999, 4.3, 0.199999999999999, 
0.100000000000001, 1.3, 1.7, 0.199999999999999, 0, 0.0999999999999996, 
0, 0.0999999999999996, 0, 0)), .Names = c("depth", "temp", "yr", 
"ar"), row.names = c(NA, 26L), class = "data.frame")

Data voidaan tietysti lukea R:ään tiedostostakin. Seuraavan koodipätkän toinen rivi lukee tiedoston sarkainerotellusta tekstitiedostosta. Komento file.choose() antaa mahdollisuuden selata tietokoneen kansiorakennetta, ja valita luettavan tiedoston tiedostolistauksesta suoraan.

Kun data on luettu sisään, loppuosa koodista piirtää lämpötilakäyrän ja merkitsee siihen myös termokliinin sijainnin. Ei haittaa, vaikka aineistossa olisi puuttuvia havaintoja tai havainnot eivät olisi tasavälisesti metrin välein, sillä koodi selviää sellaisistakin tilanteista. Ainoa ehto on, että syvyyden mittayksikkö on metri, ja lämpötilan Celsius-aste.

library(zoo)
d<-read.table(file.choose(), header=TRUE, "\t") 
 
d<-d[order(d$depth),]
d$temp<-na.locf(d$temp)
 
plot(y=rev(d$depth-max(d$depth)), x=d$temp, type="n", las=1, xlab="Lämpötila", ylab="Syvyys", xaxs="i", yaxs="i")
grid(lty=1, col="grey85")
lines(y=rev(d$depth-max(d$depth)), x=d$temp, lwd=2)
box()
 
if(any(diff(d$depth)>1)) {
   de<-diff(d$depth)
   te<-diff(d$temp)
   di<-c(0, te/de)
   d$temp2<-c(di)
   d$yr<-c(d$temp2)
   d$ar<-rev(d$temp2)
   polygon(x=c(min(d$temp), min(d$temp), max(d$temp), max(d$temp)), y=c(-rev(d$depth)[which(rev(d$ar)<=c(-1))[1]], -d$depth[which(d$yr<=c(-1))[1]+1], -d$depth[which(d$yr<=c(-1))[1]+1], -rev(d$depth)[which(rev(d$ar)<=c(-1))[1]]), col="#CC000022", border=NA)
   legend(x="topleft", fill="#CC000022", legend=paste("Termokliinin sijainti (", rev(d$depth)[which(rev(d$ar)<=c(-1))[1]], " - ", d$depth[which(d$yr<=c(-1))[1]+1], " m)", sep=""), bty="n", border="#CC000022")
} else {
   d$yr<-c(0, diff(d$temp))
   d$ar<-rev(c(0, diff(rev(d$temp))))
   polygon(x=c(min(d$temp), min(d$temp), max(d$temp), max(d$temp)), y=c(-rev(d$depth)[which(rev(d$ar)>=c(1))[1]-1], -d$depth[which(d$yr<=c(-1))[1]-1], -d$depth[which(d$yr<=c(-1))[1]-1], -rev(d$depth)[which(rev(d$ar)>=c(1))[1]-1]), col="#CC000022", border=NA)
   legend(x="topleft", fill="#CC000022", legend=paste("Termokliinin sijainti (", d$depth[which(d$yr<=c(-1))[1]-1], " - ", rev(d$depth)[which(rev(d$ar)>=c(1))[1]-1], " m)", sep=""), bty="n", border="#CC000022")
}

Lopputuloksena on esimerkiksi alla olevan kaltainen kuva, jossa vesistön pinta on ylhäällä.

Yllä olevassa ratkaisussa on kuitenkin useita hankalia piirteitä. Ensinnäkään se ei toimi, jos mittauksia ei ole tehty juuri metrin välein, mikä on usein tilanne, jos mittaukset on tehty käsin, mutta sama ongelma on myös ctd-laitteistolla tehdyissä mittausissa. Käsin tehdyissä mittauksissa mittaväli on usein epäsäännöllinen, ja voi olla kymmeniä metrejä. Ctd-datassa mittaväli voi puolestaan olla hyvin lyhyt, jopa vain muutamia senttejä. Tällaisissa tapauksissa tarvitaan jotakin muuta tapaa määrittää harppauskerroksen sijainti.

Tutkin päivän, tai oikeastaan kaksikin, useita eri tapoja, ja tehokkaimmaksi ja tarkimmaksi osoittautui sovittaa dataan gam-malli, josta sitten ennustetaan metrin välein lämpötila. Ennusteiden perusteella sitten arvioidaan harppauskerroksen ylärajan sijainti. Harppauskerros ei itse asiassa välttämättä sijaitse juuri sillä kohdin, jossa lämpötila muuttuu yhden Celsius-asteen tai enemmän metrin matkalla, vaan pikemminkin siinä, missä lämpötilan nopea muutos tapahtuu, olipa absoluuttinen muutos miten nopea tahansa.

Ylläolevaa dataa käyttäen voidaan ensin piirtää havainnoista kuva:

par(mar=c(5,5,5,5))
plot(x=d$temp, y=-(d$depth), type="l", las=1, xlab="Lämpötila (Celsius-astetta)", ylab="Syvyys (m)", ylim=c(-40, -1))
abline(h=c(0:-40), col="grey75")
axis(side=4, at=-40:0, labels=c(-40:0), las=1)
abline(v=c(1:30), col="grey75")
lines(x=d$temp, y=-(d$depth), lwd=2)

Sitten sovitetaan gam-malli ja arvioidaan harppauskerroksen ylärajan sijainti funktiolla:

termokliini7<-function(d, celsius=1) {
   kliini<-which(diff(as.numeric(predict(gam(temp~s(depth, k=6), data=d), newdata=data.frame(depth=0:60))), lag=2)<=-celsius)[1]+1
   kliini
}

Funktio termokliini7() ottaa datan, sovittaa siihen kuuden vapausasteen smooth splinen, ennustaa mallista metrin välein lämpötilan arvot, ja arvioi ennusteiden välisten erotusten perusteella harppauskerroksen yläreunan sijainnin. Ennusteiden erotusten voidaan ajatella arvioivan mallin sovitteen ensimmäistä derivaattaa, mutta hommassa on vielä pieni jippo: nyt ei verratakaan kahden peräkkäisen ennusteen erotusta, vaan mukana on kahden havainnon lag eli tässä verrataankin esimerkiksi 2 m ja 4 m ennusteiden erotusta toisiinsa. Lisäksi lämpötilaraja on kahden metrin matkalla -0.5 astetta eli tässä jo nyrkkisääntöä paljon pienempi ero merkitsee harppauskerroksen alkukohtaa.

Ylläolevaan kuvaan voidaan lisätä gam-mallin ennuste:

lines(x=predict(gam(temp~s(depth, k=6), data=d), xlim=c(0,60)), y=0:-25, col="red", lwd=2)

Funktio termokliini7() antaa seuraavan tuloksen (13 m) harppauskerroksen ylärajan sijainnille:

termokliini7(d)
[1] 13

Yläraja voidaan merkitä kuvaan vielä sinisellä vaakasuoralla viivalla:

abline(h=-13, col="blue", lwd=2)

Lopputuloksena muodostuva kuva näyttää seuraavalta:

2 Responses to “Harppauskerroksen määrittäminen”

  • Laura Partanen sanoo:

    Hei.

    Kokeilin tätä esimerkkiä omaan dataani, mutta en saanut sitä toimimaan kaikilta osin.

    Minulla on datassani samanlainen tilanne, eli yritän löytää arvojen tippumiskohtia.
    Minulle tuli seuraavia ongelmia:
    termokliini7 kaavan kohdalla R herjasi: Error in s(M3$Kup, k = 6) : unused argument (k = 6).
    Poistin tuon jälkeen k=6 ja sen jälkeen sain termokliinin ajettua. En ymmärrä mistä s tulee kohtaan temp~s? Jätin s muuttujan koodiin, toimi sillä.

    Mutta nyt saan termokliini7(d) kohdasta tulokseksi NA. Se ei ole oikein. Olen etsinyt R:llä manuaalisesti tippumiskohdat ja ensimmäisessä datassa on vain yksi tippumiskohta. Siitä pitäisi tulla 5680. Mutta toisessa datassa, jossa tuli 2 vastaukseksi niin siinä putoamiskohtia on useampia ja siinä pitäisi tulla 886, 3755, 4279 jne. Sen takia haluaisin saada koodattua tippumiskohdat automaattisesti.

    Saisiko tämän tehtyä jotenkin R:llä?

    lines(x=predict(gam(temp~s(depth, k=6), data=d), xlim=c(0,60)), y=0:-25, col=”red”, lwd=2) ei toiminut, koska R herjasi: Error in xy.coords(x, y) : ’x’ and ’y’ lengths differ
    Miksi olet laittanut xlim=c(0,60), en ymmärrä missä esimerkin aineistossa on 0-60 arvot?

    • Jarno sanoo:

      Artikkelissa esitetty funktio sovittaa aineistoon GAM-mallin, jossa voidaan käyttää spline-funktiota määrittelemään ”siloite”. Tuo kaavassa esiintyvä funktio s() määrittelee millainen spline mallissa sovitetaan. Lisäksi s()-funktion argumentilla k voidaan määritellä käytettävien solmukohtien lukumäärä.

      Aineistossa oli alunperin 60 mittaustulosta, mistä johtuen kuvan x-akselin väliksi on asetettu periaatteessa sama pituus eli tuo xlim=c(0,60).

      Jos aikasarjassa on useampia muutospisteitä, jossa aikasarjan käyttäytyminen tai suunta muuttuu, voidaan ne todennäköisesti löytää tässä käytettyä funktiiota helpommin vaikkapa strucchange-paketin funktiolla breakpoints(). Vastaavankaltaisia määrityksiä voi tehdä myös vaikkapa paketilla changepoint. Toinen vaihtoehto jota usein näkee tällaisissa tilanteissa käytettävän, on joki anomalioiden tunnistamiseen käytetty menetelmä, kuten R:n paketti tsoutliers.


Category