R/optiGrow.R

In geozoning: Zoning Methods for Spatial Data

#############################################################################
#' optiGrow grow an isolated zone by finding a bigger contour line
#'
#' @details Grow an isolated zone by finding a bigger contour line around it. A series of quantiles are tried out. The quantile sequence is generated by genQseq, and depends on LEQ and MAXP parameters. The zone is grown as much as possible, while keeping it isolated from other zones. A zone is isolated if its smaller distance to another zone is smaller than distIsoZ.
#' @param K zoning object (such as returned by calNei function)
#' @param iC index of zone to grow
#' @param qProb probability vector used to generate quantile values
#' @param refPoint xxxx
#' @param map object returned by genMap function
#' @param optiCrit criterion choice
#' @param minSize zone area threshold under which a zone is too small to be manageable
#' @param minSizeNG zone area threshold under which a zone will be removed
#' @param distIsoZ threshold distance to next zone, above which a zone is considered to be isolated
#' @param LEQ length of quantile sequence used to grow zone
#' @param MAXP quantile sequence maximum shift
#' @param simplitol tolerance for spatial polygons geometry simplification
#' @param disp 0: no info, 1: detailed info
#'
#' @return a list with components
#'\describe{
#' \item{crit}{criterion value of the new zoning}
#' \item{area}{area of the grown zone}
#' \item{Zopti}{new zoning geometry (list of SpatialPolygons)}
#' \item{qM}{quantile corresponding to new zone}
#' }
#' @export
#'
#' @examples
#' data(mapTest)
#' qProb=c(0.3,0.5)
#' criti = correctionTree(qProb,mapTest)
#' best = criti$zk[[2]][[1]]
#' Z=best$zonePolygone
#' plotZ(Z)
#' refPoint = rgeos::gCentroid(Z[[4]])
#' sp::plot(refPoint,add=TRUE,col="blue",pch=21)
#' zg=optiGrow(best,4,qProb,refPoint,mapTest) #grow zone 4
#' id=as.numeric(getZoneId(Z[[4]]))
#' linesSp(zg$Zopti[[id]],col="blue") # new zoning with grown zone 4
optiGrow = function(K,iC,qProb,refPoint,map,optiCrit=2,minSize=0.012,minSizeNG=1e-3,distIsoZ=0.075,LEQ=5,MAXP=0.1,simplitol=1e-12,disp=0)
#############################################################################
{
  ## grows current zone by trying out quantiles
   Z=K$zonePolygone
  # quantile sequence is generated by genQseq
   step=map$step
   boundary=map$boundary
   res = NULL
   area = NULL
   iE=detZoneEng(iC,Z,K$zoneNModif)

   if (iE == 0)
   {
   if(disp) print("no englobing zone")
   return(NULL)
   }
  # compute biggest envelope (to keep closest zone at dist >distIso)
  envel=calFrame(iC,Z,K$zoneNModif,distIsoZ)
  if(is.null(envel)) return(NULL)
  #
  # generate quantile values
   Qseq = genQseq(qProb,K,map,iC,iE,LEQ,MAXP,disp)
   critG=rep(0,length(Qseq))
   area=critG
   Zopt=list()

  # try each quantile of the sequence
  # and keep quantile which gives the best criterion
  for (i in 1:length(Qseq))
  {
       resi = findCinZ(iC,Z,K,map,Qseq[i],envel)
       # returns NULL if no grow
       Zopti=NULL
       if(!is.null(resi))
       {
	resp = checkContour(resi$contourSp,step,refPoint,minSizeNG)
	# if  condition not met try next contour
	if (is.null(resp)) next
	Zopti=zoneQ(resi$contourSp,iC,iE,Z,K,simplitol) # current zone is now last zone
	 if (!is.null(Zopti))
      	 {
	 # create comments for holes
	 Zopti = crComment(Zopti)
	 #
	 criti = calcDCrit(Zopti,map,optiCrit,simplitol)
	 if (!is.null(criti))
	   {
	   critG[i]=criti$resCrit
	   Zopt[[i]]=Zopti
	   area[i]=gArea(Zopti[[length(Zopti)]])
	   }
       	 } # end Zopti not null
	} #end biggest contour found for ith quantile

  } #end loop on quantiles

  if (length(Zopt)==0)
     {
     # no zoning found
     if(disp) print("no zoning found")
     return(NULL)
     }
  else # at least one zoning found
     {
       n = rev(order(critG)) # order by criterion value
       mask = area[n]>minSize
       if (any(mask)) # if zone big enough
       {
       nm = n[mask]
       qqm=Qseq[mask]
       mask2=rev(order(area[nm]))
       qqm=qqm[mask2]
       nm = nm[mask2] # sort by area
       iM = nm[1]
       } else # area condition not satisfied - take biggest area even if criterion is not the best one
	{
	 n = rev(order(area))
	 iM = n[1]
	}

  # have to decide compromise between area and criterion value
  return(list(crit=critG,area=area,Zopti=Zopt[[iM]],qM=Qseq[iM]))
  }
}