R/spatial_utils.R

In haddonm/abspatial: A set of functions to assist mapping Abalone Spatial data

# 
comptwo <- function(yr1,yr2) {  # yr1=year1; yr2=year2
  combined <- yr1[(which(yr1[,"oid"] %in% yr2[,"oid"])),]
  yr1only <- yr1[(which(yr1[,"oid"] %ni% yr2[,"oid"])),]
  yr2only <- yr2[(which(yr2[,"oid"] %ni% yr1[,"oid"])),]
  return(c(y1=dim(yr1only)[1],both=dim(combined)[1],y2=dim(yr2only)[1]))
}

# yr1=previous  yr2=current   yr3=next
# 
compthree <- function(yr1,yr2,yr3) {  # yr1=year1; yr2=year2
  pc_pcn <- yr1[(which(yr1[,"oid"] %in% yr2[,"oid"])),]
  cn_pcn <- yr2[(which(yr2[,"oid"] %in% yr3[,"oid"])),]
  pn_pcn <- yr1[(which(yr1[,"oid"] %in% yr3[,"oid"])),]
  pcn <- length(which(pn_pcn[,"oid"] %in% pc_pcn[,"oid"]))

  pc <- (dim(pc_pcn)[1] - pcn)
  cn <- (dim(cn_pcn)[1] - pcn)
  pn <- (dim(pn_pcn)[1] - pcn)

  allcn <- unique(c(yr3[,"oid"],yr2[,"oid"]))
  p <- length(which(yr1[,"oid"] %ni% allcn))

  allpc <- unique(c(yr1[,"oid"],yr2[,"oid"]))
  n <- length(which(yr3[,"oid"] %ni% allpc))

  allpn <- unique(c(yr1[,"oid"],yr3[,"oid"]))
  c <- length(which(yr2[,"oid"] %ni% allpn))

  answer <- c(p,pc,pn,pcn,c,cn,n)
  return(answer)
} # end of compthree


#' @title getblock selects the given block data from input data.frame
#'
#' @description getblock literally selects the data relating to the
#'     input blocks. The function stops with a warning if an invalid
#'     block is requested. It converts the oid column to numeric and adds
#'     a repeat column counting the years in which catches was taken
#'
#' @param inblock a vector of blocks or a block as a number
#' @param indat the abalone data.frame from which to select the data
#' @param indices identifies the column indices of the years within the
#'     data over which to summarize the data for repeated years
#'
#' @return a data.frame containing records only for the requested blocks
#' @export
#'
#' @examples
#' \dontrun{
#'  # ab <- readRDS(paste0(datadir,"G1Ha_2018_10_23.rds"))#get spatial obj
#'  # abdat <- getlatlong(ab)     # convert to data.frame + lat long
#'  # addat <- getblock(24,abdat)
#'  # addat <- getblock(124,abdat)
#' }
getblock <- function(inblock,indat,indices=c(5:11)) {
  inblock <- as.numeric(inblock)
  blks <- c(1:49,51:57)
  if (any(inblock %ni% blks))
    stop(cat("All members of ",inblock,"must be in ",blks," \n"))
  pick <- which(indat[,"blockno"] %in% inblock)
  addat <- droplevels(indat[pick,])
  addat[,"oid"] <- as.numeric(addat[,"oid"])
  addat[,"repeat"] <- apply(addat[,indices],1,countgtzero)
  return(invisible(addat))
}  # end of getblock



#' @title getlatlong converts a SpatialPolygonsDataFrame to a data.frame
#'
#' @description getlatlong takes a SpatialPolygonsDataFrame, from the sp
#'     package, and converts it into a classic data.frame and includes
#'     the long and lat of each record as the last two columns.
#'     Thanks to Craig Mundy for this function.
#'
#' @param x a SpatialPolygonsDataFrame, from the sp package
#'
#' @return a data.frame equivalent to x_at_data plus the long and lat
#'     as the last two columns
#' @export
#'
#' @examples
#' \dontrun{
#'  # ab <- readRDS(paste0(datadir,"G1Ha_2018_10_23.rds"))#get spatial obj
#'  # abdat <- getlatlong(ab)     # convert to data.frame + lat long
#' }
getlatlong <- function(x) {   #  x must be a SpatialPolygonsDataFrame
  wgs84projstring <- sp::CRS("+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs")
  sp::proj4string(x)
  LatLong <- sp::spTransform(x, wgs84projstring)
  LL <- cbind(as.data.frame(LatLong),sp::coordinates(LatLong))
  coln <- dim(LL)[2]  # explicitly name the last two columns
  colnames(LL)[(coln-1):coln] <- c("long","lat")
  return(LL)
} # end of getlatlong

#' @title getpoints selects any lat-long records within a view
#'
#' @description given a particular view (that is the leftlong, the
#'     rightlong, the uplat, and the downlat - as a vector of 4)
#'     getpoints will select those records from a data.frame
#'     being explored. This prepares the data for summarizing effort,
#'     catch, divers, etc.
#'
#' @param view the leftlong, rightlong, uplat, and the downlat as a
#'     vector of 4 which identifies the window onto Tasmania being
#'     considered.
#' @param dat a data.frame usually obtained from getlatlong applied to a
#'     SpatialPolygonsDataFrame. But it could be any that had a long and
#'     a lat field within it.
#' @param Long the name of the longitude column in the data.frame. defaults
#'     to long
#' @param Lat the name of the latitude column in the data.frame. defaults
#'     to lat
#'
#' @return a data.frame containing the specific records associated with the
#'     lat-long records within the area of the view.
#' @export
#'
#' @examples
#' \dontrun{
#'  # ab <- readRDS(paste0(datadir,"G1Ha_2018_10_23.rds"))#get spatial obj
#'  # abdat <- getlatlong(ab)     # convert to data.frame + lat long
#' }
getpoints <- function(view,dat,Long="long",Lat="lat") {
  left <- view[1]; right <- view[2]
  up <- view[3]; down <- view[4]
  pts <- NULL
  pick <- which((dat[,Long] >= left) & (dat[,Long] <= right) &
                  (dat[,Lat] <= up) & (dat[,Lat] >= down))
  if (length(pick) > 0) pts <- dat[pick,]
  return(pts)
} # end of getpoints

# pull out one year's columns of data from the total dataset.
#
getrecs <- function(pickcol,yr,adat,label,first=c(1:4),last=c(51:53)) {
  ydat <- adat[,c(first,pickcol[yr],last)]
  ydat <- droplevels(ydat[which(ydat[,label[yr]] > 0),])
  return(ydat)
}


#
getrepeats <- function(adat,indices=c(5:11)) {
  repeats <- apply(adat[,indices],1,countgtzero)
  counts <- as.matrix(table(repeats))
  answer <- cbind(counts,counts/sum(counts))
  colnames(answer) <- c("counts","props")
  return(answer)
}

#' @title getsubblock selects the given subblock data from input data.frame
#'
#' @description getsubblock literally selects the data relating to the
#'     input subblocks. The function stops with a warning if an invvalid
#'     subblock is requested.
#'
#' @param insubblk a vector of subblocks (or a single) as character
#' @param indat the abalone data.frame from which to selec the data
#' @param indices identifies the column indices of the years within the
#'     data over which to summarize the data for repeated years
#'
#' @return a data.frame containing records only for the requested subblocks
#' @export
#'
#' @examples
#' \dontrun{
#'  # ab <- readRDS(paste0(datadir,"G1Ha_2018_10_23.rds"))#get spatial obj
#'  # abdat <- getlatlong(ab)     # convert to data.frame + lat long
#'  # addat <- getsubblock(c("13C","13D","13E"),abdat)
#'  # head(addat,10)
#' }
getsubblock <- function(insubblk,indat,indices=c(5:11)) {
  sublks <- unique(indat[,"subblockno"])
  if (any(insubblk %ni% sublks))
    stop(cat("All members of ",insubblk,"must be in ",sublks," \n"))
  pick <- which(indat[,"subblockno"] %in% insubblk)
  addat <- droplevels(indat[pick,])
  addat[,"oid"] <- as.numeric(addat[,"oid"])
  addat[,"repeat"] <- apply(addat[,indices],1,countgtzero)
  return(invisible(addat))
}  # end of getblock

#' @title mapgrid includes a given variable into a plottas grpahic
#'
#' @description mapgrid is used when plotting up the distribution of
#'     different variables from a SpatialPolygonsDataFrame that has been
#'     translated into an ordinary data.frame. It uses the lat and long
#'     to place the given variable onto a plottas graphic. The variables
#'     available are blkg2012 - blkg2018, blkgtotal, days2012 -days2018,
#'     daystotal, mins2012 - mins2018, minstotal, divers2012-divers2018,
#'     kgday2012 - kgday2018. The values for each variable are all
#'     scaled relative to the maximum but can be multiplied by mult.
#'     Points will only be plotted where there is a value. An NA will
#'     result in a missing point. If only the location of an active grid
#'     is wanted then set loconly = TRUE, otherwise the relative value
#'     in x will be used to set the size of the point up to a maximum
#'     size of mult.
#'
#' @param view the length 4 vector with c(left,right,up, down) that acts
#'     as a window over Tasmania; if view = NA then the data from the input
#'     data.frame are plotted rather than using view to select the data.
#' @param dat the data.frame from which to select the 1 hectare grids
#'     within the given view
#' @param x the name of the variable selected; default = "days2012"
#' @param mult the multiplier on the size of the point plotted; default
#'     = 2
#' @param col the colour of each point; default = 2 (red)
#' @param loconly a logical value determining whether to plot only the
#'     location or a proxy for the relative value of the chosen x; the
#'     default = FALSE so mapgrid will plot the relative value
#' @param infill the background colour filling open symbols like pch=1;
#'     defaults to incol
#' @param inpch default = 21, but can be changed to any valid symbol
#'
#' @return invisibly returns the records selected by getpoints
#' @export
#'
#' @examples
#' \dontrun{
#'  # ab <- readRDS(paste0(datadir,"G1Ha_2018_10_23.rds"))#get spatial obj
#'  # abdat <- getlatlong(ab)     # convert to data.frame + lat long
#' }
mapgrid <- function(view,dat,x="days2012",mult=2.0,col=2,loconly=FALSE,
                    infill=col,inpch=21) {
  # dat=abdat; x="blkg2012"; view = maria; mult=2; col=2;loconly=FALSE; infill=col; inpch=21
  #if (is.numeric(view)) {
    addat <- getpoints(view,dat)
  #} else {
  #  addat <- dat
  #}
  pickcol <- which(colnames(addat) == x)
  if (length(pickcol) == 0)
    stop(cat("Selected variable not in data.frame  \n"))
  scalept <- mult*addat[,x]/max(addat[,x],na.rm=TRUE)
  if (loconly) scalept[scalept>0] <- 0.2
  mappoints(view,addat,inpch=inpch,incex=scalept,incol=col,infill=infill)
  return(invisible(addat))
} # end of mapgrid

#' @title mapLand - puts land mass on top of a plot; covers land overlaps
#'
#' @description mapLand - 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 mappoints can also leave points
#'   on land.
#'
#' @param col defaults to "lightgrey" but can be any colour
#' @param border the colour of the polygon outline; defaults to NULL
#'
#' @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)
#'  maptas()
#'  mapLand(incol="red")
#' }
mapLand <- function(col="lightgrey",border=NULL) {
  #ans <- c(802,2340,3732,6292,6302,6363,6686,7092,8128,8358,8488,8501)
  ans <- unique(taspoly$poly)
  for (i in ans) {
    pick <- which(taspoly$poly == i)
    polygon(taspoly$Long[pick],taspoly$Lat[pick],col=col,border=border)
  }
} # end of mapLand


#' @title mappoints includes a set of points with a set size and symbol
#'
#' @description mappoints 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 maptas
#'
#' @param view is a vector of four values defining the leftlong, the
#'     rightlong, the uplat, and the downlat. These define left, right,
#'     top, and bottom of the box plotted; defaults to c(129,155,-25,-44.6).
#'     The latitude values need to be negative (southern hemisphere).
#' @param indat a matrix or data.frame containing, at least, columns with
#'   names 'long' and 'lat'.
#' @param inpch default = 21, 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 plotsessf; this is
#'   mostly useful if the refill variable is set to 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 infill the background colour filling open symbols like pch=1;
#'     defaults to incol
#'
#' @return adds the input points to a map; returns nothing.
#' @export
#'
#' @examples
#' \dontrun{
#'  dev.new(height=6.0,width=6.0,noRStudioGD = TRUE)
#'  view=c(143.5,149,-39,-44.0)
#'  maptas(view=view)
#'  long <- c(147.4,147.8,148,148.2)
#'  lat <- c(-43.6,-43.6,-43.3,-43.2)
#'  indata <- cbind(long,lat)
#'  mappoints(view,indata,incex=1.5)
#' }
mappoints <- function(view,indat,inpch=21,incex=0.2,incol=2,intitle="",
                      refill=TRUE,Long="long",Lat="lat",infill=incol) {
  pdat <- getpoints(view,indat)
  points(pdat[,Long],pdat[,Lat],pch=inpch,cex=incex,col=incol,
         bg=infill)
  mtext(intitle,side=3,outer=F,line=-1,cex=1.0,font=4)
}  # end of mappoints

#' @title maptas sets up a schematic map outline of the SESSF
#'
#' @description maptas sets up a schematic map outline of Tasmania. This is
#'     contained within a box defined by leftlong, rightlong, uplat, and
#'     downlat, four values held in the view vector. The land areas is
#'     shaded.
#'     This all works by plotting out 12 polygons defined in sessf,
#'     although these include the whole of Australia. Each polygon has
#'     a number:
#'     802 = QLD; 2340  NT; 3732 = WA; 6292 King Island  6302 Flinders;
#'     6363 NewSouthWales;     6686 Lady Barron;    7092 Tasmania
#'     8128 Kangaroo Island    8358 South Australia 8488 Victoria
#'     8502 Bruny Island.    c(6292,6302,6686,7092,8488,8501) are
#'     used to generate a rapid plot of Tasmania. There are about 1036
#'     polygons in taspoly
#'
#' @param view is a vector of four values defining the leftlong, the
#'     rightlong, the uplat, and the downlat. These define left, right,
#'     top, and bottom of the box plotted; defaults to c(129,155,-25,-44.6).
#'     The latitude values need to be negative (southern hemisphere).
#' @param defineplot should the built in par statement be used; default
#'     = TRUE. If FALSE expects a par statement outside. This is used
#'     if you want to plot multiple maps in the same graphic
#' @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 used
#' @param border the colour of each polygons border, default = black
#' @return plots a schematic map of Tasmania, 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)
#'   maptas()
#' }
maptas <- function(view=c(143.5,149,-39,-44.0),defineplot=TRUE,
                     gridon=1.0,maintitle="",infont=7,border=1) {
  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)
  }
  leftlong <- view[1]; rightlong <- view[2]
  uplat <- view[3]; downlat <- view[4]
  plot(c(leftlong,rightlong),c(uplat,downlat),type="n",
       ylim=c(downlat,uplat),yaxs="i",
       xlim=c(leftlong,rightlong),xaxs="i",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")
    abline(h=c(downlat,uplat),lty=1,col=1)
    abline(v=c(leftlong,rightlong),lty=1,col=1)
  }
  ans <- unique(taspoly$poly)
  for (i in ans) {
    pick <- which(taspoly$poly == i)
    polygon(taspoly$Long[pick],taspoly$Lat[pick],border=border)
  }
  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 maptas


#' @title overlap calculates the oid overlaps across year triplets
#'
#' @description overlap takes a set of years and using triplets and the
#'     two and last two years it identifies the oids that are unique to
#'     each year, those that overlap across pairs of years and those that
#'     overlap over all three years.
#'
#' @param adat the data.frame of records selected from an abdat data.frame
#' @param whichvar the variable that identifies the columns od data to use
#' @param yrs which years should be included in the analysis
#'
#' @return a list of records per year, maximum known hexagons, and metadata
#' @export
#'
#' @examples
#' \dontrun{
#'  data(abdat)
#'  addat <- getsubblock(c("13C","13D","13E"),abdat)
#'  addat[,"oid"] <- as.numeric(addat[,"oid"])
#'  whchvar <- "blkg"
#'  yrs <- 2012:2018
#'  answer <- overlap(addat,whichvar=whchvar,yrs=yrs)
#'  answer
#' }
overlap <- function(adat,whichvar="blkg",yrs=c(2012:2018)) {
  #   adat=addat; whichvar="blkg"; yrs=yrs
  label <- paste0(whichvar,yrs)
  pickcol <- grep(whichvar,colnames(adat)) #columns contain blkg
  nyr <- length(yrs)
  columns <- c("Year","Records","Missed")
  records <- matrix(0,nrow=nyr,ncol=3,dimnames=list(yrs,columns))
  maxoid <- dim(adat)[1]
  for (yr in 1:nyr) {  # yr=1
    ydat <- getrecs(pickcol,yr,adat,label)
    pick <- which(ydat[,5] > 0.0)
    ydat <- droplevels(ydat[pick,])
    yrcount <- dim(ydat)[1]
    records[yr,] <- c(yrs[yr],yrcount,(maxoid - yrcount))
  }
  columns <- c("Prev","PrevCurr","PrevNext","Shared","Curr","CurrNext","Next")
  distrib <- matrix(0,nrow=nyr,ncol=length(columns),
                    dimnames=list(yrs,columns))
  # do year 1
  year1 <- getrecs(pickcol,yr=1,adat,label)
  year2 <- getrecs(pickcol,yr=2,adat,label)
  answer <- comptwo(year1,year2)
  distrib[1,] <- c(0,0,0,0,answer["y1"],answer["both"],answer["y2"])
  # do last year
  year1 <- getrecs(pickcol,yr=(nyr-1),adat,label)
  year2 <- getrecs(pickcol,yr=nyr,adat,label)
  answer <- comptwo(year1,year2)
  distrib[nyr,] <- c(answer["y1"],answer["both"],0,0,answer["y2"],0,0)
  # years 2 - (nyr-1)adat
  for (yr in 2:(nyr-1)) { # yr = 2
    yr1 <- getrecs(pickcol,yr-1,adat,label) # minus year
    yr2 <- getrecs(pickcol,yr,adat,label)   # current yr
    yr3 <- getrecs(pickcol,yr+1,adat,label) # plus year
    distrib[yr,] <- compthree(yr1,yr2,yr3)  # c(p,pc,pn,pcn,c,cn,n)
  }
  answer <- cbind(records,distrib)
  metadata <- as.matrix(c("Fishing Year",
                          "Hex in Year = Curr+CurrNext+Shared+PrevCurr",
                          "Known Hex missed in Year",
                          "Hex only in previous year",
                          "Hex in both previous and current years",
                          "Hex in both previous and next years",
                          "Hex in all three years",
                          "Hex only in current year",
                          "Hex in this year and next",
                          "Hex only in next year"),nrow=10,ncol=1)
  rownames(metadata) <- colnames(answer)
  colnames(metadata) <- "Description"
  return(list(answer=answer,totalhex=maxoid, metadata=metadata))
} # end of overlap

#' @title plotgrid plots multiple grid data sets in one graphic
#'
#' @description plotgrid simplifies the plotting of multiple plots of the
#'     grid data
#'
#' @param view the length 4 vector with c(left,right,up, down) that acts
#'     as a window over Tasmania
#' @param dat the data.frame from which to select the 1 hectare grids
#'     within the given view
#' @param plots the par statement regarding the array of plots; the default
#'     = c(2,2)
#' @param label the vector of variable names e.g. c("blkg2012","blkg2013")
#' @param index the index values for which of the label to plot
#'
#' @return nothing, but this does generate a plot of plots in one graphic
#' @export
#'
#' @examples
#' \dontrun{
#'  # ab <- readRDS(paste0(datadir,"G1Ha_2018_10_23.rds"))#get spatial obj
#'  # abdat <- getlatlong(ab)     # convert to data.frame + lat long
#' }
plotgrid <- function(view,dat,plots=c(2,2),label,index=c(1:4)) {
  par(mfcol=plots,mai=c(0.3,0.3,0.05,0.05),oma=c(1.0,1.0,0.0,0.0))
  par(cex=0.85, mgp=c(1.35,0.35,0), font.axis=7,font=7,font.lab=7)
  for (i in index) {
    maptas(view=view,gridon=0.02,defineplot=FALSE,maintitle=label[i])
    mapgrid(view=view,dat,inpch=21,x=label[i],mult=3,col=1,infill=2)
  }
  mtext("Longitude",side=1,outer=TRUE,line=0.0,cex=1.1,font=7)
  mtext("Latitude",side=2,outer=TRUE,line=0.0,cex=1.1,font=7)
} # end of plotgrid

#' @title summarygrid generates some statistics for the variable whchvar
#'
#' @description summarygrid takes the prefix of the variable in the
#'     data.frame from getlatlong and calculates summary statistics.
#'     These include the arithmetic mean, the bias-corrected geometric
#'     mean, the number of NAs, and the number of positive values
#'
#' @param dat the sub-set of the data.frame from getlatlong obtained
#'     from the function getpoints(view,abdat)
#' @param whchvar the prefix of the variable of interest. This can be
#'     "blkg" (or glkg if exploring a greenlip data set), "days","mins",
#'     "divers", "kgday".
#'
#' @return a matrix of summary statistics for each year available and
#'     the total if present
#' @export
#'
#' @examples
#' \dontrun{
#'  # ab <- readRDS(paste0(datadir,"G1Ha_2018_10_23.rds"))#get spatial obj
#'  # abdat <- getlatlong(ab)     # convert to data.frame + lat long
#'  # addat <- getsubblock(c("13C","13D","13E"),abdat)
#'  # summarygrid(addat,"blkg")
#' }
summarygrid <- function(dat,whchvar) { # dat=addat; whchvar=getvar
  pick <- grep(whchvar,colnames(dat))
  avg <- apply(dat[,pick],2,mean,na.rm=TRUE)
  geom <- apply(dat[,pick],2,geomean)
  total <- apply(dat[,pick],2,sum,na.rm=TRUE)/1000
  NumNA <- apply(dat[,pick],2,countNAs)
  NumGT0 <- apply(dat[,pick],2,countgtzero)
  ans <- cbind(avg,geom,total,NumNA,NumGT0)
  return(ans)
} # end of summarygrid