Grafická prezentace a numerické modelování geochemických dat
Jednoduché geochemické přepočty a grafy
3.1
Geochemické indexy a binární diagramy
Vybrané citace: ALBARÈDE F. 1995. Introduction to the Geochemical Modeling. — Cambridge University Press, pp. 1–543. BOYNTON, W. V., 1984. Cosmochemistry of the rare earth elements: meteoritic studies. In: Henderson, P. (ed.) Rare Earth Elements Geochemistry. Amsterdam: Elsevier, pp. 63–114. ROLLINSON H.R. 1993. Using geochemical data: Evaluation, presentation, interpretation. Longman, London, pp 1–352. R je skvělý nástroj pro jednoduché geochemické přepočty. Na ukázku použijeme data pro granitoidy sázavské suity (sazava.data).
Cvičení 3.1
•
přepočtěte všechno železo na FeOT
FeO T = FeO + 0.89981 * Fe2 O3 •
napište funkci pro výpočet mg čísla: mg # = 100
•
MgO [mol. %] FeO + MgO
(3.2)
navrhněte funkci pro výpočet indexů A/CNK a A/NK: A / CNK =
Al 2 O3 [mol. %] CaO + Na 2 O + K 2 O
A / NK =
•
(3.1)
(3.3)
Al 2 O3 [mol. %] Na 2 O + K 2 O
(3.4)
vyneste graf A/CNK – A/NK, zobrazte v něm linie A/CNK=1, A/NK = 1 a A/CNK=A/NK
Tab 3.1. Molekulové hmotnosti použité ve Cvičení 3.1 Al2O3 101.96
Fe2O3 159.69
FeO 71.85
MgO 40.31
CaO 56.08
Na2O 61.98
Varování: dělení matice vektorem probíhá podél sloupců, ne podle řádků! > > > > >
WR<-read.table("sazava.data",sep="\t") WR<-as.matrix(WR) FeOt<-WR[,"FeO"]+WR[,"Fe2O3"]*0.89981 WR<-cbind(WR,FeOt) WR[,"FeOt"]
K2O 94.20
3/2 Sa-1 Sa-2 Sa-3 Sa-4 Sa-7 SaD-1 Gbs-1 Gbs-20 6.674743 7.653495 7.726595 9.620584 6.329810 8.562721 8.587388 8.625536 Gbs-2 Gbs-3 Po-1 Po-3 Po-4 Po-5 5.691377 8.730470 2.252873 2.479924 2.461928 2.832911 > MW<-c(101.96,159.69,71.85,40.31,56.08,61.98,94.20) > oxidy<-c("Al2O3","Fe2O3","FeO“,"MgO","CaO","Na2O","K2O") > names(MW)<-oxidy > mgno<-function(x){ mol<-t(x[,oxidy])/MW[oxidy] mg<-100 * mol["MgO",]/(mol["FeO",] + mol["MgO",]) return(mg) } > mgno(WR) Sa-1 Sa-2 Sa-3 Sa-4 Sa-7 SaD-1 Gbs-1 Gbs-20 51.16987 55.42955 51.92059 54.68820 48.07043 48.91939 72.90621 67.43446 Gbs-2 Gbs-3 Po-1 Po-3 Po-4 Po-5 76.87422 47.52161 37.27046 59.80976 30.42033 27.77656 ank<-function(x){ mol<-t(x[,oxidy])/MW[oxidy] ANK<-mol["Al2O3",]/(mol["Na2O",] + mol["K2O",]) return(ANK) }
3
1
2
A/NK
4
5
6
acnk<-function(x){ mol<-t(x[,oxidy])/MW[oxidy] ACNK<-mol["Al2O3",]/(mol["Na2O",] + mol["K2O",]+mol["CaO",]) return(ACNK) } > x<-cbind(acnk(WR),ank(WR)) > colnames(x)<-c("A/CNK","A/NK") > x A/CNK A/NK Sa-1 0.8355806 2.396569 Sa-2 0.7619109 2.307463 Sa-3 0.8079150 2.562820 Sa-4 0.7194036 2.750777 Sa-7 0.9618076 3.144008 SaD-1 0.8418505 2.983742 Gbs-1 0.4370531 3.415807 Gbs-20 0.5549533 2.992973 Gbs-2 0.7507199 5.655214 Gbs-3 0.8598972 2.565055 Po-1 0.9716646 2.331685 Po-3 1.0249173 1.891822 0.4 0.6 0.8 1.0 Po-4 1.0130401 1.868079 Po-5 1.1700486 2.432997
1.2
A/CNK
Obr. 3.1 > plot(x[,1],x[,2],xlab="A/CNK", ylab="A/NK",ylim=c(0.4,6), xlim=c(0.4,1.2),pch=WR[,"Symbol"])
Diagram A/CNK – A/NK pro horniny sázavské suity
> abline(0,1); abline(v=1); abline(h=1) # Obr. 3.1
3/3 3.2
Další užitečné grafické funkce
points (x, y)
funkce points() zobrazí body o souřadnicích [x, y]. Používá se pro přidávání dat na již vynesené grafy (funkce plot() vždy smaže předchozí a začne nový diagram) III
lines (x, y)
axis (side, at, labels)
I
funkce plot() může být volána mj. s parametrem axes=FALSE, který znamená, že nebudou zobrazeny osy. Ty pak lze vyvolat příkazem axis(), s parametry: side = 1 pro x, 2 pro y. at = vektor popisovaných hodnot labels = textový vektor s jejich popisy
A
B
C
D
custom X
Obr. 3.2 Uživatelské osy vytvořené funkcí axis
> plot(1,1,xlim=c(0,3),ylim=c(1,1),axes=F,xlab="iks",ylab="ypsilon") > axis(1,0:3,c("A","B","C","D")) > axis(2,-1:1,c("I","II","III"))
3.3
II
custom Y
vykreslí přímky; x, y jsou vektory příslušných souřadnic
Spiderdiagramy
Běžně používanými grafy jsou spiderdiagramy, tedy diagramy koncentrací ve vzorku normalizovaných nějakým standardem. Nejčastěji se tento přístup aplikuje na prvky vzácných zemin, které se normalizují chondritickými obsahy. Tab. 3.2. Obsahy prvků vzácných zemin (REE) v chondritech (Boynton 1984) La
Ce
Pr
Nd
Pm
Sm
Eu
Gd
Tb
Dy
Ho
Er
Tm
Yb
Lu
0.31
0.808
0.122
0.6
NA
0.195
0.0735
0.259
0.0474
0.322
0.0718
0.21
0.0324
0.209
0.0322
Cvičení 3.2
• •
> > >
napište funkci pro normalizaci koncentrací REE chondritickými obsahy podle Boyntona (1984 – soubor boynton.data; oddělovač je čárka) spočtěte normalizované hodnoty REE koncentrací pro analýzy sázavské suity a vyneste pomocí funkcí plot(), axis(), points() a lines() spiderdiagramy vzorků Po-1 a Po-4
WR<-read.table("sazava.data",sep="\t") WR<-as.matrix(WR[,-1]) x<-read.table("boynton.data",sep=",")
> chondrit<-as.numeric(x) # konverze na num. vektor
5.0 0.1
0.5
1.0
REE/chondrite
10.0
50.0 100.0
3/4
La
Ce
Pr
Nd
Sm
Eu
Gd
Tb
Dy
Ho
Er
Tm
Yb
Lu
Obr. 3.3. Chondritem normalizované vzory vzácných zemin pro dva trondhjemity sázavské suity (Cvičení 3.2)
> names(chondrit)<-names(x)
normalizace<-function(x,chon){ norm<-t(x[,names(chon)])/chon return(norm) } > y<-normalizace(WR,chondrit) > plot(y[,"Po-1"],type="o",log="y", axes=F,xlab="",ylab="REE/chondrit",ylim=c(0.1,100)) > axis(1, 1:length(chondrit),labels=names(chondrit)) > axis(2,) > points(y[,"Po-4"]); lines(y[,"Po-4"]) # Obr. 3.3
3.4
Ternární diagramy
Dalším užitečným nástrojem jsou ternární diagramy (Obr. 3.4). Kardinálním problémem je transformace ternárních souřadnic [a, b, c] na binární [x, y], aby je bylo možno vynést standardními funkcemi jazyka R. Pokud předpokládáme, že trojúhelník je rovnostranný, se stranou jednotkové délky, potom jeho vrcholy (spodní levý, horní, a spodní pravý) mají [x, y] souřadnice A[0, 0], B[0.5, 3 /2] and C[1, 0]. Ternární koordináty [a, b, c] datového bodu X pak mohou být transformovány na binární [x, y] pomocí (Obr. 3.4):
x =1− a −
y=
b 2
3 b 2
Obr. 3.4 Ternární diagram – vynášení bodů [a, b, c] v binárních [x, y] souřadnicích
(3.5)
(3.6)
3/5 Cvičení 3.3
• •
navrhněte funkci na kreslení ternárních diagramů vyneste ternární diagram Rb–Sr–Ba pro sázavskou suitu
triplot<-function(alab,blab,clab){ suma<-WR[,alab]+WR[,blab]+WR[,clab] a<-WR[,alab]/suma b<-WR[,blab]/suma plot(1-a-b/2,sqrt(3)*b/2,xlab="",ylab="",axes=F,xlim=c(0,1), ylim=c(0,1)) # axes=F – nebudou se vůbec vykreslovat osy # xlab="", ylab="" – a nebudou ani popsány x1<-c(0,1,.5,0) y1<-c(0,0,sqrt(3)/2,0) lines(x1,y1) text(0,0,alab) text(1,0,clab) text(0.5,sqrt(3)/2,blab)
} triplot("Ba","Rb","Sr")
3.5
Směsi několika komponent
Násobení matic (pomocí operátoru “%*%”) se nechá využít pro výpočet složení směsi, pokud známe složení jejích složek a jejich vzájemné proporce. Takto lze získat třeba složení celkové horniny, pokud známe chemismus jednotlivých minerálů a jejich modální zastoupení. Cvičení 3.4
Soubor gabro.data obsahuje složení tří minerálů olivínického gabra (Albarède 1995 – str. 7):
Tab. 3.3. Složení hlavních minerálních fází olivinického gabra SiO2 Olivín
40.01
Diopsid
54.69
Plagioklas
48.07
• > > >
Al2O3 0 0 33.37
MgO
14.35
45.64
3.27
16.51
0
0
CaO 0
Na2O 0
25.52
0
16.31
2.25
spočtěte složení gabra, obsahuje-li 40 % olivínu, 30 % diopsidu a 30 % plagioklasu
mins<-read.table("gabro.data",sep="\t") mins<-as.matrix(mins) mins
SiO2 Al2O3 FeO MgO CaO Na2O ol 40.01 0.00 14.35 45.64 0.00 0.00 di 54.69 0.00 3.27 16.51 25.52 0.00 plg 48.07 33.37 0.00 0.00 16.31 2.25
> >
FeO
f<-c(0.4,0.3,0.3) f%*%mins SiO2 Al2O3 FeO MgO CaO Na2O [1,] 46.832 10.011 6.721 23.209 12.549 0.675
3/6 3.6
Statistická analýza velkých a složitých datových souborů
Statická analýza velkých a složitých souborů geochemických dat, které například zahrnují složení vzorků hornin několika intruzí, může být velmi únavné a časově náročné. Pokud není třeba rozlišovat mezi složením jednotlivých skupin analýz (např. intruzí), běžné je použití funkce apply. Mnohem složitější je však případ, pokud potřebujeme počítat s každou skupinou zvlášť. Naštěstí faktory, v kombinaci s funkcí tapply, představují elegantní řešení takového problému. Navíc, v jazyce R existují ještě další efektivní nástroje: aggregate (x, by, FUN, ...)
Aplikuje funkci FUN na každou proměnnou (= sloupec) numerické matice nebo data framu x, přičemž je respektováno rozdělení do skupin na základě faktoru (nebo listu obsahujícího několik faktorů) by. Každá proměnná v x je tak rozdělena na podmnožiny vzorků (řádků) s určitými kombinacemi faktoru/ů by. Následně je spočtena hodnota funkce FUN pro každou takovou podmnožinu, přičemž lze specifikovat další parametry funkce (...). by (x, INDICES, FUN, ...)
Data frame x je rozdělen na podmnožiny řádků na základě argumentu INDICES (což je faktor nebo jejich list). Funkce FUN je aplikována na každou takovou podmnožinu (v principu také data frame), přičemž lze specifikovat další parametry funkce (...). Cvičení 3.5
• • • •
> > > >
Spočtěte průměry pro všechny sloupce (proměnné) souboru sazava.data Pomocí funkce summary určete základní statistické parametry charakterizující distribuci SiO2 v celém souboru Jaké jsou průměrné obsahy oxidů hlavních prvků a vybraných stopových prvků (Ba, Rb, Sr, Zr) v jednotlivých intruzích? Za použití funkce by spočtěte základní statistické parametry distribuce oxidů hlavních prvků, a to v každé intruzí zvlášť.
x<-read.table("sazava.data",sep="\t") WR<-as.matrix(x[,-1]) result<-apply(WR[,-1],2,mean,na.rm=T) round(result,2) SiO2 57.95 Na2O 2.80 Ba 883.25 Cr 67.33 Tb 0.44
> >
TiO2 0.64 K2O 1.66 Rb 51.50 La 20.05 Dy 2.49
Al2O3 16.94 P2O5 0.15 Sr 443.00 Ce 42.69 Ho 0.46
FeO 4.73 CO2 0.16 Zr 94.67 Pr 4.39 Er 1.33
labels<-as.matrix(x[,1]) labels<-factor(labels)
Fe2O3 1.75 F 0.08 Nb 6.67 Nd 16.01 Tm 0.20
MnO MgO 0.14 3.57 S H2O.PLUS 0.09 1.12 Ni Co 11.17 18.80 Sm Eu 3.29 1.40 Yb Lu 1.38 0.21
CaO 8.16 H2O. 0.06 Zn 61.08 Gd 3.08
3/7 >
tapply(WR[,"SiO2"],labels,summary)
$basic Min. 1st Qu. Median 48.84 49.63 51.72
Mean 3rd Qu. 51.78 52.90
Max. 55.80
$Pozary Min. 1st Qu. Median 62.95 66.96 69.69
Mean 3rd Qu. 68.44 71.17
Max. 71.42
$Sazava Min. 1st Qu. Median 50.72 55.09 55.17
Mean 3rd Qu. 55.74 57.73
Max. 59.98
>
aggregate(WR[,2:11],list(Rock=labels),mean)
Rock SiO2 TiO2 1 basic 51.778 0.784 2 Pozary 68.440 0.290 3 Sazava 55.738 0.774
> >
Al2O3 FeO Fe2O3 MnO MgO CaO Na2O K2O 16.872 5.664 2.640 0.174 5.644 11.1200 2.250 1.236 16.360 2.075 0.480 0.050 0.840 4.4475 3.450 1.900 17.480 5.922 1.866 0.172 3.680 8.1700 2.816 1.902
trace<-c("Rb","Sr","Ba","Zr") aggregate(WR[,trace],list(Rock=labels),mean,na.rm=T)
Rock Rb Sr Ba Zr 1 basic 34.5 346.25 676.25 65.75 2 Pozary 59.5 460.75 1291.25 157.25 3 Sazava 60.5 522.00 682.25 61.00
>
by(WR[,2:11],list(Rock=labels),summary)
Rock: basic SiO2 Min. :48.84 1st Qu.:49.63 Median :51.72 Mean :51.78 3rd Qu.:52.90 Max. :55.80 MnO Min. :0.130
TiO2 Min. :0.340 1st Qu.:0.670 Median :0.760 Mean :0.784 3rd Qu.:0.800 Max. :1.350 MgO Min. :3.160
Al2O3 Min. :13.34 1st Qu.:14.17 Median :16.98 Mean :16.87 3rd Qu.:18.23 Max. :21.64 CaO Min. : 7.22
FeO Min. :2.740 1st Qu.:5.690 Median :6.220 Mean :5.664 3rd Qu.:6.430 Max. :7.240 Na2O Min. :1.67
Fe2O3 Min. :1.47 1st Qu.:2.44 Median :2.79 Mean :2.64 3rd Qu.:3.22 Max. :3.28 K2O … etc
