R/smallMultipleWDPMKorcakPlot.R

Defines functions smallMultipleWDPMKorcakPlot

smallMultipleWDPMKorcakPlot <- function(datafile, minvalue){
  # compares and plots multiple data sets from WDPM
  # names of files and titles are read from file

  # declare plotting variables
  p <- NULL
  value <- NULL
  Model <- NULL

  # read in file of datasets
  datasets <- utils::read.csv(datafile, header=TRUE, stringsAsFactors = FALSE)
  datasets.count <- nrow(datasets)
  infile <- as.character(datasets$filename[1])
  caption <- as.character(datasets$caption[1])
  d1 <- utils::read.csv(infile, header=FALSE, sep=",")
  d1.count <- nrow(d1)
  d1 <- cbind(d1,rep.int(caption,nrow(d1)))
  names(d1) <- c("num","value","p","Model")

  d1 <- subset(d1,value >= minvalue)
  d1.korcak <- korcak(d1$value)
  d1.korcak <- subset(d1.korcak,p>0)
  d1.korcak$Model <- rep.int(caption,nrow(d1.korcak))


  # re-do Korcak for values >= 1
  for (i in 2:datasets.count){
    infile <- as.character(datasets$filename[i])
    caption <- as.character(datasets$caption[i])
    d2 <- utils::read.csv(infile, header=FALSE, stringsAsFactors = FALSE)
    d2.count <- nrow(d1)
    d2 <- cbind(d2,rep.int(caption,nrow(d2)))
    names(d2) <- c("num","value","p","Model")

    d2 <- subset(d2,value>=100)
    d2.korcak <- korcak(d2$value)
    d2.korcak <- subset(d2.korcak,p>0)
    d2.korcak$Model <- rep.int(caption,nrow(d2.korcak))

    # add d2 to d1
    d1.korcak <- rbind(d1.korcak, d2.korcak)

  }

  ybreaks <- c(0.0001,0.001,0.01,0.1,0.2,0.4,0.6,0.8,1)
  xbreaks <- c(10,100,1000)

  p1 <- ggplot2::ggplot(d1.korcak) +
    ggplot2::geom_point(ggplot2::aes(value,p,colour=Model,shape=Model)) +
    ggplot2::scale_colour_manual(values = c("blue","black","red")) +
    ggplot2::xlab(expression(paste('Water area (m',''^{2}, ')', sep = "")))+
    ggplot2::ylab("Exceedance probability")+
    ggplot2::scale_y_log10() +
    ggplot2::scale_x_log10() +
    ggplot2::stat_smooth(method = "lm", formula = y ~x+ I(x^2),
                         ggplot2::aes(value,p,color=Model), size=0.25,se=F) +
    ggplot2::coord_equal(ratio = 1)

  return(p1)

}
KevinShook/WDPMr documentation built on April 23, 2022, 12:32 a.m.