-
Notifications
You must be signed in to change notification settings - Fork 20
/
plots.R
2264 lines (2037 loc) · 93.6 KB
/
plots.R
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
598
599
600
601
602
603
604
605
606
607
608
609
610
611
612
613
614
615
616
617
618
619
620
621
622
623
624
625
626
627
628
629
630
631
632
633
634
635
636
637
638
639
640
641
642
643
644
645
646
647
648
649
650
651
652
653
654
655
656
657
658
659
660
661
662
663
664
665
666
667
668
669
670
671
672
673
674
675
676
677
678
679
680
681
682
683
684
685
686
687
688
689
690
691
692
693
694
695
696
697
698
699
700
701
702
703
704
705
706
707
708
709
710
711
712
713
714
715
716
717
718
719
720
721
722
723
724
725
726
727
728
729
730
731
732
733
734
735
736
737
738
739
740
741
742
743
744
745
746
747
748
749
750
751
752
753
754
755
756
757
758
759
760
761
762
763
764
765
766
767
768
769
770
771
772
773
774
775
776
777
778
779
780
781
782
783
784
785
786
787
788
789
790
791
792
793
794
795
796
797
798
799
800
801
802
803
804
805
806
807
808
809
810
811
812
813
814
815
816
817
818
819
820
821
822
823
824
825
826
827
828
829
830
831
832
833
834
835
836
837
838
839
840
841
842
843
844
845
846
847
848
849
850
851
852
853
854
855
856
857
858
859
860
861
862
863
864
865
866
867
868
869
870
871
872
873
874
875
876
877
878
879
880
881
882
883
884
885
886
887
888
889
890
891
892
893
894
895
896
897
898
899
900
901
902
903
904
905
906
907
908
909
910
911
912
913
914
915
916
917
918
919
920
921
922
923
924
925
926
927
928
929
930
931
932
933
934
935
936
937
938
939
940
941
942
943
944
945
946
947
948
949
950
951
952
953
954
955
956
957
958
959
960
961
962
963
964
965
966
967
968
969
970
971
972
973
974
975
976
977
978
979
980
981
982
983
984
985
986
987
988
989
990
991
992
993
994
995
996
997
998
999
1000
#' @title Plot calibrated dates
#'
#' @description Plot calibrated radiocarbon dates.
#' @param x \code{CalDates} class object containing calibrated radiocarbon dates.
#' @param ind Number indicating the index value of the calibrated radiocarbon date to be displayed. Default is 1.
#' @param label (optional) Character vector to be shown on the top-right corner of the display window.
#' @param calendar Either \code{'BP'} or \code{'BCAD'}. Indicate whether the calibrated date should be displayed in BP or BC/AD. Default is \code{'BP'}.
#' @param type Either \code{'standard'} or \code{'auc'}. If set to \code{'auc'}, displays both the normalised (dashed line) and unnormalised curves. Default is \code{'standard'}.
#' @param xlab (optional) Label for the x axis. If unspecified the default setting will be applied ("Year BP" or "Year BC/AD").
#' @param ylab (optional) Label for the y axis. If unspecified the default setting will be applied ("Radiocarbon Age").
#' @param axis4 Logical value indicating whether an axis of probabilities values should be displayed. Default is TRUE.
#' @param HPD Logical value indicating whether intervals of higher posterior density should be displayed. Default is FALSE.
#' @param credMass A numerical value indicating the size of the higher posterior density interval. Default is 0.95.
#' @param customCalCurve A three column data.frame or matrix that allows you to pass and plot a custom calibration curve if you used one during calibration. You can currently only provide one such custom curve which is used for all dates.
#' @param col The primary fill color for the calibrated date distribution.
#' @param col2 The secondary colour fill color for the calibrated date distribution, used for regions outside the higher posterior interval. Ignored when \code{HPD=FALSE}.
#' @param cex.axis The magnification to be used for axis annotation relative to the current setting of cex. Default is adjusted to 0.75.
#' @param cex.xylab The magnification to be used for x and y labels relative to the current setting of cex. Default is adjusted to 0.75.
#' @param cex.label The magnification to be used for the label on the top right corner defined by the argument \code{label}. Default is adjusted to 0.75.
#' @param add if set to \code{TRUE} the calibrated date is displayed over the existing plot. Default is \code{FALSE}.
#' @param ... Additional arguments affecting the plot.
#'
#' @seealso \code{\link{calibrate}}
#'
#' @examples
#' x <- calibrate(x=c(3402,3490,4042),errors=c(20,20,30))
#' plot(x) #display the first date
#' plot(x,2) #displays the second date
#' plot(x,3, calendar="BCAD", HPD=TRUE) #display in BC/AD with higher posterior density interval
#' @method plot CalDates
#' @export
plot.CalDates <- function(x, ind=1, label=NA, calendar="BP", type="standard", xlab=NA, ylab=NA, axis4=TRUE, HPD=FALSE, credMass=0.95, customCalCurve=NA,add=FALSE,col='grey50',col2='grey82',cex.axis=0.75,cex.xylab=0.75,cex.label=0.75,...){
types <- c("standard", "simple", "auc")
if (!type %in% types){
stop("The plot type you have chosen is not currently an option.")
}
caltimeRange =c(55000,0)
if (any(x$metadata$CalCurve %in% c("intcal13","shcal13","marine13","intcal13nhpine16","shcal13shkauri16")))
{
caltimeRange =c(50000,0)
}
if (length(x$calmatrix)>1){
grd <- data.frame(calBP=as.numeric(row.names(x$calmatrix)),PrDens=x$calmatrix[,ind])
grd <- grd[grd$PrDens >0,]
yearsBP <- grd$calBP
prob <- grd$PrDens
} else {
yearsBP <- x$grids[[ind]]$calBP
prob <- x$grids[[ind]]$PrDens
}
cra <- x$metadata$CRA[ind]
error <- x$metadata$Error[ind]
calendars <- c("BP","BCAD")
if (!calendar %in% calendars){
stop("The calendar you have chosen is not currently an option.")
}
if (yearsBP[1] > yearsBP[length(yearsBP)]){
yearsBP <- rev(yearsBP)
prob <- rev(prob)
}
yvals <- c(0,prob,0,0)
if (calendar=="BP"){
plotyears <- yearsBP
xvals <- c(plotyears[1],plotyears,plotyears[length(plotyears)], plotyears[1])
if (is.na(xlab)){ xlabel <- "Years cal BP" } else { xlabel <- xlab }
} else if (calendar=="BCAD"){
plotyears <- BPtoBCAD(yearsBP)
xvals <- c(plotyears[1],plotyears,plotyears[length(plotyears)], plotyears[1])
if (is.na(xlab)){
xlabel <- "Years BC/AD"
if (all(range(plotyears)<0)) {xlabel <- "Years BC"}
if (all(range(plotyears)>0)) {xlabel <- "Years AD"}
} else { xlabel <- xlab }
} else {
stop("Unknown calendar type")
}
xrng <- c(min(xvals[yvals>0.000001])-50,max(xvals[yvals>0.000001])+50)
if (calendar=="BP"){ xlim <- rev(xrng) } else { xlim <- xrng }
xticks <- 100*(xrng%/%100 + as.logical(xrng%%100))
if (calendar=="BP"){
xticks <- seq(xticks[1], xticks[2]-100, 100)
} else {
xticks <- seq(xticks[1]-100, xticks[2], 100)
}
yrng <- c(min(yvals[yvals>0]),max(yvals[yvals>0])+(max(yvals[yvals>0])*2))
if (!add)
{
par(cex.lab=0.75)
plot(xvals,yvals, type="n", xlab=xlabel, ylab="", ylim=yrng, xlim=xlim, xaxt='n', yaxt='n', cex.lab=cex.xylab,cex.axis=cex.axis,...)
}
xticksLab <- xticks
if (calendar=="BCAD")
{
if (any(xticksLab==0)){xticksLab[which(xticksLab==0)]=1}
xticks[which(xticks>1)]=xticks[which(xticks>1)]-1
}
if(!add) {axis(1, at=xticks, labels=abs(xticksLab), las=2, cex.axis=cex.axis)}
if (!add&axis4){ axis(4, cex.axis=cex.axis) }
if (!HPD){
polygon(xvals,yvals, col=col, border=col)
} else {
polygon(xvals,yvals, col=col2, border=col2)
hdres <- hpdi(x,credMass=credMass)[[ind]]
if(calendar=="BCAD"){hdres=1950-hdres}
for (i in 1:nrow(hdres))
{
index <- which(xvals%in%hdres[i,1]:hdres[i,2])
polygon(c(xvals[index],xvals[index[length(index)]],xvals[index[1]]),c(yvals[index],0,0), col=col, border=col)
}
}
if (type=="standard" | type=="auc"){
if (type=="auc"){
lines(xvals, yvals/sum(yvals), col="black", lty="dotted")
}
if(!add)
{
par(new=TRUE)
cradf1 <- data.frame(CRA=caltimeRange[1]:0,Prob=dnorm(caltimeRange[1]:0, mean=cra, sd=error))
cradf1 <- cradf1[cradf1$Prob>0.0001,]
ylim <- c(cra-(12*error),cra+(8*error))
cradf1$RX <- reScale(cradf1$Prob, to=c(xlim[1],(xlim[1]+diff(xlim)*0.33)))
yticks <- ylim[1]:ylim[2]
yticks <- yticks[yticks %% 200 == 0]
plot(cradf1$RX,cradf1$CRA,type="l", axes=FALSE, xlab=NA, ylab=NA, xlim=xlim, ylim=ylim, col=rgb(144,238,144,120,maxColorValue=255))
polygon(c(cradf1$RX,rev(cradf1$RX)),c(cradf1$CRA,rep(xlim[1],length(cradf1$CRA))), col=rgb(144,238,144,80,maxColorValue=255), border=NA)
axis(side=2, at=yticks, labels=abs(yticks),las=2, cex.axis=cex.axis)
if (is.na(ylab)){
title(line=3, ylab="Radiocarbon Age", cex.lab=cex.xylab)
} else {
title(line=3, ylab=ylab, cex=cex.xylab)
}
}
calcurvemetadata <- x$metadata$CalCurve[ind]
calcurvecheck <- TRUE
if (calcurvemetadata == "custom" & !any(class(customCalCurve) %in% c("data.frame","matrix","array"))){
calcurvecheck <- FALSE
}
if (calcurvecheck&!add){
if (calcurvemetadata == "custom"){
cc <- as.data.frame(customCalCurve)[,1:3]
names(cc) <- c("BP","CRA","Error")
} else {
calCurveFile <- paste(system.file("extdata", package="rcarbon"), "/", calcurvemetadata,".14c", sep="")
options(warn=-1)
cc <- readLines(calCurveFile, encoding="UTF-8")
cc <- cc[!grepl("[#]",cc)]
cc.con <- textConnection(cc)
cc <- read.csv(cc.con, header=FALSE, stringsAsFactors=FALSE)
close(cc.con)
options(warn=0)
names(cc) <- c("BP","CRA","Error","D14C","Sigma")
}
if (calendar=="BCAD"){
tmp <- (xrng-1950)*-1
cc$RX <- 1950-cc$BP
} else {
tmp <- xrng
cc$RX <- cc$BP
}
cc <- cc[cc$BP >= min(tmp) & cc$BP < max(tmp),]
cc$Hi <- cc$CRA + 10
cc$Lo <- cc$CRA - 10
ccbox <- c((max(yrng)*0.2),(max(yrng)*0.9))
polygon(c(cc$RX,rev(cc$RX)),c(cc$Hi,rev(cc$Lo)), col=rgb(255,140,0,120,maxColorValue=255), border=NA)
lines(cc$RX,cc$CRA, col=rgb(255,140,0,60,maxColorValue=255))
}
}
if (!is.na(label)){
legend("topright", label, bty="n", cex=cex.label)
}
}
#' @title Plot multiple dates
#'
#' @description Plot multiple radiocarbon dates.
#' @param x A CalDates class object with length > 1.
#' @param type Whether the calibrated dates are displayed as distributions (\code{'d'}) or as horizontal bars (\code{'b'}) spanning the HPD interval. Default is \code{'d'}.
#' @param calendar Either \code{'BP'} or \code{'BCAD'}. Indicate whether the calibrated date should be displayed in BP or BC/AD. Default is \code{'BP'}.
#' @param HPD Logical value indicating whether intervals of higher posterior density should be displayed. Default is FALSE.
#' @param credMass A numerical value indicating the size of the higher posterior density interval. Default is 0.95.
#' @param decreasing Whether dates should be plotted with decreasing order of median calibrated date (i.e. old to new; TRUE) or increasing order (i.e. new to old; FALSE). If set to NULL the dates plotted in the supplied order. Default is NULL
#' @param label Whether the ID of each date should be displayed. Default is TRUE.
#' @param label.pos Relative position of the label in relation to the calibrated distribution expressed in quantiles. Default is 0.5 (median).
#' @param label.offset Horrizontal offset of label position in number of years. Default is 0.
#' @param xlim the x limits of the plot. In BP or in BC/AD depending on the choice of the parameter \code{calender}. Notice that if BC/AD is selected BC ages should have a minus sign (e.g. \code{c(-5000,200)} for 5000 BC to 200 AD).
#' @param xlab (optional) Label for the x axis. If unspecified the default setting will be applied ("Year BP" or "Year BC/AD").
#' @param ylab (optional) Label for the y axis.
#' @param col.fill A vector of primary fill color for the calibrated date distribution. Default is 'grey50'.
#' @param col.fill2 A vector of secondary secondary colour fill color for the calibrated date distribution, used for regions outside the higher posterior interval. Ignored when \code{HPD=FALSE}. Default is 'grey82'.
#' @param col.line A vector of line color (ignored when \code{type} is set to \code{'d'}. Default is 1.
#' @param lwd Line width (ignored when \code{type} is set to \code{'d'}). Default is 1.
#' @param cex.lab The magnification to be used for x and y labels relative to the current setting of cex. Default is adjusted to 1.
#' @param cex.id The magnification to be used the date labels relative to the current setting of cex. Default is adjusted to 1.
#' @param cex.axis The magnification to be used for axis annotation relative to the current setting of cex. Default is adjusted to 1.
#' @param ydisp Whether the y axis should be displayed. Ignored when \code{type} is set to \code{'b'}. Default is FALSE
#' @param gapFactor Defines spacing between calibrated dates (when \code{type} is set to \code{'d'}) or the distance between the lines and the labels (when \code{type} is set to \code{'b'}) as proportion of individual y-axis ranges. Default is 0.2.
#' @param rescale Whether probability distributions should be rescaled (applicable only when \code{type} is set to \code{'d'}, default is FALSE).
#' @seealso \code{\link{calibrate}}
#'
#' @examples
#' data("emedyd")
#' tellaswad = subset(emedyd,SiteName=='Tell Aswad')
#' x = calibrate(tellaswad$CRA,tellaswad$Error,ids=tellaswad$LabID)
#' multiplot(x,HPD=TRUE,decreasing=TRUE,label=FALSE,gapFactor = 0.1)
#' multiplot(x,type='b',calendar='BCAD',cex.id = 0.5,lwd=2,gapFactor = 0.5)
#' @import stats
#' @import grDevices
#' @import graphics
#' @import utils
#' @export
multiplot<- function(x,type='d',calendar='BP',HPD=FALSE,credMass=0.95,decreasing=NULL,label=TRUE,label.pos=0.5,label.offset=0,xlim=NULL,xlab=NA,ylab=NA,col.fill='grey50',col.fill2='grey82',col.line='black',lwd=1,cex.id=1,cex.lab=1,cex.axis=1,ydisp=FALSE,gapFactor=0.2,rescale=FALSE)
{
if(length(lwd)==1){lwd=rep(lwd,length(x))}
if(length(col.line)==1){col.line=rep(col.line,length(x))}
if(length(col.fill)==1){col.fill=rep(col.fill,length(x))}
if(length(col.fill2)==1){col.fill2=rep(col.fill2,length(x))}
if (!is.null(decreasing))
{
ord = order(medCal(x),decreasing=decreasing)
x = x[ord]
col.line = col.line[ord]
col.fill = col.fill[ord]
col.fill2 = col.fill2[ord]
}
date.pos = qCal(x,p=label.pos) + label.offset
calendars <- c("BP","BCAD")
if (!calendar %in% calendars){
stop("The calendar you have chosen is not currently an option.")
}
# Estimate general xlim
if (is.null(xlim))
{
if(anyNA(x$grids))
{
tmp = apply(x$calmatrix,1,sum)
st = as.numeric(names(tmp[which(tmp>0)[1]-1]))
en = as.numeric(names(tmp[which(tmp>0)[length(which(tmp>0))]+1]))
} else {
st = max(unlist(lapply(x$grids,function(x){max(x$calBP)})))
en = min(unlist(lapply(x$grids,function(x){min(x$calBP)})))
}
edge = 0.1*abs(st-en)
xlim = c(st+edge,en-edge)
}
yearsBP = xlim[1]:xlim[2]
if (calendar=="BP"){
plotyears <- yearsBP
xvals <- c(plotyears[1],plotyears,plotyears[length(plotyears)], plotyears[1])
if (is.na(xlab)){ xlabel <- "Years cal BP" } else { xlabel <- xlab }
} else if (calendar=="BCAD"){
plotyears <- BPtoBCAD(yearsBP)
xlim <- BPtoBCAD(xlim)
xvals <- c(plotyears[1],plotyears,plotyears[length(plotyears)], plotyears[1])
if (is.na(xlab)){
xlabel <- "Years BC/AD"
if (all(range(plotyears)<0)) {xlabel <- "Years BC"}
if (all(range(plotyears)>0)) {xlabel <- "Years AD"}
} else { xlabel <- xlab }
} else {
stop("Unknown calendar type")
}
# Plot
if (type=='b')
{
bse = hpdi(x,credMass=credMass)
plot(0,0,xlim=xlim,ylim=c(0,length(bse)+1),axes=F,xlab=xlabel,ylab="",type='n')
for (i in 1:length(bse))
{
tmp = matrix(bse[[i]][,-3],ncol=2)
if (calendar=='BP')
{
apply(tmp,1,function(x,y,lwd,col){lines(c(x),c(y,y),lwd=lwd,col=col)},y=i,lwd=lwd[i],col=col.line[i])
if(label){text(x=date.pos[i],y=i+gapFactor,label=x$metadata$DateID[i],cex=cex.id)}
}
if (calendar=='BCAD')
{
apply(tmp,1,function(x,y,lwd,col){lines(BPtoBCAD(c(x)),c(y,y),lwd=lwd,col=col)},y=i,lwd=lwd[i],col=col.line[i])
if(label){text(x=BPtoBCAD(date.pos[i]),y=i+gapFactor,label=x$metadata$DateID[i],cex=cex.id)}
}
}
}
if (type=='d')
{
if(anyNA(x$grids))
{
if (rescale)
{
x$calmatrix=apply(x$calmatrix,2,reScale)
}
tmp = apply(x$calmatrix,1,max)
ylim = as.numeric(c(0,tmp[which.max(tmp)]))
} else {
if (rescale)
{
for (i in 1:length(x))
{
x$grids[[i]]$PrDens = reScale(x$grids[[i]]$PrDens)
}
}
ylim = c(0,max(unlist(lapply(x$grids,function(x){max(x$PrDens)}))))
}
# max ylim: combine ylim giving as a space 1/7 of the original distance
gap = abs(diff(ylim))*gapFactor
# generate ylim sequences:
YLIMs = vector("list",length=length(x))
YLIMs[[1]] = ylim
for (i in 2:length(x))
{
YLIMs[[i]]=c(YLIMs[[i-1]][2]+gap,YLIMs[[i-1]][2]+gap+abs(diff(ylim)))
}
plot(0, 0, xlim=xlim, ylim=c(min(unlist(YLIMs)),max(unlist(YLIMs))+gap), type="n", ylab=ylab, xlab=xlabel, axes=F,cex.lab=cex.lab)
for (i in 1:length(x))
{
tmpYlim= YLIMs[[i]]
if (ydisp)
{axis(2,at=c(tmpYlim[1],median(tmpYlim),max(tmpYlim)),labels=round(c(min(ylim),median(ylim),max(ylim)),2),las=2,cex.axis=cex.axis)}
if(anyNA(x$grid))
{
years = as.numeric(rownames(x$calmatrix))
PrDens = x$calmatrix[,i]
ii = which(PrDens>0)[1]-1
jj = max(which(PrDens>0))+1
years = years[ii:jj]
PrDens = PrDens[ii:jj]
} else {
years=x$grid[[i]]$calBP
PrDens = x$grid[[i]]$PrDens
}
if (calendar=='BCAD'){years=BPtoBCAD(years)}
xvals = c(years,rev(years))
yvals = c(PrDens+tmpYlim[1],rep(0,length(years))+tmpYlim[1])
if (!HPD){
polygon(xvals,yvals, col=col.fill[i], border=col.fill[i])
} else {
polygon(xvals,yvals, col=col.fill2[i], border=col.fill2[i])
hdres <- hpdi(x,credMass=credMass)[[i]]
for (j in 1:nrow(hdres))
{
if (calendar=='BCAD')
{
index <- which(xvals%in%BPtoBCAD(hdres[j,1]:hdres[j,2]))
} else {
index <- which(xvals%in%hdres[j,1]:hdres[j,2])
}
polygon(c(xvals[index],xvals[index[length(index)]],xvals[index[1]]),c(yvals[index],min(tmpYlim),min(tmpYlim)), col=col.fill[i], border=col.fill[i])
}
}
if (label)
{
xx = date.pos[i]
if (calendar=='BCAD'){xx = BPtoBCAD(xx)}
ylabel=ifelse(rescale,max(yvals)-0.5,max(yvals)+gap/2)
text(x=xx,y=ylabel,labels=x$metadata$DateID[i],cex=cex.id)
}
}
}
# draw x-axis
if (calendar=="BP"){
rr <- range(pretty(plotyears))
axis(side=1,at=seq(rr[2],rr[1],-100),labels=NA,tck = -.01,cex.axis=cex.axis)
axis(side=1,at=pretty(plotyears),labels=abs(pretty(plotyears)),cex.axis=cex.axis)
} else if (calendar=="BCAD"){
yy <- plotyears
rr <- range(pretty(yy))
prettyTicks <- seq(rr[1],rr[2],100)
prettyTicks[which(prettyTicks>=0)] <- prettyTicks[which(prettyTicks>=0)]-1
axis(side=1,at=prettyTicks, labels=NA,tck = -.01,cex.axis=cex.axis)
py <- pretty(yy)
pyShown <- py
if (any(pyShown==0)){pyShown[which(pyShown==0)]=1}
py[which(py>1)] <- py[which(py>1)]-1
axis(side=1,at=py,labels=abs(pyShown),cex.axis=cex.axis)
}
}
#' @title Plot result of Monte-Carlo simulation of observed versus modelled SPDs
#'
#' @description The function visualises the observed summed probability distribution of radiocarbon dates along with a simulation envelope for the null model and regions of positive and negative deviation.
#'
#' @param x A \code{SpdModelTest} class object generated using the \code{\link{modelTest}} function.
#' @param calendar Either \code{'BP'} or \code{'BCAD'}. Indicate whether the calibrated date should be displayed in BP or BC/AD. Default is \code{'BP'}.
#' @param type Whether to display SPDs ('spd') or rates of change ('roc'). Default is 'spd'.
#' @param xlim the x limits of the plot. In BP or in BC/AD depending on the choice of the parameter \code{calender}. Notice that if BC/AD is selected BC ages should have a minus sign (e.g. \code{c(-5000,200)} for 5000 BC to 200 AD).
#' @param ylim the y limits of the plot.
#' @param col.obs Line colour for the observed SPD.
#' @param lwd.obs Line width for the observed SPD.
#' @param xaxs The style of x-axis interval calculation (see \code{\link{par}})
#' @param yaxs The style of y-axis interval calculation (see \code{\link{par}})
#' @param bbty Display options; one between \code{'b'},\code{'n'},and \code{'f'}. See details below.
#' @param bbtyRet Whether details on the intervals of positive and negative deviations are returned. Default is TRUE.
#' @param drawaxes A logical value determining whether the axes should be displayed or not. Default is TRUE.
#' @param ... Additional arguments affecting the plot
#' @details The argument \code{bbty} controls the display options of the Monte-Carlo Test. Default settings (\code{bbty='f'}) displays the observed SPD (solid black line), the simulation envelope of the fitted model (shaded grey polygon) and regions of significance positive (red semi-transparent rectangle) and negative (blue semi-transparent rectangle) deviation. The option \code{bbty='b'} removes the regions of positive/negative deviations, whilst the option \code{bbty='n'} displays the simulation envelope on existing plot.
#' @seealso \code{\link{modelTest}}
#' @import stats
#' @import grDevices
#' @import graphics
#' @import utils
#' @method plot SpdModelTest
#' @export
plot.SpdModelTest <- function(x, calendar="BP", type='spd', ylim=NA, xlim=NA, col.obs="black", lwd.obs=0.5, xaxs="i", yaxs="i", bbty="f", bbtyRet=TRUE, drawaxes=TRUE, ...){
obs <- x$result[,1:2]
envelope <- x$result[,3:4]
if (!type%in%c('spd','roc'))
{
stop("The argument 'type' should be either 'spd' or 'roc'.")
}
naindex = numeric()
if (type=='roc')
{
naindex = which(is.na(x$result.roc$roc))
obs <- x$result.roc[,1:2]
envelope <- x$result.roc[,3:4]
colnames(obs)<-c("calBP","PrDens")
colnames(envelope)<-c("lo","hi")
}
if (calendar=="BP"){
obs$Years <- obs$calBP
xlabel <- "Years cal BP"
if (any(is.na(xlim))){ xlim <- c(max(obs$Years),min(obs$Years)) }
} else if (calendar=="BCAD"){
xlabel <- 'Years BC/AD'
obs$Years <- BPtoBCAD(obs$calBP)
if (all(range(obs$Years)<0)){xlabel <- "Years BC"}
if (all(range(obs$Years)>0)){xlabel <- "Years AD"}
if (any(is.na(xlim))){xlim <- c(min(obs$Years),max(obs$Years)) }
} else {
stop("Unknown calendar type")
}
if (any(is.na(ylim)))
{
ylim <- c(0, max(envelope[,"hi"], obs$PrDens,na.rm=TRUE)*1.1)
if (type=='roc') {ylim[1] <- min(c(envelope[,"lo"],obs$PrDens),na.rm=TRUE)}
}
booms <- which(obs$PrDens>envelope[,2])
busts <- which(obs$PrDens<envelope[,1])
baseline <- rep(NA,nrow(obs))
ylab = 'Summed Probability'
if (type=='roc')
{
ylab = 'Rate of Change'
colpts = rep(col.obs,length(obs$PrDens))
colpts[booms] = 'red'
colpts[busts] = 'blue'
}
if (drawaxes & bbty != "n"){
plot(obs$Years, obs$PrDens, xlim=xlim, ylim=ylim, xlab=xlabel, ylab=ylab, type="l", col=col.obs, lwd=lwd.obs, xaxs=xaxs, yaxs=yaxs, axes=FALSE,...)
if(type=='roc')
{
abline(h=0,lty=2,lwd=1)
}
} else if (bbty != "n")
{
plot(obs$Years, obs$PrDens, xlim=xlim, ylim=ylim, xlab="", ylab="", type="l", col=col.obs, lwd=lwd.obs, xaxs=xaxs, yaxs=yaxs, axes=FALSE, ...)
if (type=='roc')
{
abline(h=0,lty=2,lwd=1)
}
}
if (drawaxes)
{
box()
axis(side=2)
}
boomPlot <- baseline
if (length(booms)>0){ boomPlot[booms]=obs[booms,2] }
bustPlot <- baseline
if (length(busts)>0){ bustPlot[busts]=obs[busts,2] }
boomBlocks <- vector("list")
counter <- 0
state <- "off"
for (i in 1:length(boomPlot)){
if (!is.na(boomPlot[i])&state=="off"){
counter <- counter+1
boomBlocks <- c(boomBlocks,vector("list",1))
boomBlocks[[counter]] <- vector("list",2)
boomBlocks[[counter]][[1]] <- boomPlot[i]
boomBlocks[[counter]][[2]] <- obs[i,"Years"]
state <- "on"
}
if (state=="on"){
if (!is.na(boomPlot[i])){
boomBlocks[[counter]][[1]] <- c(boomBlocks[[counter]][[1]],boomPlot[i])
boomBlocks[[counter]][[2]] <- c(boomBlocks[[counter]][[2]],obs[i,"Years"])
}
if (is.na(boomPlot[i])){
state <- "off"
}
}
}
bustBlocks <- vector("list")
counter <- 0
state <- "off"
for (i in 1:length(bustPlot)){
if (!is.na(bustPlot[i])&state=="off"){
counter <- counter+1
bustBlocks <- c(bustBlocks,vector("list",1))
bustBlocks[[counter]] <- vector("list",2)
bustBlocks[[counter]][[1]] <- bustPlot[i]
bustBlocks[[counter]][[2]] <- obs[i,"Years"]
state <- "on"
}
if (state=="on"){
if (!is.na(bustPlot[i])){
bustBlocks[[counter]][[1]] <- c(bustBlocks[[counter]][[1]],bustPlot[i])
bustBlocks[[counter]][[2]] <- c(bustBlocks[[counter]][[2]],obs[i,"Years"])
}
if (is.na(bustPlot[i])){
state <- "off"
}
}
}
if (length(booms)>0){
for (i in 1:length(boomBlocks)){
if (bbty=="f"){
polygon(c(boomBlocks[[i]][[2]],rev(boomBlocks[[i]][[2]])),c(rep(+100,length(boomBlocks[[i]][[1]])),rep(-100,length(boomBlocks[[i]][[1]]))),col=rgb(0.7,0,0,0.2),border=NA)
} else if (bbty %in% c("s","b","n")){
} else {
stop("Incorrect bbty argument.")
}
}
}
if (length(busts)>0){
for (i in 1:length(bustBlocks)){
if (bbty=="f"){
polygon(c(bustBlocks[[i]][[2]],rev(bustBlocks[[i]][[2]])),c(rep(+100,length(bustBlocks[[i]][[1]])),rep(-100,length(bustBlocks[[i]][[1]]))),col=rgb(0,0,0.7,0.2),border=NA)
} else if (bbty %in% c("s","b","n")){
} else {
stop("Incorrect bbty argument.")
}
}
}
if (length(naindex)>0)
{polygon(x=c(obs[-naindex,"Years"],rev(obs[-naindex,"Years"])),y=c(envelope[-naindex,1],rev(envelope[-naindex,2])),col=rgb(0,0,0,0.2),border=NA)
} else
{
polygon(x=c(obs[,"Years"],rev(obs[,"Years"])),y=c(envelope[,1],rev(envelope[,2])),col=rgb(0,0,0,0.2),border=NA)
}
if (drawaxes & bbty != "n" & calendar=="BP"){
rr <- range(pretty(obs[,"Years"]))
axis(side=1,at=seq(rr[2],rr[1],-100),labels=NA,tck = -.01)
axis(side=1,at=pretty(obs[,"Years"]),labels=abs(pretty(obs[,"Years"])))
} else if (drawaxes & bbty != "n" & calendar=="BCAD"){
yy <- obs[,"Years"]
rr <- range(pretty(yy))
prettyTicks <- seq(rr[1],rr[2],+100)
prettyTicks[which(prettyTicks>=0)] <- prettyTicks[which(prettyTicks>=0)]-1
axis(side=1,at=prettyTicks, labels=NA,tck = -.01)
py <- pretty(yy)
pyShown <- py
if (any(pyShown==0)){pyShown[which(pyShown==0)]=1}
py[which(py>1)] <- py[which(py>1)]-1
axis(side=1,at=py,labels=abs(pyShown))
}
bbp <- list(booms=boomBlocks, busts=bustBlocks)
class(bbp) <- c("BBPolygons",class(bbp))
if ((bbty %in% c("n","b")) & bbtyRet){ return(bbp) }
}
#' @title Plot the median values of calibrated radiocarbon dates or bins
#'
#' @description Plot the median values of multiple calibrated radiocarbon dates or bins in a barcode-like strip.
#' @param x A vector containing median values obtained from \code{\link{medCal}} or \code{\link{binMed}}
#' @param yrng y-axis range of the bars.
#' @param width width of the bars (optional)
#' @param col color of the bars
#' @param border the color to draw the border. Use border = NA to omit borders.
#' @param ... Additional arguments affecting the plot
#'
#'
#' @seealso \code{\link{medCal}}; \code{\link{binMed}}
#' @import stats
#' @import grDevices
#' @import graphics
#' @import utils
#' @examples
#'\dontrun{
#' #Load EUROEVOL Data
#' data(euroevol)
#'
#' #Subset Danish Dates
#' denmark <- subset(euroevol,Country=="Denmark")
#'
#' #Calibrate and Bin
#' denmarkDates <- calibrate(x=denmark$C14Age,errors=denmark$C14SD)
#' denmarkBins <- binPrep(sites=denmark$SiteID,ages=denmark$C14Age,h=200) #200 years bin size
#'
#' #Compute median date for each bin
#' bm <- binMed(x=denmarkDates,bins=denmarkBins)
#'
#' #Compute median date for each date
#' dm <- medCal(denmarkDates)
#'
#' #Compute SPD
#' denmarkSPD <- spd(x=denmarkDates,bins=denmarkBins,timeRange=c(10000,4000))
#'
#' #Plot SPD and barCodes of median dates
#' plot(denmarkSPD,runm=200)
#' barCodes(dm,yrng=c(0,0.01))
#'
#' #Plot SPD and barCodes of median bins in BC/AD
#' plot(denmarkSPD,runm=200,calendar="BCAD")
#' barCodes(BPtoBCAD(bm),yrng=c(0,0.01))
#'}
#' @import stats
#' @import grDevices
#' @import graphics
#' @import utils
#' @export
barCodes <- function(x, yrng=c(0,0.03), width=20, col=rgb(0,0,0,25,maxColorValue=255), border=NA, ...){
barcodes <- x
halfbw <- width/2
for (a in 1:length(barcodes)){
polygon(x=c(barcodes[a]-halfbw,barcodes[a]-halfbw,barcodes[a]+halfbw,barcodes[a]+halfbw,barcodes[a]-halfbw),y=c(yrng[1],yrng[2],yrng[2],yrng[1],yrng[1]), border=border, col=col, ...)
}
}
#
# crossHairs <- function(x, pch.pts=19, cex.pts=1, fixXorder=FALSE, rescaleY=FALSE,...){
#
# if (!"quickMarks" %in% class(x)){
# stop("Input must be of class \"quickMarks\"")
# }
# if (rescaleY){
# cra <- reScale(x$CRA)
# error <- x$Error / (max(x$CRA)-min(x$CRA))
# } else {
# cra <- x$CRA
# error <- x$Error
# }
# if (fixXorder){
# xstart <- x$q68s *-1
# xend <- x$q68e *-1
# xmed <- x$qMed *-1
# } else {
# xstart <- x$q68s
# xend <- x$q68e
# xmed <- x$qMed
# }
# for (a in 1:nrow(x)){
# lines(c(xstart[a],xend[a]), c(cra[a],cra[a]), ...)
# lines(c(xmed[a],xmed[a]), c(cra[a]-error[a],cra[a]+error[a]), ...)
# points(xmed[a],cra[a], pch=pch.pts, cex=cex.pts, ...)
# }
# }
#
#' @title Plot a summed probability distribution
#'
#' @description Plot a summed probability distribution (SPD) of radiocarbon dates
#' @param x A \code{CalSPD} class object.
#' @param runm A number indicating the window size of the moving average to smooth the SPD. If set to \code{NA} no moving average is applied. Default is NA
#' @param calendar Either \code{'BP'} or \code{'BCAD'}. Indicate whether the calibrated date should be displayed in BP or BC/AD. Default is \code{'BP'}.
#' @param type Either \code{'standard'} or \code{'simple'}. The former visualise the SPD as an area graph, while the latter as line chart.
#' @param xlim the x limits of the plot. In BP or in BC/AD depending on the choice of the parameter \code{calender}. Notice that if BC/AD is selected BC ages should have a minus sign (e.g. \code{c(-5000,200)} for 5000 BC to 200 AD).
#' @param ylim the y limits of the plot.
#' @param ylab (optional) Label for the y axis. If unspecified the default setting will be applied ("Summed Probability")
#' @param spdnormalised A logical variable indicating whether the total probability mass of the SPD is normalised to sum to unity.
#' @param rescale A logical variable indicating whether the SPD should be rescaled to range 0 to 1.
#' @param fill.p Fill colour for the SPD
#' @param border.p Border colour for the SPD
#' @param xaxt Whether the x-axis tick marks should be displayed (\code{xaxt='s'}, default) or not (\code{xaxt='n'}).
#' @param yaxt Whether the y-axis tick marks should be displayed (\code{xaxt='s'}, default) or not (\code{xaxt='n'}).
#' @param cex.axis The magnification to be used for axis annotation relative to the current setting of cex. Default is 1.
#' @param add Whether or not the new graphic should be added to an existing plot.
#' @param ... Additional arguments affecting the plot
#'
#'
#'
#' @seealso \code{\link{spd}}; \code{\link{plot.CalGrid}}
#' @examples
#' \dontrun{
#' data(emedyd)
#' levant <- emedyd[emedyd$Region=="1"|emedyd$Region=="2",]
#' bins <- binPrep(levant$SiteName, levant$CRA, h=50)
#' x <- calibrate(levant$CRA, levant$Error, normalised=FALSE)
#' spd.levant <- spd(x, bins=bins, timeRange=c(17000,8000))
#' spd.northernlevant <- spd(x[levant$Region=="2"], bins=bins[levant$Region=="2"],
#' timeRange=c(17000,8000))
#' plot(spd.levant, runm=50, xlim=c(16000,9000))
#' plot(spd.northernlevant, runm=50, add=TRUE, fill.p="black")
#' legend("topleft", legend=c("All Levant dates","Northern Levant only"),
#' fill=c("grey75","black"), border=NA)
#' plot(spd.levant, runm=50, xlim=c(16000,9000), type="simple")
#' plot(spd.northernlevant, runm=50, col="red", type="simple", add=TRUE)
#'}
#' @import stats
#' @import grDevices
#' @import graphics
#' @import utils
#' @method plot CalSPD
#' @export
plot.CalSPD <- function(x, runm=NA, calendar="BP", type="standard", xlim=NA, ylim=NA, ylab="Summed Probability", spdnormalised=FALSE, rescale=FALSE, fill.p="grey75", border.p=NA, xaxt='s', yaxt='s', add=FALSE, cex.axis=1, ...){
types <- c("standard","simple")
if (!type %in% types){
stop("The plot type you have chosen is not currently an option.")
}
spdvals <- x$grid$PrDens
if (!is.na(runm)){ spdvals <- runMean(spdvals, runm, edge="fill") }
if (spdnormalised){ spdvals <- spdvals/sum(spdvals) }
if (rescale){ spdvals <- reScale(spdvals) }
if (any(is.na(ylim))){ ylim <- c(0,max(spdvals)*1.1) }
if (calendar=="BP"){
plotyears <- x$grid$calBP
xlabel <- "Years cal BP"
if (any(is.na(xlim))){ xlim <- c(max(plotyears),min(plotyears)) }
} else if (calendar=="BCAD"){
plotyears <- BPtoBCAD(x$grid$calBP)
xlabel <- "Years BC/AD"
if (all(range(plotyears)<0)){xlabel <- "Years BC"}
if (all(range(plotyears)>0)){xlabel <- "Years AD"}
if (any(is.na(xlim))){ xlim <- c(min(plotyears),max(plotyears)) }
} else {
stop("Unknown calendar type")
}
if (xaxt=='n'){ xlabel <- "" }
if (yaxt=='n'){ ylabel <- "" } else { ylabel <- ylab }
if (type=="standard"){
par(xaxs="i")
par(yaxs="i")
if (!add){
plot(plotyears, spdvals, xlim=xlim, ylim=ylim, type="l", col="white", ylab=ylabel, xlab=xlabel, xaxt="n", yaxt=yaxt, ...)
}
polygon(c(plotyears,rev(plotyears)),c(spdvals,rep(0,length(spdvals))),border=border.p, col=fill.p)
} else if (type=="simple"){
if (!add){
plot(plotyears, spdvals, xlim=xlim, ylim=ylim, type="l", ylab="", xlab=xlabel, xaxt="n", yaxt=yaxt, ...)
} else {
lines(plotyears, spdvals, xlim=xlim, ylim=ylim, type="l", ylab="", xlab=xlabel, ...)
}
}
box()
if (calendar=="BP" & xaxt!="n"){
rr <- range(pretty(plotyears))
axis(side=1,at=seq(rr[2],rr[1],-100),labels=NA,tck = -.01,cex.axis=cex.axis)
axis(side=1,at=pretty(plotyears),labels=abs(pretty(plotyears)),cex.axis=cex.axis)
} else if (calendar=="BCAD" & xaxt!="n"){
yy <- plotyears
rr <- range(pretty(yy))
prettyTicks <- seq(rr[1],rr[2],+100)
prettyTicks[which(prettyTicks>=0)] <- prettyTicks[which(prettyTicks>=0)]-1
axis(side=1,at=prettyTicks, labels=NA,tck = -.01)
py <- pretty(yy)
pyShown <- py
if (any(pyShown==0)){pyShown[which(pyShown==0)]=1}
py[which(py>1)] <- py[which(py>1)]-1
axis(side=1,at=py,labels=abs(pyShown))
}
}
#' @title Plot a summed probability distribution (from a CalGrid object)
#'
#' @description Plot a summed radiocarbon probability distribution. This is a basic function for plotting SPDs that have been constructed manually or by calibrating a summed or otherwise irregular CRA grid. In most instances, it is sensible to use \code{plot.CalSPD} instead.
#'
#' @param x A "CalGrid" class object of summed probabilities per calendar year BP.
#' @param runm A number indicating the window size of the moving average to smooth the SPD. If set to \code{NA} no moving average is applied. Default is NA
#' @param calendar Either \code{'BP'} or \code{'BCAD'}. Indicate whether the calibrated date should be displayed in BP or BC/AD. Default is \code{'BP'}.
#' @param fill.p Fill colour of the polygon depicting the summed probability distribution.
#' @param border.p Border colour of the polygon depicting the summed probability distribution.
#' @param xlim the x limits of the plot. In BP or in BC/AD depending on the choice of the parameter \code{calender}. Notice that if BC/AD is selected BC ages should have a minus sign (e.g. \code{c(-5000,200)} for 5000 BC to 200 AD).
#' @param ylim Adjust y axis limits (otherwise sensible default).
#' @param cex.lab Size of label text.
#' @param cex.axis Size of axis text.
#' @param mar Adjust margins around plot.
#' @param add Whether or not the new graphic should be added to an existing plot.
#' @param ... Additional arguments affecting the plot
#'
#' @seealso \code{\link{spd}}; \code{\link{plot.CalSPD}}
#' @examples
#' data(euroevol)
#' mycaldates <- calibrate(euroevol[1:10,"C14Age"], euroevol[1:10,"C14SD"], normalised=FALSE)
#' myspd <- spd(mycaldates, timeRange=c(8000,2000))
#' plot(myspd) #ordinary plot using \code{plot.CalSPD}
#' plot(myspd$grid) #working plot using the internal CalGrid object
#' @import stats
#' @import grDevices
#' @import graphics
#' @import utils
#' @method plot CalGrid
#' @export
plot.CalGrid <- function(x, runm=NA, calendar="BP", fill.p="grey50", border.p=NA, xlim=NA, ylim=NA, cex.lab=0.75, cex.axis=cex.lab, mar=c(4,4,1,3), add=FALSE,...){
yearsBP <- x$calBP
prob <- x$PrDens
calendars <- c("BP","BCAD")
if (!calendar %in% calendars){
stop("The calendar you have chosen is not currently an option.")
}
if (yearsBP[1] > yearsBP[length(yearsBP)]){
yearsBP <- rev(yearsBP)
prob <- rev(prob)
}
if (calendar=="BP"){
plotyears <- yearsBP
xvals <- c(plotyears[1],plotyears,plotyears[length(plotyears)], plotyears[1])
xlabel <- "Years cal BP"
} else if (calendar=="BCAD"){
plotyears <- BPtoBCAD(yearsBP)
xvals <- c(plotyears[1],plotyears,plotyears[length(plotyears)], plotyears[1])
xlabel <- "Years BC/AD"
if (all(range(plotyears)<0)){xlabel<-"Years BC"}
if (all(range(plotyears)>0)){xlabel<-"Years AD"}
} else {
stop("Unknown calendar type")
}
yvals <- c(0,prob,0,0)
if (calendar=="BP"){
xrng <- c(max(plotyears)+50, min(plotyears)-50)
xticks <- 200*(xrng%/%100 + as.logical(xrng%%200))
xticks <- seq(xticks[1]-200, xticks[2], -200)
} else {
xrng <- c(min(plotyears)-50, max(plotyears)+50)
xticks <- 200*(xrng%/%100 + as.logical(xrng%%200))
xticks <- seq(xticks[1]-200, xticks[2], 200)
}
if (is.na(ylim[1])){ ylim <- c(0,max(yvals*1.1)) }
if (is.na(xlim[1])){ xlim <- xrng }
if (!add){
par(mar=mar) #c(bottom, left, top, right)
par(cex.lab=cex.lab)
plot(xvals,yvals, type="n", xlab=xlabel, ylab="", xlim=xlim, ylim=ylim, xaxt='n', yaxt='n',axes=F, cex.axis=cex.axis,...)
}
xticksLab <- xticks
if (calendar=="BCAD"){
if (any(xticksLab==0)){xticksLab[which(xticksLab==0)]=1}
xticks[which(xticks>1)]=xticks[which(xticks>1)]-1
}
if (!add){
axis(1, at=xticks, labels=abs(xticksLab), las=2, cex.axis=cex.axis)
axis(2, cex.axis=cex.axis)
}
if (!is.na(runm)){ yvals <- runMean(yvals, runm, edge="fill") }
polygon(xvals,yvals, col=fill.p, border=border.p)
box()
}
#' @import stats
#' @import grDevices
#' @import graphics
#' @import utils
plot.UncalGrid <- function(x, type="adjusted", fill.p="grey50", border.p=NA, xlim=NA, ylim=NA, cex.lab=0.75, cex.axis=cex.lab, mar=c(4,4,1,3),...){
years <- x$CRA
if (type=="adjusted"){
prob <- x$PrDens
} else if (type=="raw"){
prob <- x$Raw
} else {
stop("Currently, type must be 'raw' or 'adjusted'.")
}
if (years[1] > years[length(years)]){
years <- rev(years)
prob <- rev(prob)
}
plotyears <- years
xvals <- c(plotyears[1],plotyears,plotyears[length(plotyears)], plotyears[1])
xlabel <- "C14 Years"
yvals <- c(0,prob,0,0)
xrng <- c(max(plotyears)+50, min(plotyears)-50)
if (dist(xrng)>10000){
xticks <- 1000*(xrng%/%1000 + as.logical(xrng%%1000))
xticks <- seq(xticks[1]-1000, xticks[2], -1000)
} else {
xticks <- 100*(xrng%/%100 + as.logical(xrng%%100))
xticks <- seq(xticks[1]-100, xticks[2], -100)
}
if (is.na(xlim[1])){ xlim <- xrng }
if (is.na(ylim[1])){ ylim <- c(0,max(yvals*1.1)) }
par(mar=mar) #c(bottom, left, top, right)
par(cex.lab=cex.lab)
plot(xvals,yvals, type="n", xlab=xlabel, ylab="", xlim=xlim, ylim=ylim, xaxt='n', yaxt='n', cex.axis=cex.axis,...)
axis(1, at=xticks, labels=abs(xticks), las=2, cex.axis=cex.axis)
axis(4, cex.axis=cex.axis)
polygon(xvals,yvals, col=fill.p, border=border.p)
}
#' @title Plot a stack of SPDs
#'
#' @description Visualises multiple SPDs grouped as a \code{stackCalSPD} object.
#'
#' @param x A \code{stackCalSPD} class object. Result of \code{\link{stackspd}} function.
#' @param type How to display the SPDs.Current options are \code{'stacked'},\code{'lines'}, '\code{'proportion'}. and \code{'multipanel'}. Default is \code{'stacked'}.
#' @param calendar Either \code{'BP'} or \code{'BCAD'}. Indicate whether the calibrated date should be displayed in BP or BC/AD. Default is \code{'BP'}.
#' @param spdnormalised A logical variable indicating whether the total probability mass of the SPDs are normalised to sum to unity. Default is FALSE.
#' @param rescale A logical variable indicating whether the summed probabilities values should be rescaled to range 0 to 1. Default is FALSE.Notice that this is different from setting \code{spdnormalised} to TRUE.
#' @param runm A number indicating the window size of the moving average to smooth the SPD. If set to \code{NA} no moving average is applied. Default is NA
#' @param xlim the x limits of the plot. In BP or in BC/AD depending on the choice of the parameter \code{calender}. Notice that if BC/AD is selected BC ages should have a minus sign (e.g. \code{c(-5000,200)} for 5000 BC to 200 AD).
#' @param ylim the y limits of the plot.
#' @param xaxt Whether the x-axis tick marks should be displayed (\code{xaxt='s'}, default) or not (\code{xaxt='n'}).
#' @param yaxt Whether the y-axis tick marks should be displayed (\code{xaxt='s'}, default) or not (\code{xaxt='n'}).
#' @param gapFactor Defines spacing between SPDs as proportion of the y-axis range for multipanel plots. Default is 0.2.
#' @param col.fill Vector of fill color for the observed SPDs. The default color scheme is based on the Dark2 pallette of RColorBrewer package.
#' @param col.line Line colour for the observed SPDs.The default color scheme is based on the Dark2 palette of RColorBrewer package.
#' @param lwd.obs Line width for the observed SPDs. Default is 1.
#' @param lty.obs Line type for the observed SPDs. Default is 1.
#' @param cex.lab The magnification to be used for x and y labels relative to the current setting of cex. Default is adjusted to 1.
#' @param cex.axis The magnification to be used for axis annotation relative to the current setting of cex. Default is adjusted to 1.
#' @param legend Whether legend needs to be displayed. Item names will be retrieved from the values supplied in the argument \code{group} in \code{\link{stackspd}}. Default is TRUE.
#' @param legend.arg list of additional arguments to pass to \code{\link{legend}}; names of the list are used as argument names. Only used if \code{legend} is set to TRUE. If supplied legend position must be given (e.g. \code{legend.arg=list(x='bottomright')}.
#' @param ylab a title for the y axis
#' @param ymargin multiplier for the maximum value on ylim range. Default is 1.1.
#' @param rnd integer indicating the number of decimal places to be displayed in the y-axis for when \code{type} is set "multitype".