R/SurveySim.R

In DIGSS: Determination of Intervals Using Georeferenced Survey Simulation

#'Survey Simulations
#'
#'Master function that runs survey simulations
#'
#'@details `surveySim()` is the main function that will conduct the simulations
#'of survey grids based on the Survey Parameters created by user.
#'@param survey.parameters List of class `surveySim`. See package's vignette or the help of `parametersExample` for details.
#'The following items must be included in list of class `surveySim`\itemize{
#'\item `col.width` the space between columns in the grid IN METERS
#'\item `grid.type` options are: "square","rectangle","staggered","hexagonal","arbitrary.staggered", following Kintigh 1988
#'\item `simulations` number of random maps to be created and contrasted with the grids
#'\item `area` vector with horizontal and vertical size of area surveyed in km
#'OBS: Sites will all be ellipses with radii not too different and random angles
#'\item `site.density` measured as number of sites/\ifelse{html}{\out{km<sup>2</sup>}}{\eqn{km^2}}. Can be either one value or a vector with 2 values (min and max) to create a range of densities
#'\item `site.area` can be one of two options: 1. one value indicating the area of all sites, in \ifelse{html}{\out{m<sup>2</sup>}}{\eqn{m^2}};
#' or 2. a vector with 4 values: min, max, mean (or median), and standard deviation in \ifelse{html}{\out{m<sup>2</sup>}}{\eqn{m^2}}
#'\item `overlap` maximum overlap of site area, ranging from 0 = no overlap allowed to 1 = complete overlap possible
#'\item `obj.density` artifacts per \ifelse{html}{\out{m<sup>2</sup>}}{\eqn{m^2}}. Can be a single value (uniform for all sites) or a range of values defined as min and max
#'\item `obj.distribution` type of cloud distribution for artifacts inside sites. Choose from: 'uniform', 'linear', 'spherical', 'sinusoidal'
#'\item `survey.radius` the radius of the survey pit (assumed to be a circle)
#'}
#'@param artifact.analysis Includes the analysis of artifacts in sites. Can be very time-consuming for mid-high artifact densities. Default = `TRUE`
#'@param plot If `TRUE` the last iteration of the simulations will be plotted.
#'@param plot.artifacts If `TRUE` the artifacts in each site will be plotted. Default = `FALSE` because it is resource intensive for mid-high densities.
#'@param areaprecision Area covered by sites is approximated by a cookie-cutter approach. Default precision = `1000`
#'gives approximate area within 1% error. Higher values reduce error but are more resource intensive.
#'
#'@return A list with four objects: \tabular{ll}{
#'    \code{Summary} \tab A matrix with summary statistics about number of surveys, frequency of site founds, artifacts presence, and success rate of simulations.\cr
#'    \tab \cr
#'    \code{BySite} \tab A matrix with results of the analyses by site from each of the survey areas created in the simulation. \cr
#'    \tab \cr
#'    \code{ByArtifact} \tab A matrix with results of the analyses by artifact presence in survey pits from each of the survey areas created in the simulation. \cr
#'    \tab \cr
#'    \code{Parameters} \tab A list with the parameters used to run the simulation. \cr
#'    }
#'#'@references
#'Kintigh (1988) The Effectiveness of Subsurface Testing: A Simulation Approach.
#'American Antiquity, 53:686-707. \doi{10.2307/281113}
#'@examples
#'\donttest{
#'#Runs simulations based on parametersExample
#'surveyExample<-surveySim(parametersExample)
#'
#'#Edit parametersExample to have 50 simulations and then run simulations
#'tmp_parameters <-parametersExample
#'tmp_parameters$simulations<-50
#'survey50<-surveySim(tmp_parameters)
#'}
#' @importFrom grDevices rainbow rgb
#' @importFrom graphics axis box lines plot.new plot.window points polygon text title
#' @importFrom stats quantile rnorm runif sd
#' @importFrom utils data flush.console
#'
#'@export

surveySim<-function(survey.parameters, artifact.analysis=TRUE, plot=TRUE, plot.artifacts=FALSE,areaprecision=1000){

  #We start with some checks and trying to make the function more user friendly
  if(class(survey.parameters)!="surveySim"){
    stop(cat("ERROR: Wrong Type of Object.","Consider using parametersExample as a starting point.",sep="\n"))
  }

  #Here we extract variables from the list (easier than call the list everytime)
  col.width<-survey.parameters$col.width
  grid.type<-survey.parameters$grid.type
  simulations<-survey.parameters$simulations
  Area<-survey.parameters$area
  site.density<-survey.parameters$site.density
  site.area<-survey.parameters$site.area
  overlap<-survey.parameters$overlap
  obj.density<-survey.parameters$obj.density
  obj.distribution<-survey.parameters$obj.distribution
  #I'm bringing this back to the scale of km2 already here
  survey.radius<-survey.parameters$survey.radius/1e3
  #this is for plotting the survey circle
  circle.x<-seq(-survey.radius,survey.radius,length=40)
  circle.y<-sqrt(survey.radius^2-circle.x^2)

  circle.x<-c(circle.x,circle.x[length(circle.x):1])
  circle.y<-c(circle.y,-circle.y)

  #this is for the progress bar
  char="*" #this is the character to be plotted in the progressbar
  charwidth<-nchar(char,"width")

  totalwidth<-getOption("width")-5
  info.width<-nchar("Simulation: XXX% | Sites: XXX%")
  width<-trunc(totalwidth/charwidth)-info.width
  totalwidth<-trunc(totalwidth/charwidth)+5

  GRID<-c("square","rectangle","staggered","hexagonal","arbitrary.staggered")

  grid.type<-pmatch(grid.type,GRID)

  if (is.na(grid.type)){
    stop(cat("ERROR: invalid grid type.","Valid grids are: square,rectangle,staggered,hexagonal,arbitrary.staggered",sep="\n"))
  }

  #1. Build the grid
  ##We will get the number of rows t, i (interval between units in a row),s(space between columns),
  ## and e (Distance from transect to survey area edge)
  s<-col.width/1000

  t<-floor(Area[1]/s)+1

  e<-(Area[1]%%(s))/2

  if(grid.type==1){ #square

    i=s
  }
  if(grid.type==2){ #rectangle

    tmp<-as.numeric(readline(prompt = "Type  the ratio between length and width of the rectangles in the grid:"))

    i=s*tmp
  }
  if(grid.type==3){ #staggered

    i=s
  }
  if(grid.type==4){ #hexagonal

    i=(3^0.5/2)*s
  }
  if(grid.type==5){ #arbitrary staggered

    tmp<-as.numeric(readline(prompt = "Type  the ratio between length and width of the rectangles in the grid:"))
    i=s*tmp
  }

  #2. We create the grid (or rather the points of intersection)

  xseq<-seq(0,Area[1],by=s)
  yseq<-seq(0,Area[2],by=i)

  #.this centralizes the grid, so that the edges are even in the field
  xseq<-xseq+(Area[1]-max(xseq))/2
  yseq<-yseq+(Area[2]-max(yseq))/2

  #i.this is the object that will store final summary stats
  SurveySummary<-list("Summary"=matrix(NA,5,6),
                      "BySite"=matrix(0,simulations,10),
                      "ByArtifact"=matrix(NA,simulations,7),
                      "Parameters"=list("col.width"=col.width,
                                        "grid.type"=grid.type,
                                        "simulations"=simulations,
                                        "Area"=Area,
                                        "site.density"= site.density,
                                        "site.area"=site.area,
                                        "overlap"=overlap,
                                        "obj.density"=obj.density,
                                        "obj.distribution"=obj.distribution,
                                        "survey.radius"=survey.radius*1000))

  colnames(SurveySummary$Summary)<-c("Mean","StDev","Min","Max","Quantile 2.5%","Quantile 97.5%")
  row.names(SurveySummary$Summary)<-c("SurveysPerSim","SitesFound%","SitesFoundOnArtifacts%","ArtifactsPerSurvey","SuccessRateIndex")

  colnames(SurveySummary$BySite)<-c("SurveyHits","TotalSurveys","SitesFound",
                                    "TotalSites","SitesFound%","SurveyFind%",
                                    "MaxSurveys/Site","AvgSurveys/Site","RealSiteArea",
                                    "SuccessRateIndex")
  colnames(SurveySummary$ByArtifact)<-c("SurveyHits","AvgArtifacts/Survey","SitesFound","TotalSites","SitesFound%","SurveyFind%","AvgArtifactDensity")

  #this is the number of surveys in the grid, to use in later calculations
  n.pits<-length(xseq)*length(yseq)
  if(grid.type==3||grid.type==4||grid.type==5){
    n.pits=n.pits-floor(length(yseq)/2)
  }

  #3. We create now the loop that will create all the maps and grids
  for (a in 1:simulations){
    sitemap<-fieldMap(Area,site.density,site.area,overlap,plot=FALSE,areaprecision=areaprecision)
    site.frame<-sitemap$site.frame

    ###This is to set up the progress bar
    value<-round(width*((a-1)/simulations),0)
    perc.ite<-round((a-1)/simulations*100,0)

    #3a. We will do the analysis per site, so we don't need to store artifacts in the memory.
    ##i. create a matrix to store how many sites each survey hits

    SurveyHitsMatrix<-matrix(0,length(xseq),length(yseq))
    SurveysPerSite<-rep(0,nrow(site.frame))

    if(artifact.analysis==TRUE){
      ArtiHitsMatrix<-matrix(0,length(xseq),length(yseq))
      SurveysWithArtiPerSite<-rep(0,nrow(site.frame))
      #this will create the vector for unique object densities per site, if use selects a range
      obj.density.vector<-rep(NA,nrow(site.frame))
    }

    #Plot1 - here we start the plot
    if(plot==TRUE & a==simulations){
      plot.new()
      plot.window(c(0,Area[1]),c(0,Area[2]),asp=Area[1]/Area[2])

      axis(1,at=round(xseq,2))
      axis(2,at=round(yseq,2))
      box()
      title(main=paste("Sites from simulation",a,"of",simulations),xlab="km",ylab="km")
    }

    for (b in 1:nrow(site.frame)){
      cat("\r",rep.int(" ",totalwidth),sep="")
      cat("\r   |",rep.int(char,value),rep.int(" ",charwidth*(width-value)),
          "| Simulation:",perc.ite,"% | Sites: ",round(b/nrow(site.frame)*100,0),"%",sep="")
      flush.console()

      #ii.create a tmp matrix to see which surveys hit this specific site
      tmpHitMatrix<-matrix(FALSE,length(xseq),length(yseq))

      #iii. test matrices.
      for (c in 1:length(yseq)){

        #iv. this if is to fix the length of the x positions for hexagonal grids.
        if(grid.type==3||grid.type==4||grid.type==5){
          if(b%%2==0){
            xvals=xseq+0.5*s
            xvals<-xvals[-length(xvals)]
          } else {
            xvals=xseq
          }
        }else{
          xvals=xseq
        }
        ##I'm adding radius of pit to a and b, to increase the site accordingly and test if any
        ## part of the pit hits the site.
        tmpHit<-((xvals-site.frame[b,5])*cos(site.frame[b,4])+(yseq[c]-site.frame[b,6])*sin(site.frame[b,4]))^2/(site.frame[b,7]+survey.radius)^2+
          ((xvals-site.frame[b,5])*sin(site.frame[b,4])-(yseq[c]-site.frame[b,6])*cos(site.frame[b,4]))^2/(site.frame[b,8]+survey.radius)^2<=1

        tmpHitMatrix[1:length(xvals),c]<-tmpHit
      }

      #v.here we add the stats for this site
      SurveysPerSite[b]<-sum(tmpHitMatrix)
      SurveyHitsMatrix<-SurveyHitsMatrix+tmpHitMatrix

      ##vi. if you want to look at artifact counts, here is where we start the madness.
      if(artifact.analysis==TRUE){
        ##. we will only generate the cloud of artifacts if strictly necessary

        obj.density.vector[b]<-runif(1,obj.density[1],obj.density[length(obj.density)])

        if(sum(tmpHitMatrix)>0|plot.artifacts==TRUE & a==simulations){
          if(plot.artifacts==TRUE & a==simulations){
            plot.now=TRUE
          }else{
            plot.now=FALSE
          }

          artifacts<-cloudGenerator(obj.density.vector[b],site.frame[b,7],site.frame[b,8],site.frame[b,4],
                                    site.frame[b,5],site.frame[b,6],type=obj.distribution,plot=plot.now)
        }

        ##this is the matrix that will store the total number of artifacts per survey
        tmpArtiCount<-matrix(0,length(xseq),length(yseq))

        ##vii. now we go and see how many survey pits find how many pieces.
        for(c in 1:length(yseq)){
          for(d in 1:length(xvals)){
            if(sum(tmpHitMatrix)>0){
              for (e in 1:length(artifacts$coords)){
                if(nrow(artifacts$coords[[e]])!=0){
                  tmpBelongIn<-((xvals[d]-site.frame[b,5])*cos(site.frame[b,4])+(yseq[c]-site.frame[b,6])
                                *sin(site.frame[b,4]))^2/(site.frame[b,7]*(e/length(artifacts$coords))+survey.radius)^2+
                    ((xvals[d]-site.frame[b,5])*sin(site.frame[b,4])-(yseq[c]-site.frame[b,6])
                     *cos(site.frame[b,4]))^2/(site.frame[b,8]*(e/length(artifacts$coords))+survey.radius)^2<=1

                  if(e>1){
                    tmpBelongOut<-((xvals[d]-site.frame[b,5])*cos(site.frame[b,4])+(yseq[c]-site.frame[b,6])
                                   *sin(site.frame[b,4]))^2/(site.frame[b,7]*((e-1)/length(artifacts$coords))-survey.radius)^2+
                      ((xvals[d]-site.frame[b,5])*sin(site.frame[b,4])-(yseq[c]-site.frame[b,6])
                       *cos(site.frame[b,4]))^2/(site.frame[b,8]*((e-1)/length(artifacts$coords))-survey.radius)^2>1
                  }else{
                    tmpBelongOut<-TRUE
                  }

                  if(tmpBelongIn==TRUE&tmpBelongOut==TRUE){
                    artsIN<-((artifacts$coords[[e]][,1]-xvals[d])^2+(artifacts$coords[[e]][,2]-yseq[c])^2)^0.5<=survey.radius
                    artN<-sum(artsIN)

                    tmpArtiCount[d,c]<-tmpArtiCount[d,c]+artN
                  }
                }
              }
            }
            if(tmpArtiCount[d,c]>0){
              SurveysWithArtiPerSite[b]<-SurveysWithArtiPerSite[b]+1
            }
          }
        }

        ArtiHitsMatrix<-ArtiHitsMatrix+tmpArtiCount
      }

      #Plot2: here we continue the plot, byt plotting the sites. Artifacts should have been plotted by now above
      if(plot==TRUE & a==simulations){

        #1.Get the angles, x and y to plot the ellypses
        angles<-seq(0,2*pi,length=72)

        #x= h + a cos(t)*cos(c)-b*sin(t)*sin(c)
        xcoords<- site.frame[b,5]+site.frame[b,7]*cos(angles)*cos(site.frame[b,4])-site.frame[b,8]*sin(angles)*sin(site.frame[b,4])
        #y= k + b sin(t)*cos(c)+a*cos(t)*sin(c)
        ycoords<- site.frame[b,6]+site.frame[b,8]*sin(angles)*cos(site.frame[b,4])+site.frame[b,7]*cos(angles)*sin(site.frame[b,4])

        polygon(xcoords,ycoords,border=rgb(0,0,1,),lty=3)

        text(site.frame[b,5],site.frame[b,6],b,cex=1)
      }
    }

    #4. Get the summary stats
    #i. By Site
    SurveySummary$BySite[a,1]<-length(which(SurveyHitsMatrix>0)) # #surveyhits
    SurveySummary$BySite[a,2]<-n.pits # #Totalsurveys
    SurveySummary$BySite[a,3]<-length(which(SurveysPerSite>0)) #Sites found
    SurveySummary$BySite[a,4]<-nrow(site.frame) #Total sites
    SurveySummary$BySite[a,5]<-SurveySummary$BySite[a,3]/SurveySummary$BySite[a,4] # %sites found
    SurveySummary$BySite[a,6]<-SurveySummary$BySite[a,1]/SurveySummary$BySite[a,2] # % survey find
    SurveySummary$BySite[a,7]<-max(SurveysPerSite) # Max surveys/site
    SurveySummary$BySite[a,8]<-mean(SurveysPerSite) # avg surveys / site
    SurveySummary$BySite[a,9]<-sitemap$actualArea #Real Site area
    SurveySummary$BySite[a,10]<-SurveySummary$BySite[a,6]/SurveySummary$BySite[a,9] #SuccessRateIndex

    #ii. By Artifact
    if(artifact.analysis==TRUE){
      SurveySummary$ByArtifact[a,1]<-length(which(ArtiHitsMatrix>0)) #survey hits
      SurveySummary$ByArtifact[a,2]<-ifelse(SurveySummary$ByArtifact[a,1]>0,sum(ArtiHitsMatrix)/SurveySummary$ByArtifact[a,1],0)#Avg Artifacts per survey
      SurveySummary$ByArtifact[a,3]<-length(which(SurveysWithArtiPerSite>0)) #sitesfound
      SurveySummary$ByArtifact[a,4]<-nrow(site.frame) #Total Sites
      SurveySummary$ByArtifact[a,5]<-SurveySummary$ByArtifact[a,3]/SurveySummary$ByArtifact[a,4] #%Sites found
      SurveySummary$ByArtifact[a,6]<-SurveySummary$ByArtifact[a,1]/n.pits #%Survey finds
      SurveySummary$ByArtifact[a,7]<-mean(obj.density.vector)
    }

    #Plot 3: Finish plot by adding the grid:
    if(plot==TRUE & a==simulations){
      for(r in 1:length(yseq)){
        xvals = xseq
        nextxvals = xseq

        if(grid.type==3||grid.type==4||grid.type==5)
          if(r%%2==0){
            xvals=xseq+0.5*s
            xvals<-xvals[-length(xvals)]
          }else{
            #this is for me to plot the lines that connect the survey dots
            nextxvals=nextxvals+0.5*s
            nextxvals<-nextxvals[-length(xvals)]
          }
        points(xvals,rep(yseq[r],length(xvals)),pch=16,cex=0.5)
        for (q in 1:length(xvals)){
          polygon(circle.x+xvals[q],circle.y+yseq[r],lwd=2)
        }

        #here I plot lines for the grids, to make them look nicer
        lines(c(xvals[1],xvals[length(xvals)]),c(yseq[r],yseq[r]),lty=2)

        if(r!=length(yseq)){
          for(b in 1: length(xvals)){
            if(is.na(nextxvals[b])!=TRUE){
              lines(c(xvals[b],nextxvals[b]),c(yseq[r],yseq[r+1]),lty=2)
            }
          }
        }
        #and this is just to add lines to the hexagon to make the grid look like triangles
        if(grid.type==4){
          if(r%%2==1){
            for (b in 2:length(xvals)){
              lines(c(xvals[b],nextxvals[b-1]),c(yseq[r],yseq[r+1]),lty=2)
            }

          }else{
            for (b in 1:length(xvals)){
              lines(c(xvals[b],nextxvals[b+1]),c(yseq[r],yseq[r+1]),lty=2)
            }
          }
        }
      }
    }
  }

  #Here we get the final Summary stats:
  #iii. The Total summaries
  ##SurveysPersim
  SurveySummary$Summary[1,1]<-mean(SurveySummary$BySite[,2])

  ##SitesFound%
  SurveySummary$Summary[2,1]<-mean(SurveySummary$BySite[,5])
  SurveySummary$Summary[2,2]<-sd(SurveySummary$BySite[,5])
  SurveySummary$Summary[2,3]<-min(SurveySummary$BySite[,5])
  SurveySummary$Summary[2,4]<-max(SurveySummary$BySite[,5])
  SurveySummary$Summary[2,5]<-quantile(SurveySummary$BySite[,5],0.025)
  SurveySummary$Summary[2,6]<-quantile(SurveySummary$BySite[,5],0.975)

  ##SitesFoundon Artifacts%
  SurveySummary$Summary[3,1]<-mean(SurveySummary$ByArtifact[,5])
  SurveySummary$Summary[3,2]<-sd(SurveySummary$ByArtifact[,5])
  SurveySummary$Summary[3,3]<-min(SurveySummary$ByArtifact[,5])
  SurveySummary$Summary[3,4]<-max(SurveySummary$ByArtifact[,5])
  SurveySummary$Summary[3,5]<-quantile(SurveySummary$ByArtifact[,5],0.025)
  SurveySummary$Summary[3,6]<-quantile(SurveySummary$ByArtifact[,5],0.975)

  ##ArtifactsPerSurvey
  SurveySummary$Summary[4,1]<-mean(SurveySummary$ByArtifact[,2])
  SurveySummary$Summary[4,2]<-sd(SurveySummary$ByArtifact[,2])
  SurveySummary$Summary[4,3]<-min(SurveySummary$ByArtifact[,2])
  SurveySummary$Summary[4,4]<-max(SurveySummary$ByArtifact[,2])
  SurveySummary$Summary[4,5]<-quantile(SurveySummary$ByArtifact[,2],0.025)
  SurveySummary$Summary[4,6]<-quantile(SurveySummary$ByArtifact[,2],0.975)

  ##SuccessRateIndex
  SurveySummary$Summary[5,1]<-mean(SurveySummary$BySite[,10])
  SurveySummary$Summary[5,2]<-sd(SurveySummary$BySite[,10])
  SurveySummary$Summary[5,3]<-min(SurveySummary$BySite[,10])
  SurveySummary$Summary[5,4]<-max(SurveySummary$BySite[,10])
  SurveySummary$Summary[5,5]<-quantile(SurveySummary$BySite[,10],0.025)
  SurveySummary$Summary[5,6]<-quantile(SurveySummary$BySite[,10],0.975)

  #this finishes the progress bar
  cat("\r   |",rep.int(char,width),
      "| Simulation:",100,"% | Sites: ",100,"%",sep="")
  cat("\n")
  flush.console()

  return(SurveySummary)
}