R/ausmaps.R

In haddonm/rforcpue: Functions to Assist With The Analysis of CPUE Data

#' @title addpoints includes a set of points with a set size and symbol
#'
#' @description addpoints given a matrix or data.frame, which must contain
#'   fields names 'Lat' and 'Long', these value pairs will be added to the plot
#'   generated by plotaus
#' @param indat a matrix or data.frame containing, at least, columns with the
#'   names 'Long' nad 'Lat'.
#' @param inpch default = 20, but can be changed to any valid symbol
#' @param incex default 0.2, but can be a constant or variable
#' @param incol default 2 (red) but can be any colour or rgb definition
#' @param intitle main top center title, the same as in plotaus; this is
#'   mostly useful if the refill variable is set to TRUE
#' @param jit default FALSE. If TRUE it will jitter the outcome by the
#'   'wobble' valriable amount. However, for overlapping points it is often
#'   better to use an rgb colour with the last value set to a fraction
#'   e.g. incol=rgb(9)0,0,0,1/5)
#' @param wobble the amount to jitter a point id jit = TRUE
#' @param refill default TRUE, refill the land after plotting points.
#' @param Long used to define the longitude field; defaults to 'Long'
#' @param Lat  used to define the Latitude field; default to 'Lat'
#' @param txtout logical determines whether the number of points added and \
#'     their range will be printed to the screen; defaults to TRUE
#' @param namecatch is the name of the catch variable; defaults to 'catch_kg'
#' @return adds the input points to a map; returns nothing.
#' @export addpoints
#' @examples
#' dev.new(height=6.0,width=7.5,noRStudioGD = TRUE)
#' plotaus()
#' Long <- c(140,141,142,143)
#' Lat <- c(-41.1,-41.1,-41.1,-41.1)
#' indata <- cbind(Long,Lat)
#' addpoints(indata,incex=1.5)
addpoints <- function(indat,inpch=20,incex=0.2,incol=2,intitle="",jit=F,
                      wobble=1, refill=T,Long="Long",Lat="Lat",txtout=TRUE,
                      namecatch="catch_kg") {
   if (length(grep("catch_kg",colnames(indat))) > 0)
      if (txtout)
         cat(length(indat[,"catch_kg"]),range(indat[,"catch_kg"],na.rm=T),"\n")
   if (jit) {
      points(jitter(indat[,Long],wobble),jitter(indat[,Lat],wobble),
             pch=inpch,cex=incex,col=incol)
   } else {
      points(indat[,Long],indat[,Lat],pch=inpch,cex=incex,col=incol)
   }
   ans <- c(802,2340,3732,6292,6302,6363,6686,7092,8128,8358,8488,8501)
   if (refill) {
      for (i in ans) {
         pick <- which(aus$poly == i)
         polygon(aus$Long[pick],aus$Lat[pick],col="lightblue")
      }
      mtext(intitle,side=3,outer=F,line=-1,cex=1.0,font=4)
   }
}  # end of addpoints

#' @title maps - lists the important functions with their syntax ready for use
#'
#' @description maps - lists the important functions with their syntax ready
#'   for use. Just copy the printed text and modify to suit your own purposes.
#' @return prints the useable functions and their syntax ready for use.
#' @export
#' @examples
#' maps()
maps <- function() {
   cat("leftlong <- 129;  rightlong <- 155 \n")
   cat("uplat <- -25;  downlat <- -44.6 \n")
   cat("plotaus(leftlong,rightlong,uplat,downlat,defineplot,gridon,infont) \n")
   cat("plotpolys(intmp,leftlong,rightlong,uplat,downlat,leg,mincount,intitle,gridon,hot,map,scaler) \n")
   cat("addpoints(indat,inpch,incex,incol,intitle,jit,wobble) \n")
   cat("plotLand()  \n")
}  # end of maps function

#' @title plotaus sets up a schematic map outline of the SESSF
#'
#' @description plotaus sets up a schematic map outline of Australia. Only the
#'     landmass contained within a box defined by leftlong, rightlong, uplat,
#'     and downlat. is included. The land areas is shaded.
#'
#' @param leftlong the longitude of the left hand side of the box plotted;
#'   defaults to 129.0
#' @param rightlong the longitude of the right hand size of the box plotted;
#'   defaults to 155.0
#' @param uplat the latitude (negative) of the top of the box; default -25.0
#' @param downlat nagative latitude of the bottom of the box; default -44.6
#' @param defineplot default = TRUE; declare the par values for a plot. If this
#'   is set to FALSE it would be possible to plot out multiple maps in one plot.
#' @param gridon plot lines at the value of 'gridon'. If set to 0.5 it would
#'   include lines every half degree.
#' @param maintitle default = "", adds a title to the plot at top center.
#' @param infont default = 7, bold Times. Used to define the font for axes, etc
#' @return plots a schematic map of the SESSF, If defineplot=TRUE it will alter
#'   the default par values for plots.
#' @export
#' @examples
#' \dontrun{
#'  dev.new(height=6.0,width=7.5,noRStudioGD = TRUE)
#'  plotaus()
#' }
plotaus <- function(leftlong=129,rightlong=155,uplat=-25,downlat=-44.6,
                      defineplot=T,gridon=0,maintitle="",infont=7) {
   ans <- c(802,2340,3732,6292,6302,6363,6686,7092,8128,8358,8488,8501)
   # xbound <- c(138.1333,138.1333,138.1333,144.0,144.0,144.0,146.0,147.0,147.0,
   #             147.0,148.0,148.0,157.0,148.0,148.0,148.31666,148.7333,148.7333,
   #             160.0,149.5,149.5,160.0)
   # ybound <- c(-34.2,-44.75,-40.0,-40.0,-38.5,-40.0,-42.0,-42.0,-44.75,
   #             -42.0,-41.0,-40.75,-40.75,-40.75,-38.4,-38.4,-37.8,-37.26,
   #             -37.26,-37.26,-33.58333,-33.58333)
   if (defineplot) {
      par(mfrow= c(1,1))
      par(mai=c(0.5,0.5,0.15,0.15), oma=c(0,0,0,0), xaxs="i",yaxs="i")
      par(cex=0.85, mgp=c(1.5,0.35,0), font.axis=infont,font=infont)
   }
   plot(c(leftlong,rightlong),c(uplat,downlat),type="n",ylim=c(downlat,uplat),
        xlim=c(leftlong,rightlong),xlab="",ylab="")
   if (gridon > 0) {
      abline(v=seq(trunc(leftlong),trunc(rightlong+1),gridon),lty=2,col="grey")
      abline(h=seq(trunc(downlat-1),trunc(uplat),gridon),lty=2,col="grey")
   }
   # lines(xbound,ybound,col=4)
   # segments(136.0,-34.5,136.0,-40.5,col=4)    # west of 85 east of 84
   # segments(133.0,-32.0,133.0,-39.0,col=4)  # west of 84 east of 83
   # segments(129.0,-31.0,129.0,-37.0,col=4)  # west bound of 83 east of 82
   # segments(121.0,-33.0,121.0,-37.0,col=4)
   # segments(152.0,-28.2,157.0,-28.2,col=4)
   # segments(157.0,-21,157.0,-33.58333,col=4)
   # segments(157.0,-21,157.0,-33.58333,col=4)
   # segments(157.0,-21,157.0,-33.58333,col=4)
   # segments(160.0,-33.58333,160.0,-35,col=4)

   for (i in ans) {
      pick <- which(aus$poly == i)
      polygon(aus$Long[pick],aus$Lat[pick],col="lightgrey")
   }
   if (defineplot) {
      title(ylab=list("Latitude S", cex=1.2, col=1, font=infont),
            xlab=list("Longitude E", cex=1.2, col=1, font=infont))
   }
   mtext(maintitle,side=3,outer=F,line=-1.0,cex=1.0,font=infont)
}  # end of plotaus

#' @title plotLand - puts land mass on top of a plot; covers land overlaps
#'
#' @description plotLand - puts land mass on top of a plot; covers land overlaps
#'   Used primarily after using plotpolys, which often leaves some oblongs
#'   sitting over parts of the coast. Using addpoints can also leave points
#'   on land. However, if you want to illustrate things on land then now you
#'   can use RGB = TRUE, and if you set alpha to 0.1 this will be transluscent
#'   andpoints in rivers, etc, will still be visible. If you use this option
#'   then incol becomes the red intensity and must lie between 0 and 1.
#' @param incol defaults to "lightgrey" but can be any single number representing
#'    a colour.
#' @param RGB should we use an rgb colour; default = FALSE
#' @param green the green intensity default = NA
#' @param blue the blue intensity
#' @param alpha the standard density 1 = full colour 1/10 = one tenth
#'
#' @return Nothing returned but the plot of land is repeated in light blue
#' @export
#' @examples
#' \dontrun{
#'   dev.new(height=6.0,width=7.5,noRStudioGD = TRUE)
#'   plotaus()
#'   plotLand(incol="red")
#'   plotLand(incol=1,RGB=TRUE,green=0,blue=0,alpha=(1/10))
#' }
plotLand <- function(incol="lightgrey",RGB=FALSE,
                     green=NA,blue=NA,alpha=NA) {
   ans <- c(802,2340,3732,6292,6302,6363,6686,7092,8128,8358,8488,8501)
   for (i in ans) {  # aus is an internal dataset to cede
      pick <- which(aus$poly == i)
      if (RGB) {
         polygon(aus$Long[pick],aus$Lat[pick],
                 col=rgb(incol,green,blue,alpha))
      } else {
         polygon(aus$Long[pick],aus$Lat[pick],col=incol)
      }
   }
} # end of plotLand




#' @title plotpoly adds a 'gridon' box to the map
#'
#' @description plotpoly adds a 'gridon' box to the map t is purely
#'     used within plotpolys and so is not exported
#' @param inlong top left longitude
#' @param inlat  top left latitude
#' @param col  the colour of the box polygon, defaults to black
#' @param gridon the size of the box, defaults to 0.1 degree
#'
#' @return nothing, the function adds a rectangular polygon to map
#'     it is not exported
plotpoly <- function(inlong,inlat,col=1,gridon=0.1) {
   #inlong and inlat represent the top left of a 6x6 box
   x <- c(inlong,inlong,inlong+gridon,inlong+gridon,inlong)
   y <- c(inlat,inlat-gridon,inlat-gridon,inlat,inlat)
   polygon(x,y,col=col,lty=1)
}

#' @title plotpolys - plots a grid and fills it with relative catches
#'
#' @description plotpolys - uses the same input object as plotaus plus the
#'    leftlong, rightlong, uplat, and downlat variables to which it adds a grid
#'    defined by 'gridon'. It then sums the catches found in each grid and then
#'    plots this as a patchwork of colours representing the relative density.
#'    Uses the function 'plotpoly'.
#' @param intmp - a matrix or data.frame containing, at least, columns with the
#'    names 'Long', 'Lat', and 'catch_kg'.
#' @param leftlong the longitude of the left hand side of the box plotted
#' @param rightlong the longitude of the right hand size of the box plotted
#' @param uplat the latitude (negative) of the top of the box
#' @param downlat nagative latitude of the bottom of the box
#' @param leg - legend location, default is 'topleft', alternatives are
#'    'topright', 'center', 'left', and 'right'
#' @param mincount - the minimum number of observations in a grod cell before it
#'    is included in the colour plotting; default - 1.0
#' @param intitle - default="". Provides for a title top center.
#' @param gridon - default 0.1 = 0.1 degree square. could be 0.2, 0.25, 0.5, etc
#' @param hot - default FALSE. If true then instead of 8 levels of catch
#'    intensity there are only 4 and these differ each by an order of magnitude.
#' @param map - default TRUE, use FALSE is only the tabulation of catch data into grid cells
#'    is required without a plot.
#' @param scaler - default to NA, otherwise used to select which scale to use
#'    valid values are 1 - 9. This is used to set the same scale in each of
#'    multiple plots.
#' @param textout default = TRUE; determines whether the summary information
#'    from the tabulation is pronted to screen or not
#' @param namecatch the name of the variable which to sum into each polygon.
#'     this defaults to 'catch_kg', but could be 'Effort', etc.
#' @return invisibly returns a matrix of grids with the total catch in each cell
#'    in units of tonnes. Also, if map = TRUE, adds the grid and a colour rendition
#'    to the map.
#' @export plotpolys
#' @examples
#' \dontrun{
#' plotaus()
#' }
plotpolys <- function(intmp,leftlong,rightlong,uplat,downlat,leg="topleft",
                      mincount=1,intitle=" ",gridon=0.1,hot=F,map=T,scaler=NA,
                      textout=TRUE,namecatch="catch_kg") {
   colnames(intmp) <- tolower(colnames(intmp))
   pick <- which((intmp$long < leftlong) | (intmp$long > rightlong))
   if (length(pick > 0)) { intmp <- intmp[-pick,] }
   pick <- which((intmp$lat > uplat) | (intmp$lat < downlat))
   if (length(pick > 0)) { intmp <- intmp[-pick,] }
   divis <- 1.0/gridon
   intmp$tlong <- trunc(intmp$long*divis)/divis
   intmp$tlat <- trunc(intmp$lat*divis)/divis
   east <- trunc(min(intmp$tlong,na.rm=T))
   west <- trunc(max(intmp$tlong,na.rm=T)+1)
   north <- trunc(max(intmp$tlat,na.rm=T))
   south <- trunc(min(intmp$tlat,na.rm=T)-1)
   cols <- seq(east,west,gridon)
   rows <- seq(north,south,-gridon)
   nlats <- length(rows)
   nlongs <- length(cols)
   totcat <- matrix(0,nrow=nlats,ncol=nlongs,dimnames=list(rows,cols))
   # polycol<- matrix(0,nrow=nlats,ncol=nlongs,dimnames=list(rows,cols))
   count  <- matrix(0,nrow=nlats,ncol=nlongs,dimnames=list(rows,cols))
   for (y in 1:nlats) {
      for (x in 1:nlongs) {
         long <- cols[x]
         lat <- rows[y]
         pick <- which((intmp$tlong == long) & (intmp$tlat == lat))
         if (length(pick) > 0) {
            count[y,x] <- length(pick)
            totcat[y,x] <- sum(intmp[pick,namecatch],na.rm=T)
         }
      }
   }
   totcat <- totcat/1000  # convert to tonnes
   maxcat <- max(totcat)
   subdivdat <- c(1000,800,700,500,250,100,50,1,
                  3000,2000,1000,750,500,250,100,1,
                  5000,3000,2000,1000,500,250,100,1,
                  10000,5000,2500,1000,500,250,100,1,
                  20000,10000,5000,2000,1000,500,100,1,
                  50000,25000,10000,5000,1000,500,100,1,
                  200000,100000,50000,20000,5000,1000,100,1,
                  500000,250000,100000,50000,10000,5000,1000,1,
                  1000000,500000,100000,50000,10000,5000,1000,1)
   subdivdat <- subdivdat/1000
   numclass <- length(subdivdat)/8  ## 8 possible scales
   subdivinfo <- matrix(subdivdat,nrow=numclass,ncol=8,byrow=TRUE)
   if (is.na(scaler)) {
      pick <- which(subdivinfo[,1] > maxcat)
   } else {
      if ((scaler <= numclass) & (scaler > 0)) {
         pick <- scaler
      } else {
         warning("Incorrect scale selected  \n")
         pick <- 9
      }
   }
   if (length(pick)) { subdiv <- subdivinfo[pick[1],] }
   else { subdiv <- subdivinfo[9,] }
   pick <- which(totcat > 0)
   lower <- quantile(totcat[pick],probs=(0.1))
   pick <- which(totcat > lower)
   upper <- round(quantile(totcat[pick],probs=(0.9)),3)
   subdiv2 <- c(upper,upper/10,upper/100,upper/1000)
   counthot <- numeric(4)
   hotcols <- c(1,2,3,0)
   countpoly <- numeric(8)
   for (y in 1:nlats) {
      for (x in 1:nlongs) {
         if (count[y,x] >= mincount) {
            if (hot) {
               pick2 <- which(subdiv2 < totcat[y,x])
               pickcol <- hotcols[pick2[1]]
               counthot[pick2[1]] <- counthot[pick2[1]] + 1 }
            else { pick2 <- which(subdiv < totcat[y,x])
            pickcol <- pick2[1]
            countpoly[pick2[1]]  <-  countpoly[pick2[1]] + 1
            }

            plotpoly(cols[x],rows[y],pickcol,gridon)
         }
      }
   }
   if (hot) {
      label <- paste(">",subdiv2,sep="")
      legcol <- hotcols
   } else {
      label <- paste(">",subdiv,sep="")
      legcol <- seq(1,8,1)
   }
   if (map) {
      ans <- c(802,2340,3732,6292,6302,6363,6686,7092,8128,8358,8488,8501)
      for (i in ans) {
         pick <- which(aus$poly == i)
         polygon(aus$Long[pick],aus$Lat[pick],col="lightgreen")
      }
   }
   plotLand()
   x <- leftlong; y <- uplat  # default topleft
   if (leg != "off") {
      if (leg == "topleft") { x <- leftlong + (rightlong-leftlong)/8; y <- uplat }
      if (leg == "topright") { x <- rightlong-2; y <- uplat }
      if (leg == "center") { x <- ((leftlong+rightlong)/2)-1; y <- uplat }
      if (leg == "left") { x <- leftlong; y <- ((uplat+downlat)/2) }
      if (leg == "right") { x <- rightlong-(rightlong-leftlong)/3; y <- downlat-((downlat-uplat)/2) }
      legend(x,y,label,col=legcol,lwd=6,bty="n",cex=0.65)
   }
   mtext(intitle,side=3,line=-1.0,outer=F,font=7,cex=1.0)
   if (textout) {
     cat("subdiv    ",subdiv2,"\n")
     cat("counthot  ",counthot,"\n")
     cat(maxcat,"   ",subdiv,"\n")
     cat("countpoly ",countpoly,"\n")
   }
   return(invisible(totcat))
} # end of plotpolys