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:


Vastaa

Sähköpostiosoitettasi ei julkaista. Pakolliset kentät on merkitty *