Join One or Two PolySets using a Logic Operation

Description

Join one or two PolySets using a logic operation.

Usage

joinPolys(polysA,polysB=NULL,operation="INT")

Arguments

Details

This function interfaces with the Clipper library, specifically version 6.2.1 released 2014-10-31 (http://www.angusj.com/delphi/clipper.php), developed by Angus Johnson. Prior to 2013-03-23, joinPolys used the General Polygon Clipper library (http://www.cs.man.ac.uk/aig/staff/alan/software/) by Alan Murta at the University of Manchester. We keep this historic reference to GPC because joinPolys remains faithful to Murta's definition of a generic polygon, which we describe below.

Murta (2004) defines a generic polygon (or polygon set) as zero or more disjoint boundaries of arbitrary configuration. He relates a boundary to a contour, where each may be convex, concave or self-intersecting. In a PolySet, the polygons associated with each unique PID loosely correspond to a generic polygon, as they can represent both inner and outer boundaries. Our use of the term generic polygon includes the restrictions imposed by a PolySet. For example, the polygons for a given PID cannot be arranged arbitrarily.

If polysB is NULL, this function sequentially applies the operation between the generic polygons in polysA. For example, suppose polysA contains three generic polygons (A, B, C). The function outputs the PolySet containing ((A op B) op C).

If polysB is not NULL, this function applies operation between each generic polygon in polysA and each one in polysB. For example, suppose polysA contains two generic polygons (A, B) and polysB contains two generic polygons (C, D). The function's output is the concatenation of A C, B op C, A op D, B op D, with PIDs 1 to 4, respectively. Generally there are n times m comparisons, where n = number of polygons in polysA and m = number of polygons in polysB. If polysB contains only one generic polygon, the function maintains the PIDs from polysA. It also maintains them when polysA contains only one generic polygon and the operation is difference. Otherwise, if polysA contains only one generic polygon, it maintains the PIDs from polysB.

Value

If polysB is NULL, the resulting PolySet contains a single generic polygon (one PID), possibly with several components (SIDs). The function recalculates the PID and SID columns.

If polysB is not NULL, the resulting PolySet contains one or more generic polygons (PIDs), each with possibly several components (SIDs). The function recalculates the SID column, and depending on the input, it may recalculate the PID column.

References

Murta, A. (2004) A General Polygon Clipping Library. Accessed: July 29, 2004.
http://www.cs.man.ac.uk/aig/staff/alan/software/gpc.html

See Also

addPolys, appendPolys, clipPolys, closePolys, fixBound, fixPOS, locatePolys, plotMap, plotPoints, thickenPolys, thinPolys.

Examples

local(envir=.PBSmapEnv,expr={
  oldpar = par(no.readonly=TRUE)
  #--- load the data (if using R)
  if (!is.null(version$language) && (version$language=="R"))
    data(nepacLL,envir=.PBSmapEnv)

  ### Example 1. Cut a triangle out of Vancouver Island
  par(mfrow=c(1,1))
  #--- create a triangle to use in clipping
  polysB <- data.frame(PID=rep(1, 3), POS=1:3,
    X=c(-127.5, -124.5, -125.6), Y = c(49.2, 50.3, 48.6))
  #--- intersect nepacLL with the single polygon, and plot the result
  plotMap(joinPolys(nepacLL, polysB), col="cyan")
  #--- add nepacLL in a different line type to emphasize the intersection
  addPolys(nepacLL, border="purple", lty=3, density=0)
  box()

  ### Example 2. Cut Texada and Lasqueti Islands out of Boxes
  xlim = list(box1=c(-124.8,-124),box2=c(-124,-123.9))
  ylim = list(box1=c(49.4,49.85), box2=c(49.85,49.9))
  Xlim = extendrange(xlim); Ylim=extendrange(ylim)
  polyA = as.PolySet(data.frame(
    PID = rep(1:2,each=4), POS = rep(1:4,2),
    X = as.vector(sapply(xlim,function(x){x[c(1,1,2,2)]})),
    Y = as.vector(sapply(ylim,function(x){x[c(1,2,2,1)]}))
    ), projection="LL")
  data(nepacLLhigh,envir=.PBSmapEnv)
  polyB = nepacLLhigh[is.element(nepacLLhigh$PID,c(736,1912)),]
  polyC = joinPolys(polyA, polyB, "DIFF")
  par(mfrow=c(2,2),cex=1,mgp=c(2,0.5,0))
  plotMap(polyA,col="lightblue",xlim=Xlim,ylim=Ylim)
  addPolys(polyB,col="gold");
  text(mean(Xlim)-0.05,Ylim-0.04,"Boxes (A,B)  and  Isles (C,D)")
  labs = calcCentroid(polyA)
  labs[1,c("X","Y")] = labs[2,c("X","Y")]+c(-0.1,-0.05)
  text(labs[,"X"],labs[,"Y"],c("A","B"),font=2)
  plotMap(polyC[is.element(polyC$PID,1),],col="pink",xlim=Xlim,ylim=Ylim)
  text(mean(Xlim)-0.05,Ylim-0.04,"Box A  \"DIFF\"  Isle C")
  plotMap(polyC[is.element(polyC$PID,3),],col="green",xlim=Xlim,ylim=Ylim)
  text(mean(Xlim)-0.05,Ylim-0.04,"Box A  \"DIFF\"  Isle D")
  plotMap(polyC[is.element(polyC$PID,c(1,3)),],col="cyan",xlim=Xlim,ylim=Ylim)
  text(mean(Xlim)-0.05,Ylim-0.04,"Box A  \"DIFF\"  Isles (C,D)")
  par(oldpar)
})