R/PBSmapping.r

In PBSmapping: Mapping Fisheries Data and Spatial Analysis Tools

##==============================================================================
## Copyright (C) 2003-2022  Fisheries and Oceans Canada
## Nanaimo, British Columbia
## This file is part of PBS Mapping.
##
## PBS Mapping is free software; you can redistribute it and/or modify
## it under the terms of the GNU General Public License as published by
## the Free Software Foundation; either version 2 of the License, or
## (at your option) any later version.
##
## PBS Mapping is distributed in the hope that it will be useful,
## but WITHOUT ANY WARRANTY; without even the implied warranty of
## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
## GNU General Public License for more details.
##
## You should have received a copy of the GNU General Public License
## along with PBS Mapping; if not, write to the Free Software
## Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
##==============================================================================
## PBS Mapping
## File:
##   PBSmapping.r
## Version:
##   See DESCRIPTION.
## Description:
##   Mapping software for R.
## Authors:
##   Specifications: Jon Schnute, Nicholas Boers, Rowan Haigh
##   Code:           Nicholas Boers (and others as documented)
## Change Log:
##    See PBSmapping/ChangeLog
##==============================================================================

##------------------------------------------------
## Functions:
## =========
##  addLabels............Add labels to an exist map
##  addLines.............Add polylines to an existing map plot
##  addPoints............Add points to an existing map plot
##  addPolys.............Add polygons to an existing map plot
##  addStipples..........Add stipples to an existing map plot
##  appendPolys..........Append a two-column matrix to a PolySet
##  as.EventData.........Coerce a data frame to an object with class EventData
##  as.LocationSet.......Coerce a data frame to an object with class LocationSet
##  as.PolyData..........Coerce a data frame to an object with class PolyData
##  as.PolySet...........Coerce a data frame to an object with class PolySet
##  calcArea.............Calculate the areas of polygons found in a PolySet
##  calcCentroid.........Calculate the centroids of polygons found in a PolySet
##  calcLength...........Calculate the length of polylines found in a PolySet
##  calcMidRange.........Calculate the midpoint of the X/Y ranges of polygons found in a PolySet
##  calcSummary..........Apply functions to polygons in a PolySet
##  calcVoronoi..........Calculate the Voronoi (Dirichlet) tesselation for a set of points
##  clipLines............Clip a PolySet, where each unique (PID, SID) describes a polyline
##  clipPolys............Clip a PolySet, where each unique (PID, SID) describes a polygon
##  closePolys...........Close a PolySet of polylines to form polygons
##  combineEvents........Combine measurements of events
##  combinePolys.........Combine several polygons into a single polygon by modifying the PID and SID indices
##  convCP...............Convert output from contourLines into a PolySet
##  convDP...............Convert EventData/PolyData into a PolySet
##  convLP...............Convert two polylines into a polygon
##  convUL...............Convert coordinates between UTM and Lon/Lat
##  dividePolys..........Divide a single polygon into several polygons
##  extractPolyData......Extract PolyData from a PolySet
##  findCells............Find events in grid cells
##  findPolys............Find events in polygons
##  fixBound.............Fix a PolySet's vertices by moving vertices near a boundary to the actual boundary
##  fixPOS...............Fix the POS column of a PolySet by recalculating it using sequential integers
##  importEvents.........Import a text file and convert into EventData
##  importGSHHS..........Import a GSHHG binary file provided by Paul Wessel
##  importLocs...........Import a text file and convert into a LocationSet
##  importPolys..........Import a text file and convert into a PolySet with optional PolyData attribute
##  importShapefile......Import an ArcGIS shapefile (temporarily removed due to maptools disappearing)
##  is.EventData.........Returns TRUE if its argument is of class EventData.
##  is.LocationSet.......Returns TRUE if its argument is of class LocationSet
##  is.PolyData..........Returns TRUE if its argument is of class PolyData
##  is.PolySet...........Returns TRUE if its argument is of class PolySet
##  isConvex.............Determine whether polygons found in a PolySet are convex
##  isIntersecting.......Determine whether polygons found in a PolySet are self-intersecting
##  joinPolys............Join one or two PolySets using a logic operation
##  locateEvents.........Locate events on the current plot (using the locator function)
##  locatePolys..........Locate polygons on the current plot (using the locator function)
##  makeGrid.............Make a grid of polygons, using PIDs and SIDs according to the input arguments
##  makeProps............Append column for polygon property to PolyData based on measurements in PolyData's Z column
##  plotLines............Plot a PolySet as polylines
##  plotMap..............Plot a PolySet as a map, using the correct aspect ratio
##  plotPoints...........Plot EventData/PolyData, where each unique EID or (PID, SID) describes a point
##  plotPolys............Plot a PolySet as polygons
##  print.EventData......Displays information about an EventData object
##  print.LocationSet....Displays information about a LocationSet object
##  print.PolyData.......Displays information about a PolyData object
##  print.PolySet........Displays information about a PolySet object
##  print.summary.PBS....Displays information about a summary.PBS object
##  refocusWorld.........Refocus the 'worldLL'/'worldLLhigh' data sets.
##  summary.EventData....summary method for EventData class
##  summary.LocationSet..summary method for LocationSet class
##  summary.PolyData.....summary method for PolyData class
##  summary.PolySet......summary method for PolySet class
##  thickenPolys.........Thicken a PolySet, where each unique (PID, SID) describes a polygon
##  thinPolys............Thin a PolySet, where each unique (PID, SID) describes a polygon
##------------------------------------------------

PBSprint <- FALSE

## addLabels----------------------------2008-08-25
##  Add labels to an exist map
##  '...' contains arguments for the 'text()' function
##-------------------------------------------NB|RH
addLabels <- function(data, xlim=NULL, ylim=NULL, polyProps=NULL,
   placement="DATA", polys=NULL, rollup=3,
   cex=NULL, col=NULL, font=NULL, ...)
{
	## performs the validation below when we know what type of data we have

	## load the limits, so we can clip before plotting
	if (is.null(xlim))
		xlim <- options()$map.xlim;
	if (is.null(ylim))
		ylim <- options()$map.ylim;
	if (is.null(xlim) || is.null(ylim)) {
		stop("You must create a plot before you can add to one.\n");
	}

	projectionPlot <- options()$map.projection;
	if (placement == "DATA") {
		projectionData <- attr(data, "projection");
	} else {
		projectionData <- attr(polys, "projection");
	}
	.checkProjection(projectionPlot, projectionData);
	zoneData <- attr(data, "zone");

	## assume PolyData
	EventData <- FALSE;
	type <- "p";

	## test for PolyData/EventData
	if (!any(is.element(names(data), c("PID", "EID")))) {
		stop("'data' must be PolyData or EventData.\n");
	}
	## test for 'label' column
	if (!is.element("label", names(data))) {
		stop("'data' is missing a column named 'label'.\n");
	}

	## test for EventData
	if (is.EventData (data)) {
		## ensure placement="DATA"
		if (placement != "DATA") {
			stop("When 'data' is EventData, the 'placement' attribute must equal \"DATA\".\n");
		}

		data <- .validateEventData(data);
		if (is.character(data))
			stop(paste("Invalid EventData 'data'.\n", data, sep=""));
		EventData <- TRUE;
		type <- "e";
	}
	## otherwise, PolyData...
	else {
		data <- .validatePolyData(data);
		if (is.character(data))
			stop(paste("Invalid PolyData 'data'.\n", data, sep=""));
	}

	## if placement is 'DATA', test for 'X'/'Y' columns (they might possibly
	## be missing at this point if 'data' is PolyData)
	if (placement == "DATA") {
		if (any(!is.element(c("X", "Y"), names(data)))) {
			stop("'data' must contain columns 'X' and 'Y' when placement=\"DATA\"'.\n");
		}

	## at this point, placement is DATA and 'data' has label, X, and Y columns
	## fall through for general handling...

	}
	## otherwise, add 'X'/'Y' columns (only possible for PolyData)
	else {
		if (is.null(polys)) {
			stop("Must supply 'polys' when placement is not equal to \"DATA\".\n");
		}
		polys <- .validatePolySet(polys);
		if (is.character(polys))
			stop(paste("Invalid PolySet 'polys'.\n", polys, sep=""));

		## if 'data' contains SID, 'polys' must as well
		if ((is.element("SID", names(data)) && !is.element("SID", names(polys)))) {
			stop("Since 'data' contains an SID field, 'polys' must as well.\n");
		}

		## at this point, placement != DATA, 'data' has label column, and 'polys'
		## is a valid PolySet

		## reduce 'polys' if possible to avoid unnecessary processing (we only
		## need to calculate X/Y for PIDs/SIDs that appear in 'data')
		polys <- polys[is.element(paste(polys$PID, polys$SID), unique(paste(data$PID, data$SID))), ];

		## 'polys': PolySet >> PolyData
		if (placement == "CENTROID") {
			polys <- calcCentroid(polys, rollup=rollup);
		} else if (placement == "MEAN_RANGE") {
			polys <- calcMidRange(polys, rollup=rollup);
		} else if (placement == "MEAN_XY") {
			polys <- calcSummary(polys, rollup=rollup, FUN=mean);
		} else {
			stop(paste("Unknown 'placement'.	Must be one of \"DATA\", \"CENTROID\",",
			"\"MEAN_RANGE\", or \"MEAN_XY\".\n", sep="\n"));
		}

		## at this point, 'polys' is PolyData with X/Y columns; merge label column
		## from 'data' with PolyData 'polys'

		data <- merge(data, polys, by=intersect(intersect(names(data), names(polys)), c("PID", "SID")), all.x=TRUE);

		## check if any entries were not found in PolyData 'polys'
		n <- nrow(data);
		data <- na.omit(data);
		if (nrow(data) != n) {
			warning("Omitted some polygons from 'data' because they did not exist in 'polys'.\n");
		}

		## we can continue as if placement was DATA, as 'data' now has
		## PID, (optional SID), X, Y, and label columns
		## fall through for general handling...
	}

	## clip the data
	data <- data[data$X >= xlim[1] & data$X <= xlim[2] & data$Y >= ylim[1] & data$Y <= ylim[2], ];

	## add the labels
	.addFeature(feature="labels", data=data, polyProps=polyProps, isEventData=EventData, cex=cex, col=col, font=font, ...);

	## prepare the return value
	if (!is.null(data)) {
		## add projection and zone
		attr(data, "projection") <- projectionData;
		attr(data, "zone") <- zoneData;

		## add to the class attribute
		if (EventData && !is.EventData(data, fullValidation=FALSE)) {
			class(data) <- c("EventData", class(data));
		} else if (!EventData && !is.PolyData(data, fullValidation=FALSE)) {
			class(data) <- c("PolyData", class(data));
		}
	}
	invisible (data);
}
##~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~addLabels


## addLines-----------------------------2012-03-01
##  Add polylines to an existing map plot.
##  '...' contains arguments for the 'lines()' function
## ---------------------------------------------NB
addLines <- function(polys, xlim=NULL, ylim=NULL, polyProps=NULL,
   lty=NULL, col=NULL, arrows=FALSE, ...)
{
	## check 'polys'
	polys <- .validatePolySet(polys);
	if (is.character(polys))
		stop(paste("Invalid PolySet 'polys'.\n", polys, sep=""));

	## load/check 'xlim'/'ylim'
	if (is.null(xlim))
		xlim <- options()$map.xlim;
	if (is.null(ylim))
		ylim <- options()$map.ylim;
	if (is.null(xlim) || is.null(ylim)) {
		stop("You must create a plot before you can add to one.\n");
	}

	## load/check 'projection'
	projectionPlot <- options()$map.projection;
	projectionPoly <- attr(polys, "projection");
	.checkProjection(projectionPlot, projectionPoly);

	## clip
	polys <- .clip(polys, xlim, ylim, isPolygons=FALSE, keepExtra=FALSE);
	if (is.null(polys)) {
		warning("All vertices were clipped: no polylines to plot.\n");
		return (NULL);
	}

	## valid columns for PolyProps (described in ?lines under ... and in ?arrows)
	if (arrows)
		pPropsCols <- c("lty", "lwd", "col", "pch", "lend", "ljoin", "lmitre", "angle", "length", "code")
	else
		pPropsCols <- c("lty", "lwd", "col", "pch", "lend", "ljoin", "lmitre", "type")
	## properties that aren't included in par
	nonParProps <- c("angle", "length", "code", "type")
	
	## validate the polyProps argument
	polyProps <- .validatePolyProps(polyProps, parCols=pPropsCols);
	if (is.character(polyProps))
		stop(paste("Invalid PolyData 'polyProps'.\n", polyProps, sep=""));
	if (is.element("SID", names(polyProps))
			&& !is.element("SID", names(polys))) {
		stop("No lines to plot since 'polyProps' contains SIDs but 'polys' does not.\n");
	}

	## obtain default values for 'polyProps': !par values from ?lines/?arrows
	propDefaults <- append(par(setdiff(pPropsCols, nonParProps)),			## in par
		list(type="l",length=0.25,angle=30,code=2)) # !in par
	## don't replace attributes already in 'polyProps' with default values
	propDefaults <- propDefaults[setdiff(names(propDefaults), c(names(polyProps)))]
	## use this function's arguments to override defaults obtained from par
	if (!is.null(col)) propDefaults[["col"]] <- col;
	if (!is.null(lty)) propDefaults[["lty"]] <- lty;
	## look into ..., too, and evaluate statements like
	##	 if (!is.null(list(...)$lwd)) propDefaults[[\"lwd\"]] <- list(...)$lwd;
	pPropsColsTmp <- setdiff(pPropsCols, c("col", "lty"));
	expr <- paste("if (!is.null(list(...)$", pPropsColsTmp, ")) propDefaults[[\"", pPropsColsTmp,
		"\"]] <- list(...)$", pPropsColsTmp, sep="", collapse="; ");
	eval(parse(text=expr))

	## adds an SID column to 'polyProps' if one exists in 'polys'
	polyProps <- .preparePolyProps(polys$PID, polys$SID, polyProps);

	## flesh out 'polyProps' columns (add propDefaults to polyProps, cycling by PID)
	if (length(propDefaults) > 0)
		polyProps <- .addProps(type="p", polyProps=polyProps, propDefaults);
	polyPropsReturn <- polyProps;

	## create an index for the properties
	## use an expression like
	##   polyProps$props <- paste(polyProps$lty, polyProps$lwd, ..., sep="-");
	expr <- paste("polyProps$props <- paste(",
		paste("polyProps$", pPropsCols, sep="", collapse=", "), ", sep=\"-\")", sep="");
	eval(parse(text=expr));
	
	## determine if we can use "fast IDs" (i.e., no 'paste')
	fastIDdig <- .createFastIDdig(polysA=polys, polysB=polyProps, cols=c("PID", "SID"));
	## if fastIDdig > 0, we can use "fast IDs"; must divide SID by 10^fastIDdig

	polyProps$IDs <- .createIDs(polyProps, cols=c("PID", "SID"), fastIDdig);
	polyPropsIdx <- split(polyProps$IDs, polyProps$props);

	## create an 'IDs' column for 'polys'
	polys$IDs <- .createIDs(polys, cols=c("PID", "SID"), fastIDdig);

	## reduce before split
	toSplit <- polyProps[!duplicated(polyProps$props), c(pPropsCols, "props")];
	res <- split(toSplit[, pPropsCols], toSplit$props);
	
	## create lists of X/Y vertices
	polyx <- split(polys$X, polys$IDs);
	polyy <- split(polys$Y, polys$IDs);

	for (prop in names(polyPropsIdx)) {
		toSet <- as.vector(res[prop][[1]]);

		## separate lines with NAs for plotting;
		## for arrows, a start/end of NA does not produce a warning/error
		new.x <- .insertNAs(polyx, polyPropsIdx[[prop]]);
		new.y <- .insertNAs(polyy, polyPropsIdx[[prop]]);

		## when we call lines, set par values like
		##   col=as.vector(toSet$col), lwd=as.vector(toSet$lwd), ...
		if (arrows) {
			eval(parse(text=paste("arrows(",
				"x0=new.x[1:(length(new.x)-1)], ",
				"y0=new.y[1:(length(new.y)-1)], ",
				"x1=new.x[2:length(new.x)], ",
				"y1=new.y[2:length(new.y)], ",
				paste(pPropsCols, "=as.vector(toSet$", pPropsCols, ")", sep="", collapse=", "), ")", sep="")))
		} else {
			eval(parse(text=paste("lines(x=new.x, y=new.y, ",
				paste(pPropsCols, "=as.vector(toSet$", pPropsCols, ")", sep="", collapse=", "), ")", sep="")))
		}
	}

	## add to the class attribute
	if (!is.null(polyPropsReturn) &&
			!is.PolyData(polyPropsReturn, fullValidation=FALSE))
		class(polyPropsReturn) <- c("PolyData", class(polyPropsReturn));

	invisible(polyPropsReturn);
}

## addPoints----------------------------2010-11-10
##	Add points to an existing map plot.
##	'...' contains arguments for the '.addFeature()' function,
##	which will then pass those arguments to 'points' function.
## ---------------------------------------------NB
addPoints <- function(data, xlim=NULL, ylim=NULL, polyProps=NULL,
   cex=NULL, col=NULL, pch=NULL, ...)
{
	## performs the validation below when we know what type of data we have

	if (is.null(xlim))
		xlim <- options()$map.xlim;
	if (is.null(ylim))
		ylim <- options()$map.ylim;
	if (is.null(xlim) || is.null(ylim)) {
		stop("You must create a plot before you can add to one.\n");
	}

	projectionPlot <- options()$map.projection;
	projectionData <- attr(data, "projection");
	.checkProjection(projectionPlot, projectionData);

	## perform some additional validation
	## convert matrix to data frame
	if (is.matrix(data)) {
		data <- .mat2df(data);
	}
	if (!any(is.element(names(data), c("PID", "EID")))) {
		stop("'data' must be either PolyData or EventData.\n");
	}

	## if it contains neither an EID nor PID column, add an EID column
	## and produce a warning
	if (!any(is.element(c("EID", "PID"), names(data)))) {
		data$EID <- 1:nrow(data);
		warning(paste("'data' contains neither an 'EID' nor 'PID' column.	Added an 'EID' column",
			"for plotting only.\n", sep="\n"));
	}

	## is it event data? -- validate it here
	isEventData <- FALSE;
	type <- "p";
	if (is.element("EID", names(data)) && !is.element("PID", names(data))) {
		isEventData <- TRUE;
		type <- "e";
		data <- .validateEventData(data);
		if (is.character(data))
			stop(paste("Invalid EventData 'data'.\n", data, sep=""));
		## we know the EventData has X and Y columns
	} else {
		data <- .validatePolyData(data);
		if (is.character(data))
			stop(paste("Invalid PolyData 'data'.\n", data, sep=""));
		## does the PolyData have X and Y columns?
		if (any(!is.element(c("X", "Y"), names(data))))
			stop("'data' must contain columns 'X' and 'Y'.\n");
	}

	## clip
	data <- data[data$X >= xlim[1] & data$X <= xlim[2]
		& data$Y >= ylim[1] & data$Y <= ylim[2], ];

	polyPropsReturn <-
		.addFeature(feature="points", data=data, polyProps=polyProps,
		isEventData=isEventData, cex=cex, col=col, pch=pch, ...);

	if (!is.null(polyPropsReturn) && !is.PolyData(polyPropsReturn, fullValidation=FALSE))
		class(polyPropsReturn) <- c("PolyData", class(polyPropsReturn));

	invisible (polyPropsReturn);
}
##~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~addPoints


## addPolys-----------------------------2018-06-05
##  Add polygons to an existing map plot.
## ------------------------------------------NB|RH
addPolys <- function(polys, xlim=NULL, ylim=NULL, polyProps=NULL,
   border=NULL, lty=NULL, col=NULL, colHoles=NULL, density=NA, angle=NULL, ...)
{
	## check 'polys'
	polys <- .validatePolySet(polys);
	if (is.character(polys))
		stop(paste("Invalid PolySet 'polys'.\n", polys, sep=""));

	## load/check 'xlim'/'ylim'
	if (is.null(xlim))
		xlim <- options()$map.xlim;
	if (is.null(ylim))
		ylim <- options()$map.ylim;
	if (is.null(xlim) || is.null(ylim)) {
		stop("You must create a plot before you can add to one.\n");
	}

	## load/check 'projection'
	projectionPlot <- options()$map.projection;
	projectionPoly <- attr(polys, "projection");
	.checkProjection(projectionPlot, projectionPoly);

	## clip
	polys <- .clip(polys, xlim, ylim, isPolygons=TRUE, keepExtra=FALSE);
	if (is.null(polys)) {
		warning("All vertices were clipped: no polygons to plot.\n");
		return (NULL);
	}

	## "polyProps": valid columns
	pPropsCols <- c("angle", "border", "col", "colHoles", "lty", "density")

	## validate the polyProps argument
	pProps <- .validatePolyProps(polyProps, parCols=pPropsCols);
	if (is.character(pProps))
		stop(paste("Invalid PolyData 'polyProps'.\n", pProps, sep=""));
	if (is.element("SID", names(pProps)) && !is.element("SID", names(polys))) {
		stop("No polygons to plot since 'polyProps' contains SIDs but 'polys' does not.\n");
	}

	## default values
	pPropsDefs <- list(angle=45, border=1, col=0, colHoles=NULL, lty=par("lty"), density=NA);
	## ... used only when not already in "pProps"
	pPropsDefs <- pPropsDefs[!is.element(names(pPropsDefs), names(pProps))];

	## allow arguments to override defaults
	for (p in setdiff(pPropsCols, "density")) {
		## create an expression like
		##   if (!is.null(angle)) pPropsDefaults[["angle"]] <- angle;
		## for each of the properties and evaluate it
		expr <- paste("if (!is.null(", p, ")) pPropsDefs[[\"", p, "\"]] <- ", p, sep="");
		eval(parse(text=expr));
	}
	## handle "density" with care: in R (1.8.1+), density == NULL is equivalent to
	## density == 0, and density == NA equivalent to density < 0; make the
	## conversion here to match S-PLUS
	if (is.null(density)) {
		pPropsDefs[["density"]] <- 0;
	} else if (is.na(density)) {
		pPropsDefs[["density"]] <- -1;
	} else {
		pPropsDefs[["density"]] <- density;
	}

	## add an SID column to 'polyProps' (if one exists in 'polys')
	pProps <- .preparePolyProps(polys$PID, polys$SID, pProps);

	## add all of the new attributes to 'polyProps'; if "colHoles" not already
	## in "pProps", this function will not add it (since def. value is NULL)
	if (length(pPropsDefs) > 0)
		pProps <- .addProps(type="p", polyProps=pProps, pPropsDefs);

	## only play with "pProps$colHoles" if it's a character vector; otherwise, cannot possibly
	## contain the string "transparent"
	if (is.element("colHoles", names(pProps)) && is.character (pProps$colHoles)) {
		pProps$colHoles[pProps$colHoles == "transparent"] <- NA;
	}

	## save a copy to return
	polyPropsReturn <- pProps;

	## create an index for the properties
	## use an expression like
	##	 pProps$props <- paste(pProps$angle, pProps$border, ..., sep="-");
	## that includes each of the properties and evaluate it
	expr <- paste("pProps$props <- paste(",
		paste("pProps$", pPropsCols, sep="", collapse=", "), ", sep=\"-\")", sep="");
	eval(parse(text=expr));

	## attempt to use fast IDs (i.e., no "paste"); calculate it here so we can
	## reuse it and also use consistent indices between calls to .createIDs
	fastIDdig <- .createFastIDdig(polysA=polys, polysB=pProps, cols=c("PID", "SID"));

	#----------------------------------------------------------------------------
	## Create a data structure that links each plotting style with a list
	## of IDs to plot using that style.
	polys$IDs <- .createIDs(polys, cols=c("PID", "SID"), fastIDdig);
	idxS <- which(!duplicated(polys$IDs))
	idxE <- c((idxS-1)[-1], length(polys$IDs))
	holes <- (polys$POS[idxS] > polys$POS[idxE])
	holesIDs <- (polys[idxS, "IDs"])[holes];

	pProps$IDs <- .createIDs(pProps, cols=c("PID", "SID"), fastIDdig)
	holes <- is.element(pProps$IDs, holesIDs)

	pPropsOuter <- pProps[!holes, c("IDs", "props")]
	pPropsOuter <- split(pPropsOuter$IDs, pPropsOuter$props)

	## the only case where we (completely) do not plot holes as when
	## all of the "colHole"s are trasparent
	if (is.element("colHoles", names(pProps)) && all(is.na(pProps$colHoles))) {
		pPropsHoles <- NA
	} else if (any(holes)) {
		pPropsHoles <- pProps[holes, c("IDs", "props")]
		pPropsHoles <- split(pPropsHoles$IDs, pPropsHoles$props)
	} else {
		pPropsHoles <- NULL
	}

	integrateHoles <- !is.element("colHoles", names(pProps));
	#----------------------------------------------------------------------------

	#----------------------------------------------------------------------------
	## Create a data structure that links each plotting style with the style's
	## associated settings.
	pProps <- pProps[!duplicated(pProps$props), ];
	pProps <- split(pProps[, intersect(pPropsCols, names(pProps))], pProps$props);
	#----------------------------------------------------------------------------

	## create 'polylines' and 'polys'
	polylines <- polys;

	## "integrate holes" for 'polys' if necessary
	## note that "polylines" is not "rolled up"; it contains original PIDs/SIDs
	if (integrateHoles && length(unique(polys$SID)) > 1) {
		polys <- .rollupPolys(polys, rollupMode=2, exteriorCCW=-1, closedPolys=1, addRetrace=TRUE);
	}

	## create an 'IDs' column for both 'polys' and 'polylines'
	polylines$IDs <- .createIDs(polylines, cols=c("PID", "SID"), fastIDdig);
	polys$IDs <- .createIDs(polys, cols=c("PID", "SID"), fastIDdig);

	polylinex <- split(polylines$X, polylines$IDs);
	polyliney <- split(polylines$Y, polylines$IDs);
	polyx <- split(polys$X, polys$IDs);
	polyy <- split(polys$Y, polys$IDs);

	#- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
	## WORK-AROUND:
	## Problem: the R polygon command behaves different from the S polygon
	##  command. The R command draws filling lines according to the 'lty' par
	##  parameter, whereas the S polygon command ignores 'lty' for filling lines.
	#- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
	## Problem:
	##  S-PLUS draws retrace lines when 'lty' != 0.
	##  R refused to draw shading lines (i.e., density) when 'lty' == 0.
	## To stop S-PLUS from plotting retrace lines, the polygon command needs:
	##   1) lty=0 *AND*
	##   2) border=F
	## To stop R from plotting retrace lines, the polygon command needs:
	##   1) lty=0 *OR*
	##   2) border=F
	#- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
	if (!is.null(version$language) && (version$language == "R")) {
		fillBorder <- NA;	 ## in R, NA => omit borders
		polyLty <- 1;
	} else {
		fillBorder <- FALSE;
		polyLty <- 0;
	}

	## walk through the properties, plotting the appropriate fills and edges as we go
	mode <- "o"
	pPropsIdx <- pPropsOuter

	## A) if "pPropsHoles" == NA, all holes were transparent and we have
	##    absolutely nothing to plot for holes; in this instance, range need to include the holes
	## B) if "pPropsHoles" != NA, we either plot borders (if we integrated holes)
	##    or a filled polygon; in these instances, range much include the holes

	## pPropsHoles is NULL when no holes
	## pPropsHoles is NA when "colHoles" are all transparent
	if (!is.null(pPropsHoles) && !any(is.na(pPropsHoles))) {
		range <- c(1:length(pPropsOuter),NA,1:length(pPropsHoles))
	} else {
		range <- c(1:length(pPropsOuter),NA,NA)
	}

	for (propIdx in range) {
		if (is.na(propIdx)) {
			if (mode == "o") {
				## switch modes
				mode <- "h"
				pPropsIdx <- pPropsHoles
				if (is.null(pPropsIdx))
					break
			}
			next
		}

		## set up the drawing parameters
		prop <- names(pPropsIdx)[propIdx]
		toSet <- as.vector(pProps[prop][[1]]);

		## need to as.vector function to remove levels (remove the factors!)
		border <- as.vector(toSet$border);
		lty <- as.vector(toSet$lty);
		if (mode == "o" || is.null(toSet$colHoles)) {
			col <- as.vector(toSet$col);
		} else if (mode == "h") {
			col <- as.vector(toSet$colHoles);
		}
		density <-	as.vector(toSet$density);
		angle <- as.vector(toSet$angle);

		if (!is.na(col)) {
			## plot the fill
			new.polyx <- .insertNAs(polyx, pPropsIdx[[propIdx]])
			new.polyy <- .insertNAs(polyy, pPropsIdx[[propIdx]])

			## Subtle differences exist between the S-PLUS and R "polygon" command.
			## The S-PLUS command fails if not given at least three points, whereas the
			## R command draws nothing for 0-1 point and a line for 2 points.
			## The following tests for 0 points...
			## RH (180605) added 'xpd=TRUE' (see line 1636)
			if (!all(is.na(new.polyx)))
				polygon(x=new.polyx, y=new.polyy, col=col, lty=polyLty,
					border=fillBorder, density=density, angle=angle, xpd=TRUE, ...);

			## plot the edges
			if (is.na(as.logical(border)) || as.logical(border)) {
				if (is.character(border) && !is.na(as.logical(border)))
					border <- as.logical(border);

				new.polylinex <- .insertNAs(polylinex, pPropsIdx[[propIdx]])
				new.polyliney <- .insertNAs(polyliney, pPropsIdx[[propIdx]])

				## RH (180605) added 'xpd=TRUE' (see line 1636)
				polygon(x=new.polylinex, y=new.polyliney, lty=lty,
					border=border, density=0, xpd=TRUE, ...);
			}
		}
	}

	if (!is.null(polyPropsReturn) && !is.PolyData(polyPropsReturn, fullValidation=FALSE))
		class(polyPropsReturn) <- c("PolyData", class(polyPropsReturn));

	invisible(polyPropsReturn);
}
##~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~addPolys


## addStipples--------------------------2004-06-16
##  side=-1 implies outside
##  side= 0 implies both sides
##  side=+1 implies inside
## Note:
##  It is up to the caller to replot the polygons if the stipples cover part of it.
##  'distance' measured as a percentage of the abs(diff(xlim))
## ---------------------------------------------NB
addStipples <- function(polys, xlim=NULL, ylim=NULL, polyProps=NULL,
   side=1, density=1, distance=4, ...)
{
	## initialization
	polys <- .validatePolySet(polys);
	if (is.character(polys))
		stop(paste("Invalid PolySet 'polys'.\n", polys, sep=""));

	if (is.null(xlim))
		xlim <- options()$map.xlim;
	if (is.null(ylim))
		ylim <- options()$map.ylim;
	if (is.null(xlim) || is.null(ylim)) {
		stop("You must create a plot before you can add to one.\n");
	}

	projectionPlot <- options()$map.projection;
	projectionPoly <- attr(polys, "projection");
	.checkProjection(projectionPlot, projectionPoly);

	## always thicken polys to a distance 2% of the plot region (xlim)
	## since thickenPolys() wants a tolerance in kilometers, convert degrees
	## to kilometers if necessary
	if (projectionPlot == "LL") {
		## Equatorial radius 6,378.14 km
		## Polar radius 6,356.78 km
		## Mean radius 6,371.3 km
		## Sources:
		##  http://en.wikipedia.org/wiki/Earth
		##  http://en.wikipedia.org/wiki/Earth_radius
		R <- 6371.3;

		## Use arc length along the line of latitude (mean(ylim))
		##  (NOT a Great Circle distance in this case)

		## This is the radius of the line of latitude, perpedicular to the axis of rotation
		radiusAtLat <- R * sin(pi/2.0 - (abs(mean(ylim))*pi/180.0));
		arcLen <- radiusAtLat * (abs(diff(xlim))*pi/180.0);
		tol <- arcLen * 0.02;
	} else {
		tol <- abs(diff(xlim)) * 0.02;
	}

	## Distance is measured as a percentage abs(diff(xlim)); adjust for latitude
	distance <- distance / 100.0;
	distX <- abs(diff(xlim)) * distance;
	if (projectionPlot == "LL") {
		xyRatio <- cos((mean(ylim) * pi) / 180);
		distY <- distX * xyRatio;
	} else {
		distY <- distX;
	}

	## 1) clip polys to the window to avoid excess computations
	## 2) extend the limits to avoid plotting stipples along the border where
	##    a border only artificially exists because of the clipping; 1.1 fudge-
	##    factor added to ensure points will be off the plot region
	polys <- .clip(polys,
		xlim=xlim + c(-distX, distX)*1.1,
		ylim=ylim + c(-distY, distY)*1.1,
		isPolygons=TRUE, keepExtra=FALSE);
	thick <- thickenPolys(polys, tol=tol, keepOrig=FALSE);
	x <- rep(thick$X, each=(15 * density));
	y <- rep(thick$Y, each=(15 * density));
	events <- data.frame(EID=1:length(x),
		X=x + runif(length(x), -distX, distX),
		Y=y + runif(length(y), -distY, distY));

	## keep only those events that fall within the proper limits
	events <- events[events$X > xlim[1] & events$X < xlim[2]
		& events$Y > ylim[1] & events$Y < ylim[2], ];

	link <- findPolys(events, polys)
	if (side == -1) {
		## if outside
		stip <- events[!is.element(events$EID, unique(link$EID)), ]
	} else if (side == 1) {
		## if inside
		link <- link[link$Bdry == FALSE, ]
		stip <- events[is.element(events$EID, unique(link$EID)), ]
	} else {
		## else both sides
		stip <- events;
	}

	stip <- merge(stip, link[, intersect(names(link), c("EID", "PID", "SID"))], all.x=TRUE, by="EID");
	stip <- stip[, !is.element(names(stip), "EID")];

	.addFeature(feature="points", data=stip, polyProps=polyProps, isEventData=FALSE, ...);

	if (!is.null(polyProps) && !is.PolyData(polyProps, fullValidation=FALSE))
		class(polyProps) <- c("PolyData", class(polyProps));
	invisible(polyProps);
}
##~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~addStipples


##==============================================================================
appendPolys <- function(polys, mat, PID=NULL, SID=NULL, isHole=FALSE)
{
	## validate the arugments
	if (missing(polys) || missing(mat)) {
		stop(paste("Requires arguments 'polys' and 'mat'.	To create a new PolySet, pass",
		"NULL as the 'polys' argument.\n", sep="\n"));
	}
	## validate 'polys'
	if (!is.null(polys)) {
		polys <- .validatePolySet(polys);
		if (is.character(polys))
			stop(paste("Invalid PolySet 'polys'.\n", polys, sep=""));
	}
	## validate 'mat'
	if (!is.matrix(mat) || (ncol(mat) != 2) || !is.numeric(mat)) {
		stop(paste("'mat' must be a matrix with two numeric columns.	The first column contains",
		"the X values and the second contains the Y values.\n", sep="\n"));
	}
	if ((!is.null(PID) && (!is.numeric(PID) || length(PID) != 1)) ||
		(!is.null(SID) && (!is.numeric(SID) || length(SID) != 1)) ||
		(!is.logical(isHole) || length(isHole) != 1)) {
		stop(paste("If supplied as arguments, the length of 'PID', 'SID', and 'isHole' must",
			"equal one.	'PID' and 'SID' must be numerics, and 'isHole' must be a",
			"logical.\n", sep="\n"));
	}

	## save the attributes of the data frame (.validatePolySet returns a data frame);
	## in this case, SAVE ALL ATTRIBUTES since we are "modifying" an existing data set
	attrNames <- setdiff(names(attributes(polys)), c("names", "row.names", "class"));
	attrValues <- attributes(polys)[attrNames];

	## get PID
	if (is.null(PID)) {
		if (is.null(polys)) {
			PID <- 1;
		} else {
			PID <- max(polys$PID) + 1;
		}
	}

	## get SID; only set it if there exists an SID column
	if (is.null(SID) && !is.null(polys) && is.element("SID", names(polys))) {
		curSIDs <- polys[polys$PID == PID, "SID"];
		if (!is.null(curSIDs) && (length(curSIDs) > 0)) {
			SID <- max(curSIDs) + 1;
		} else {
			SID <- 1;
		}
	}

	## build the POS column
	POS <- 1:nrow(mat);
	if (isHole)
		POS <- rev(POS);

	## build the new rows
	if (!is.null(SID)) {
		newRows <- data.frame(PID=rep(PID, nrow(mat)), SID=rep(SID, nrow(mat)), POS=POS, X=mat[, 1], Y=mat[, 2]);
	} else {
		newRows <- data.frame(PID=rep(PID, nrow(mat)), POS=POS, X=mat[, 1], Y=mat[, 2]);
	}

	## if SID not in 'polys', but it is specified here, add it to 'polys'
	if (!is.null(polys) && !is.element("SID", names(polys)) && !is.null(SID))
		polys$SID <- rep(1, nrow(polys));

	## create ordering for 'poly' columns
	extraCols <- vector();
	if (!is.null(polys)) {
		extraCols <- setdiff(names(polys), c("PID", "SID", "POS", "X", "Y"));
		cols <- c("PID", "POS", "X", "Y", extraCols);
	} else {
		cols <- c("PID", "POS", "X", "Y");
	}
	if (is.element("SID", names(polys)) || !is.null(SID))
		cols <- c(cols[1], "SID", cols[-1]);

	## append NA for extraCols in our 'newRows'
	if (length(extraCols) > 0)
		newRows <- cbind(newRows, matrix(NA, nrow=nrow(newRows), ncol=length(extraCols)));

	## bind 'em
	if (!is.null(polys)) {
		ret <- (rbind(polys[, cols], newRows));
	} else {
		ret <- (newRows);
	}

	if (!is.null(ret)) {
		## restore the attributes (including projection and zone)
		attributes(ret) <- c(attributes(ret), attrValues);

		if (!is.PolySet(ret, fullValidation=FALSE))
			class(ret) <- c("PolySet", class(ret));
	}

	return (ret);
}

##==============================================================================
as.EventData <- function(x, projection=NULL, zone=NULL)
{
	## if a data frame, remove all other classes
	if (is.data.frame(x))
		class(x) <- "data.frame"

	x <- .validateEventData(x);
	if (is.character(x))
		stop(paste("Cannot coerce into EventData.\n", x, sep=""));;

	if (!is.EventData(x, fullValidation=FALSE))
		class(x) <- c("EventData", class(x));
	if (!is.null(projection))
		attr(x, "projection") <- projection;
	if (!is.null(zone))
		attr(x, "zone") <- zone;

	return (x);
}

##==============================================================================
as.LocationSet <- function(x)
{
	## if a data frame, remove all other classes
	if (is.data.frame(x))
		class(x) <- "data.frame"

	x <- .validateLocationSet(x);
	if (is.character(x))
		stop(paste("Cannot coerce into a LocationSet.\n", x, sep=""));;

	if (!is.LocationSet(x, fullValidation=FALSE))
		class(x) <- c("LocationSet", class(x));

	return (x);
}

##==============================================================================
as.PolyData <- function(x, projection=NULL, zone=NULL)
{
	## if a data frame, remove all other classes
	if (is.data.frame(x))
		class(x) <- "data.frame"

	x <- .validatePolyData(x);
	if (is.character(x))
		stop(paste("Cannot coerce into PolyData.\n", x, sep=""));;

	if (!is.PolyData(x, fullValidation=FALSE))
		class(x) <- c("PolyData", class(x));
	if (!is.null(projection))
		attr(x, "projection") <- projection;
	if (!is.null(zone))
		attr(x, "zone") <- zone;

	return (x);
}

##==============================================================================
as.PolySet <- function(x, projection=NULL, zone=NULL)
{
	## if a data frame, remove all other classes
	if (is.data.frame(x))
		class(x) <- "data.frame"

	x <- .validatePolySet(x);
	if (is.character(x))
		stop(paste("Cannot coerce into a PolySet.\n", x, sep=""));;

	if (!is.PolySet(x, fullValidation=FALSE))
		class(x) <- c("PolySet", class(x));
	if (!is.null(projection))
		attr(x, "projection") <- projection;
	if (!is.null(zone))
		attr(x, "zone") <- zone;

	return (x);
}

##==============================================================================
calcArea <- function(polys, rollup=3)
{
	polys <- .validatePolySet(polys);
	if (is.character(polys))
		stop(paste("Invalid PolySet 'polys'.\n", polys, sep=""));

	## copied the following code to calcCentroid()
	msg <- ""
	## check if projection is LL
	if (!is.null(attr(polys, "projection"))
		&& !is.na(attr(polys, "projection"))
		&& attr(polys, "projection") == "LL") {
		msg <- paste("To calculate the areas of polygons defined by longitude-latitude",
			"coordinates, this routine first projects them in UTM.\n", sep="\n");
		## check if zone needs to be set
		## SIMILAR code to convUL()
		if (is.null(attr(polys, "zone")) || is.na(attr(polys, "zone"))) {
			m <- mean(polys$X)
			if (m <= -180 || m > 180) {
				msg <- paste(msg, paste(
					"\nAttempted to automatically calculate the missing 'zone' attribute, but",
					"that failed because the mean longitude falls outside the range",
					"-180 < x <= 180.\n", sep="\n"), sep="");
				stop(msg);
			}
			attr(polys, "zone") <- ceiling((m + 180) / 6);
			msg <- paste(msg, paste("\nFor the UTM conversion, automatically detected zone ",
				attr(polys, "zone"), ".\n", sep=""), sep="");
		} else {
			msg <- paste(msg, paste( "\nUsing 'zone' attribute (", 
				attr(polys, "zone"), ").\n", sep=""), sep="");
		}
		warning(msg);
		## now the zone is set
		polys <- convUL(polys);
	}

	## rollup 'polys' to fix CW/CCW and rollup polys
	## since the calcArea function returns signed areas, this is important
	polys <- .rollupPolys(polys, rollupMode=rollup, exteriorCCW=1, closedPolys=1, addRetrace=TRUE);

	inRows <- nrow(polys);
	outCapacity <- inRows;

	## create the data structure that the C functions expect
	if (!is.element("SID", names(polys))) {
		inID <- c(polys$PID, integer(length=inRows), polys$POS);
	} else {
		inID <- c(polys$PID, polys$SID, polys$POS);
	}
	inXY <- c(polys$X, polys$Y);

	## call the C function
	results <- .C("calcArea",
		inID=as.integer(inID),
		inXY=as.double(inXY),
		inVerts=as.integer(inRows),
		outID=integer(2 * outCapacity),
		outArea=double(outCapacity),
		outRows=as.integer(outCapacity),
		outStatus=integer(1),
		PACKAGE="PBSmapping");
	## note: outRows is set to how much space is allocated -- the C function should consider this

	if (results$outStatus == 1) {
		stop("Insufficient physical memory for processing.\n");
	}
	if (results$outStatus == 2) {
		stop(paste("Insufficient memory allocated for output.	Please upgrade to the latest",
			"version of the software, and if that does not fix this problem, please",
			"file a bug report.\n", sep="\n"));
	}

	## determine the number of rows in the result
	outRows <- as.vector(results$outRows);

	## extract the data from the C function results
	if (outRows > 0) {
		d <- data.frame(PID=results$outID[1:outRows],
			SID=results$outID[(outCapacity+1):(outCapacity+outRows)],
			area=results$outArea[1:outRows]);

		if (!is.element("SID", names(polys)) || rollup == 1)
			d$SID <- NULL;

		if (!is.PolyData(d, fullValidation=FALSE))
			class(d) <- c("PolyData", class(d));

		return(d);
	} else {
		return(NULL);
	}
}

##------------------------------------------------------------------------------
## calcCentroid:
##
## Capabilities:
##   [Y] Holes
##   [-] Projection
## Legend:
##   [Y] yes   [N] no   [-] N/A   [?] good question
##
## Returns:
##   PolyData with centroids of all polygons.
##   Values for 'rollup' include:
##     1 - rollup to the PID level
##     2 - rollup to the hole level; i.e., integrate holes
##     3 - do not rollup
##
##   Setting onlyPID=TRUE results in only one centroid per PID.
##------------------------------------------------------------------------------
calcCentroid <- function(polys, rollup=3)
{
	polys <- .validatePolySet(polys);
	if (is.character(polys))
		stop(paste("Invalid PolySet 'polys'.\n", polys, sep=""));

	## to ensure correct orientation of the exterior/interior vertices, run a
	## fixPOS() on the PolySet
	if ((rollup == 1) || (rollup == 2)) {
		polys <- .rollupPolys(polys, rollupMode=rollup, exteriorCCW=1, closedPolys=1, addRetrace=TRUE);
		## As for the overall clockwise/counter-clockwise orientation of the
		## merged PolySet, it shouldn't matter.	The retrace lines don't change
		## the area, and hence shouldn't change the centroid.
	}

	inRows <- nrow(polys);
	outCapacity <- inRows;

	## create the data structure that the C functions expect
	if (!is.element("SID", names(polys))) {
		inID <- c(polys$PID, integer(length=inRows), polys$POS);
	} else {
		inID <- c(polys$PID, polys$SID, polys$POS);
	}
	inXY <- c(polys$X, polys$Y);

	## call the C function
	results <- .C("calcCentroid",
		inID=as.integer(inID),
		inXY=as.double(inXY),
		inVerts=as.integer(inRows),
		outID=integer(2 * outCapacity),
		outXY=double(2 * outCapacity),
		outRows=as.integer(outCapacity),
		outStatus=integer(1),
		PACKAGE="PBSmapping");
	## note: outRows is set to how much space is allocated -- the C function should consider this

	if (results$outStatus == 1) {
		stop("Insufficient physical memory for processing.\n");
	}
	if (results$outStatus == 2) {
		stop(paste("Insufficient memory allocated for output.	Please upgrade to the latest",
			"version of the software, and if that does not fix this problem, please",
			"file a bug report.\n", sep="\n"));
	}

	## determine the number of rows in the result
	outRows <- as.vector(results$outRows);

	## extract the data from the C function results
	if (outRows > 0) {
		d <- data.frame(PID=results$outID[1:outRows],
			SID=results$outID[(outCapacity+1):(outCapacity+outRows)],
			X=results$outXY[1:outRows],
			Y=results$outXY[(outCapacity+1):(outCapacity+outRows)]);

		if (!is.element("SID", names(polys)))
			d$SID <- NULL;

		## add projection and zone
		attr(d, "projection") <- attr(polys, "projection");
		attr(d, "zone") <- attr(polys, "zone");

		if (!is.PolyData(d, fullValidation=FALSE))
			class(d) <- c("PolyData", class(d));

		return(d);
	} else {
		return(NULL);
	}
}

##==============================================================================
# New calcConvexHull using chull() in package grDevices
calcConvexHull <- function(xydata,keepExtra=FALSE)
{
	xydata <- .validateXYData(xydata);
	if (is.character(xydata))
		stop(paste("Invalid X/Y data 'xydata'.\n", xydata, sep=""));

	inRows <- nrow(xydata);
	ii     <- grDevices::chull(xydata$X,xydata$Y); n <- length(ii)
	hull0  <- xydata[ii,]
	zXY    <- is.element(names(hull0),c("X","Y"))
	zExtra <- !zXY
	hull   <- cbind(PID=rep(1,n),POS=1:n,hull0[,zXY])
	if (keepExtra && any(zExtra==TRUE)) {
			hull <- cbind(hull,hull0[,zExtra])
			names(hull)[5:ncol(hull)] <- names(hull0)[zExtra] }

	## add projection and zone attributes
	attr(hull, "inRows") <- inRows;
	attr(hull, "projection") <- attr(xydata, "projection");
	attr(hull, "zone") <- attr(xydata, "zone");

	if (!is.PolySet(hull, fullValidation=FALSE))
		class(hull) <- c("PolySet", class(hull));
	hull <- rbind(hull,hull[1,]); hull[n+1,2] <- n+1; # close polygon
	return(hull)
}

##==============================================================================
calcLength <- function (polys, rollup=3, close=FALSE)
{
	polys <- .validatePolySet(polys);
	if (is.character(polys))
		stop(paste("Invalid PolySet 'polys'.\n", polys, sep=""));

	## according to documentation
	if (rollup == 2) {
		rollup <- 3;
		warning("'rollup=2' has no meaning in 'calcLength'; reset it to 'rollup=3'.\n");
	}

	## STEP 1: compute a distance column
	if (!is.null(attr(polys, "projection")) && !is.na(attr(polys, "projection"))
		&& ((attr(polys, "projection") == "UTM") || (attr(polys, "projection") == "LL")
		|| (attr(polys, "projection") == 1))) {
		## if necessary, close polygons; don't roll them up until after the distance calculation
		if (close) {
			polys <- .rollupPolys(polys, rollupMode=3, exteriorCCW=-1, closedPolys=1, addRetrace=FALSE);
		}
		polys$D <- .calcDist(polys);
	}
	else {
		stop(paste("Invalid projection attribute.	Supported projections include \"LL\",", "\"UTM\", and 1.\n"));
	}

	## STEP 2: prepare the result
	if (is.element("D", names(polys))) {
		## simplify polys
		polys <- polys[, which(is.element(names(polys), c("PID", "SID", "D")))]

		## build index and remove "unnecessary" distances (i.e., distances between polygons)
		polys$IDX <- .createIDs(polys, cols=c("PID", "SID"));
		polys$D[which(!duplicated(polys$IDX)[c(2:nrow(polys), 1)])] <- 0;

		## split and sum
		## if rollup equals 1, split by PID rather than IDX
		if (rollup == 1) {
			idx <- "PID";
		} else {
			idx <- "IDX";
		}
		result <- split(polys$D, polys[[idx]])
		result <- lapply(result, sum);
		polys <- polys[!duplicated(polys[[idx]]), ]
		toMerge <- data.frame(names(result), length=unlist(result))
		names(toMerge)[1] <- idx;
		polys <- merge(polys, toMerge, by=idx)
	}

	if (is.element("SID", names(polys)) && (rollup != 1)) {
		d <- data.frame(PID=polys$PID, SID=polys$SID, length=polys$length);
	} else {
		d <- data.frame(PID=polys$PID, length=polys$length);
	}

	if (!is.null(d) && !is.PolyData(d, fullValidation=FALSE))
		class(d) <- c("PolyData", class(d));

	return (d);
}

##==============================================================================
calcMidRange <- function(polys, rollup=3)
{
	polys <- .validatePolySet(polys);
	if (is.character(polys))
		stop(paste("Invalid PolySet 'polys'.\n", polys, sep=""));

	## calculate range
	polys <- calcSummary (polys, FUN="range", rollup=rollup);

	## calculate mean range
	polys$X <- (polys$X1 + polys$X2) / 2;
	polys$Y <- (polys$Y1 + polys$Y2) / 2;

	## drop obsolete columns
	polys$X1 <- NULL;
	polys$X2 <- NULL;
	polys$Y1 <- NULL;
	polys$Y2 <- NULL;

	## no need to assign 'projection' and 'zone' attributes; 'calcSummary'
	## maintains them

	## assign class if necessary
	if (!is.null(polys) &&
			!is.PolyData(polys, fullValidation=FALSE))
		class(polys) <- c("PolyData", class(polys));

	return (polys);
}

#calcSummary======================================2012-03-01
#	Apply functions to polygons in a PolySet.
#---------------------------------------------------------NB
calcSummary <- function(polys, rollup=3, FUN, ...)
	## Maintains extra attributes.
{
	polys <- .validatePolySet(polys);
	if (is.character(polys))
		stop(paste("Invalid PolySet 'polys'.\n", polys, sep=""));

	## roll up the PolySet if necessary; don't change the number of points
	if (rollup == 1 || rollup == 2) {
		polys <- .rollupPolys(polys, rollupMode=rollup, exteriorCCW=-1,
													closedPolys=-1, addRetrace=FALSE);
	} else if (rollup != 3) {
		stop(
"Invalid value for argument 'rollup'.	Must equal 1, 2 or 3.\n");
	}

	## create an index
	idx <- .createIDs(polys, cols=c("PID", "SID"), fastIDdig=NULL);
	uniqueIdx <- unique(idx);

	## get data value (run FUN on each vector of X and vector of Y)
	dataX <- lapply(split(polys$X, idx), FUN, ...);
	dataY <- lapply(split(polys$Y, idx), FUN, ...);

	nCol <- length(dataX[[1]]);

	## as soon as you bind in a string with numbers in a matrix, the whole
	## matrix becomes characters (all items in matrix have same type)

	## cbind() can result in calling the data frame method if the arguments
	## are all either data frames or vectors, and this will result in the
	## coercion of character vectors to factors. (source: R cbind() help file)

	dataX <- cbind(data.frame(IDX=names(dataX)),
		as.data.frame(matrix(unlist(dataX, use.names=FALSE), ncol=nCol, byrow=TRUE)));
	dataX$IDX <- as.character(dataX$IDX); # remove factor
	dataY <- cbind(data.frame(IDX=names(dataY)),
		as.data.frame(matrix(unlist(dataY, use.names=FALSE), ncol=nCol, byrow=TRUE)));
	dataY$IDX <- as.character(dataY$IDX); # remove factor

	if (nCol > 1) {
		names(dataX) <- c("IDX", paste("X", seq(1, nCol), sep=""));
		names(dataY) <- c("IDX", paste("Y", seq(1, nCol), sep=""));
	} else {
		names(dataX) <- c("IDX", "X");
		names(dataY) <- c("IDX", "Y");
	}

	polys <- polys[!duplicated(idx), ];

	## build PolyData result
	result <- data.frame(IDX=idx[!duplicated(idx)]);
	result$PID <- polys$PID;
	result$SID <- polys$SID;
	result <- merge(result, dataX, by="IDX");
	result <- merge(result, dataY, by="IDX");
	result <- result[, -1];

	if (!is.null(result)) {
		## copy the projection and zone attribute from input
		attr(result, "projection") <- attr(polys, "projection");
		attr(result, "zone") <- attr(polys, "zone");

		if (!is.PolyData(result, fullValidation=FALSE))
			class(result) <- c("PolyData", class(result));
	}

	return (result);
}

#calcVoronoi----------------------------2017-06-27
# Calculate the Voronoi (Dirichlet) tesselation for a set of points.
#-------------------------------------------JTS/RH
calcVoronoi <- function(xydata, xlim=NULL, ylim=NULL, eps=1e-09, frac=0.0001)
{
	.checkRDeps ("calcVoronoi", c("deldir"))

	xydata <- .validateXYData(xydata)
	if (is.character(xydata))
		stop(paste("Invalid X/Y data 'xydata'.\n", xydata, sep=""))
	if (nrow(xydata) < 2)
		stop("This function requires two or more input points.\n")

	if (!is.null(xlim) && (length(xlim) != 2 || diff(xlim) < 0))
		stop("Invalid 'xlim' argument.\n")
	if (!is.null(ylim) && (length(ylim) != 2 || diff(ylim) < 0))
		stop("Invalid 'xlim' argument.\n")

	## if limits null, calculate them by increasing each extent by 10% of its range
	if (is.null(xlim))
		xlim <- range(xydata$X) + (rep(diff(range(xydata$X)) * 0.10, 2) * c(-1, 1))
	if (is.null(ylim))
		ylim <- range(xydata$Y) + (rep(diff(range(xydata$Y)) * 0.10, 2) * c(-1, 1))

	if (!is.null(attr(xydata, "projection"))
			&& attr(xydata, "projection") == "LL")
		warning(paste("When the projection is \"LL\", this function naively treats the data as",
			"a one-to-one projection. Consider converting it to UTM before running", "this function.", sep="\n"));

	dd <- deldir(xydata$X, xydata$Y, rw=c(xlim, ylim), eps=eps, frac=frac, digits=7) ## package deldir
	if (is.null(dd))
		stop(paste("The call to \"deldir\", this function's dependency, failed.	The function",
			"does not support linear (horizontal/vertical) data and that could be the", "cause.", sep="\n"));

	## dd$dirsgs:
	##  (x1, y1) -> (x2, y2): a line segment in a Dirichlet tile
	##  ind1 and ind2: indices of the two points separated by this line segment
	dirsgs <- dd$dirsgs

	## dd$summary:
	##  (x, y): points used in the tesselation; may be less than input if some
	##    points judged to be the same (see "frac" argument)
	##  n.tside: number of side of the Dirichlet tile surrounding the point
	summary <- dd$summary

	## we will use "summary" below AND in ".addCorners" to determine closest polygon for adding corners

	## each line segment belongs to two tiles (ind1 and ind2); create a data
	## structure with two copies of each point (i.e., a copy for each tile)
	pt1tile1 <- dirsgs[, c("x1", "y1", "ind1", "bp1")]
	pt1tile2 <- dirsgs[, c("x1", "y1", "ind2", "bp1")]
	pt2tile1 <- dirsgs[, c("x2", "y2", "ind1", "bp2")]
	pt2tile2 <- dirsgs[, c("x2", "y2", "ind2", "bp2")]
	names(pt1tile1) <- names(pt1tile2) <- names(pt2tile1) <- names(pt2tile2) <- c("X", "Y", "tile", "bp")
	res <- rbind(pt1tile1, pt1tile2, pt2tile1, pt2tile2)
	res <- res[order(res$tile), ]

	## the points in "res" are ordered by "tile" -- the same order of the points in "summary"

	## add the "origin" point (twice) for each tile to the data structure
	res$Xo <- rep(summary[, "x"], times=2*summary[, "n.tside"])
	res$Yo <- rep(summary[, "y"], times=2*summary[, "n.tside"])
	## compute offset from each "origin"
	res$Xdiff <- res$X - res$Xo
	res$Ydiff <- res$Y - res$Yo
	## compute the arctan for each offset, so that we know the order when building polygons
	res$atan2 <- atan2(res$Ydiff, res$Xdiff)

	## NOTE: the order set below is paramount to the next section of code.
	## DO NOT change the ordering unless you revise the next section as well.

	## order them and then eliminate duplicates
	res <- res[order(res$tile, res$atan2), ]
	res <- res[!duplicated(paste(res$tile, res$atan2)), ]

	## number of vertices in each polygon
	nVerts <- diff(c(which(!duplicated(res$tile)), length(res$tile)+1))

	##---------------------------------------------------------------------------
	## At this point, the order of points in each tile is sometimes inappropriate
	## for plotting. To obtain the best results, tiles with boundary points should
	## have one boundary point appear first in the tile and one last.
	##
	## Here is the idea:
	## We have two types of tiles:
	##   (1) "correct" tiles (completely interior OR boundary points in the first and last position) and
	##   (2) "incorrect" tiles (boundary points adjacent to each other in a given tile).
	## We only need to reorder vertices of the "incorrect" tiles. We can think of
	## the reordering as "rotating" the vertices until the first and last are the boundary points.
	##
	## If F is a non-boundary point and T is a boundary point, a tile with seven
	## vertices may initially look like FFFTTFF.	If we rotate these vertices right three times
	##   FFFFTTF (1st rotation)
	##   FFFFFTT (2nd rotation)
	##   TFFFFFT (3rd rotation)
	## we achieve the desired ordering.
	##
	## To calculate the rotation, we look at the index of the first boundary point
	## in the tile.	If we subtract its index from the number of vertices and add
	## one, we obtain the number of right rotations (e.g., 7 - 5 + 1=3 in the above example).
	##
	## Given the number of right rotations, we can build a new index vector that
	## will "rotate" points within the desired tiles.

	## divide 'res' into two data frames, one where tiles touch two or fewer
	## boundary points and one where tiles touch more than two boundary points
	bpByTile <- split(res$bp, res$tile)
	noRotate <- rep((lapply(bpByTile, "sum") > 2), times=lapply(bpByTile, "length"))
	resNoRotate <- res[noRotate, ]   ## touch more than two
	res <- res[!noRotate, ]          ## touch two or fewer

	## vertices in each tile
	tileLen <- diff(c(which(!duplicated(res$tile)), nrow(res)+1))
	## find the tiles to rotate
	bpByTile <- split(res$bp, res$tile)
	bpByTileWhich <- lapply(bpByTile, "which")
	toRotate <- (lapply(bpByTileWhich, "diff") == 1)
	toRotate[is.na(toRotate)] <- FALSE

	## focus on the tiles to fix
	toFix <- bpByTileWhich[toRotate]
	toFixLen <- tileLen[toRotate]
	toFixMin <- lapply(toFix, "min")

	## setup a "rotation" vector
	rotR <- unlist(toFixLen) - unlist(toFixMin)
	rotR <- -1 - rotR

	## a "-1" rotation means no change, "-2" means one to the right, ...
	rot <- rep(-1, length(tileLen))
	rot[toRotate] <- rotR

	## start with an index of 1 .. len for each tile
	newIDX <- unlist(lapply(tileLen, "seq"))

	## rotate those indices
	rot <- rep(rot, times=unlist(tileLen))
	newIDX <- newIDX + rot
	mod <- rep(unlist(tileLen), times=unlist(tileLen))
	newIDX <- (newIDX %% mod) + 1

	## adjust indices from 1 .. len to reflect their position in the data frame
	base <- rep(which(!duplicated(res$tile)) - 1, times=unlist(tileLen))
	newIDX <- newIDX + base

	## reorder
	res <- res[newIDX, ]

	## combine the two 'res' data frames (the no rotate and rotate one)
	if (any(noRotate)) {
		res <- rbind(resNoRotate, res)
		## it is safe to 'order' by tile: according to the R help on 'order',
		## "any unresolved ties will be left in their original ordering."
		res <- res[order(res$tile), ]
	}
	##---------------------------------------------------------------------------

	## build the POS column
	maxPos <- max(nVerts);
	POS <- rep(seq(from=1, to=maxPos), times=length(nVerts));
	n <- vector();
	n[seq(1, by=2, length.out=length(nVerts))] <- nVerts;
	n[seq(2, by=2, length.out=length(nVerts))] <- maxPos - nVerts;
	b <- rep(c(TRUE, FALSE), length.out=length(n));
	POS <- POS[rep(b, n)];

	names(res)[is.element(names(res), "tile")] <- "PID"
	res$POS <- POS

	res <- res[, c("PID", "POS", "X", "Y")]
	class(res) <- c("PolySet", class(res))

	## polygons that straddle corners will be missing the corner points;
	## add the missing corner points now
	res <- .addCorners(res, summary)

	## ensure edges reach the proper extents
	res <- .expandEdges(res, data.frame(X=summary[, "x"], Y=summary[, "y"]), xlim, ylim)

	## set attributes appropriate to the result
	attr(res, "projection") <- 1;
	attr(res, "zone") <- NULL;

	return (res);
}
##~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~calcVoronoi

##==============================================================================
clipLines <- function(polys, xlim, ylim, keepExtra=FALSE)
{
	polys <- .validatePolySet(polys);
	if (is.character(polys))
		stop(paste("Invalid PolySet 'polys'.\n", polys, sep=""));

	ret <- .clip(polys, xlim, ylim, isPolygons=FALSE, keepExtra);

	## .clip retains extra attributes

	if (!is.null(ret) &&
			!is.PolySet(ret, fullValidation=FALSE))
		class(ret) <- c("PolySet", class(ret));

	return (ret);
}

##==============================================================================
clipPolys <- function(polys, xlim, ylim, keepExtra=FALSE)
{
	polys <- .validatePolySet(polys);
	if (is.character(polys))
		stop(paste("Invalid PolySet 'polys'.\n", polys, sep=""));

	ret <- .clip(polys, xlim, ylim, isPolygons=TRUE, keepExtra);

	## .clip retains extra attributes

	if (!is.null(ret) &&
			!is.PolySet(ret, fullValidation=FALSE))
		class(ret) <- c("PolySet", class(ret));

	return (ret);
}

##==============================================================================
closePolys <- function(polys)
{
	polys <- .validatePolySet(polys);
	if (is.character(polys))
		stop(paste("Invalid PolySet 'polys'.\n", polys, sep=""));

	## save attributes, so they can be added to the resulting data frame
	attrNames  <- setdiff(names(attributes(polys)), c("names", "row.names", "class"));
	attrValues <- attributes(polys)[attrNames];

	xlim <- range(polys$X);
	ylim <- range(polys$Y);

	inRows <- nrow(polys);
	outCapacity <- as.integer(2 * inRows);

	## create the data structure that the C function expects
	if (!is.element("SID", names(polys))) {
		inID <- c(polys$PID, integer(length=inRows), polys$POS);
	} else {
		inID <- c(polys$PID, polys$SID, polys$POS);
	}
	inXY <- c(polys$X, polys$Y);
	limits <- c(xlim, ylim);

	## call the C function
	results <- .C("closePolys",
		inID=as.integer(inID),
		inXY=as.double(inXY),
		inVerts=as.integer(inRows),
		limits=as.double(limits),
		outID=integer(3 * outCapacity),
		outXY=double(2 * outCapacity),
		outRows=as.integer(outCapacity),
		outStatus=integer(1),
		PACKAGE="PBSmapping");
	## note: outRows is set to how much space is allocated -- the C function should consider this

	if (results$outStatus == 1) {
		stop("Insufficient physical memory for processing.\n");
	}
	if (results$outStatus == 2) {
		stop(paste("Insufficient memory allocated for output.	Please upgrade to the latest",
			"version of the software, and if that does not fix this problem, please",
			"file a bug report.\n", sep="\n"));
	}

	## determine the number of rows in the result
	outRows <- as.vector(results$outRows);

	## extract the data from the C function results
	if (outRows > 0) {
		d <- data.frame(PID=results$outID[1:outRows],
			SID=results$outID[(outCapacity+1):(outCapacity+outRows)],
			POS=results$outID[(2*outCapacity+1):(2*outCapacity+outRows)],
			X=results$outXY[1:outRows],
			Y=results$outXY[(outCapacity+1):(outCapacity+outRows)]);

		if (!is.element("SID", names(polys)))
			d$SID <- NULL;

		## restore the attributes, including 'projection' and 'zone'
		attributes(d) <- c(attributes(d), attrValues);

		if (!is.PolySet(d, fullValidation=FALSE))
			class(d) <- c("PolySet", class(d));

		return(d);
	} else {
		return(NULL);
	}
}

#combineEvents====================================2012-03-01
#	Combine measurements of events.
#---------------------------------------------------------NB
combineEvents <- function(events, locs, FUN, ..., bdryOK=TRUE)
{
	events <- .validateEventData(events);
	if (is.character(events))
		stop(paste("Invalid EventData 'events'.\n", events, sep=""));
	if (!is.element("Z", names(events))) {
		stop ("EventData is missing required column 'Z'.\n");
	}

	## filter the boundary points
	if (!bdryOK)
		locs <- locs[locs$Bdry == 0, ];

	## make the list...
	if (is.element("SID", names(locs))) {
		colIDs <- 2;
		colNames <- c("PID", "SID", "Z");
		locs <- split(locs$EID, paste(locs$PID, locs$SID, sep="-"));
	} else {
		colIDs <- 1;
		colNames <- c("PID", "Z");
		locs <- split(locs$EID, locs$PID);
	}

	summary <- lapply(locs, function(x, FUN, events, ...){
		FUN(events[is.element(events$EID, x), "Z"], ...);
	}, FUN, events, ...);

	## The as.numeric function calls (below) are very important.
	## They prevent elements in the PID and SID columns from becoming factors.
	output <- as.numeric(unlist(strsplit(names(summary), "-")));
	output <- as.data.frame(matrix(output, byrow=TRUE, ncol=colIDs));

	output <- cbind(output, unlist(summary));
	names(output) <- colNames;

	if (!is.null(output) &&
			!is.PolyData(output, fullValidation=FALSE))
		class(output) <- c("PolyData", class(output));

	return(output);
}

##==============================================================================
combinePolys <- function(polys)
{
	polys <- .validatePolySet(polys)
	if (is.character(polys))
		stop(paste("Invalid PolySet 'polys'.\n", polys, sep=""))

	## save the attributes for the data frame (.validatePolySet returns a data frame)
	attrNames <- setdiff(names(attributes(polys)), c("names", "row.names", "class"));
	attrValues <- attributes(polys)[attrNames];

	IDs <- .createIDs(polys, cols=c("PID", "SID"))

	## create/insert new PID/SID columns
	nVerts <- diff(c(which(!duplicated(IDs)), length(IDs) + 1))
	nPolys <- length(nVerts)
	polys$SID <- rep(1:nPolys, times=nVerts)
	polys$PID <- 1

	## reorder columns
	columns <- c("PID", "SID", "POS", "X", "Y")
	polys <- polys[, c(columns, setdiff(names(polys), columns))]

	## restore the attributes
	attributes(polys) <- c(attributes(polys), attrValues);

	return (polys)
}

##==============================================================================
convCP <- function(data, projection=NULL, zone=NULL)
{
	## do a preliminary test on the 'data' to see if it resulted from a call to
	## contourLines()
	if ((length(data) < 1) || (length(data[[1]]) != 3) || (any(!is.element(names(data[[1]]), c("level", "x", "y"))))) {
		stop(paste("convCP() accepts data from R's contourLines() function.	The format of",
			"'data' does not match the expected format.\n", sep="\n"));
	}

	## 'vData': 'data' converted to a vector
	vData <- unlist(data);

	levelsIdx <- which(names(vData) == "level");
	levelsIdx <- c(levelsIdx, length(vData) + 1);
	## 'nVerts': number of vertices for each polyline
	nVerts <- (diff(levelsIdx) - 1) / 2;
	## 'levelsIdx': indices to each 'level' element in vector 'vData'
	levelsIdx <- levelsIdx[-length(levelsIdx)];

	## create the X column
	## 'n': number of times to repeat each Boolean
	## 'b': vector of Booleans
	n <- vector();
	n[seq(1, by=2, length=length(nVerts))] <- nVerts;
	n[seq(2, by=2, length=length(nVerts))] <- nVerts + 1;
	n <- c(1, n[-length(n)], n[length(n)-1]);
	b <- rep(c(FALSE, TRUE), length.out=length(n));
	x <- vData [rep(b, n)];

	## create the Y column
	n <- vector();
	n[seq(1, by=2, length=length(nVerts))] <- nVerts + 1;
	n[seq(2, by=2, length=length(nVerts))] <- nVerts;
	b <- rep(c(FALSE, TRUE), length.out=length(n));
	y <- vData [rep(b, n)];

	## 'dupLevel': Boolean vector; TRUE when a level is a duplicate
	dupLevel <- duplicated(vData[levelsIdx]);

	## create the PID column
	pid <- vector()
	pid[which(!dupLevel)] <- 1:length(unique(vData[levelsIdx]));
	pid[which(dupLevel)] <- rep(pid[!dupLevel], diff(c(which(!dupLevel), length(levelsIdx)+1))-1);
	PolyData <- data.frame(PID=pid);
	pid <- rep(PolyData$PID, nVerts);

	## create the SID column
	d <- diff(c(which(!dupLevel), length(levelsIdx)+1));
	maxSid <- max(d);
	sid <- rep(seq(from=1, to=maxSid), times=length(d));
	n <- vector();
	n[seq(1, by=2, length.out=length(d))] <- d;
	n[seq(2, by=2, length.out=length(d))] <- maxSid - d;
	b <- rep(c(TRUE, FALSE), length.out=length(n));
	PolyData$SID <- sid[rep(b, n)];
	sid <- rep(PolyData$SID, nVerts);

	## create the POS column
	maxPos <- max(nVerts);
	pos <- rep(seq(from=1, to=maxPos), times=length(nVerts));
	n <- vector();
	n[seq(1, by=2, length.out=length(nVerts))] <- nVerts;
	n[seq(2, by=2, length.out=length(nVerts))] <- maxPos - nVerts;
	b <- rep(c(TRUE, FALSE), length.out=length(n));
	pos <- pos[rep(b, n)];

	## add level column to 'PolyData'
	PolyData$level <- vData[levelsIdx];

	## put the data frames together in a list
	PolySet <- data.frame(PID=pid, SID=sid, POS=pos, X=x, Y=y);

	if (!is.null(PolySet)) {
		## add attributes from arguments; we cannot add them from 'data' since it isn't a PBS Mapping type
		if (!is.null(projection))
			attr(PolySet, "projection") <- projection;
		if (!is.null(zone))
			attr(PolySet, "zone") <- zone;

		## add class attributes as necessary
		if (!is.PolySet(PolySet, fullValidation=FALSE))
			class(PolySet) <- c("PolySet", class(PolySet));
	}
	if (!is.null(PolyData) && !is.PolyData(PolyData, fullValidation=FALSE))
		class(PolyData) <- c("PolyData", class(PolyData));

	return (list(PolySet=PolySet, PolyData=PolyData));
}

##==============================================================================
convDP <- function(data, xColumns, yColumns)
{
	## validate all the arguments
	if (!any(is.element(names(data), c("PID", "EID")))) {
		stop("'data' must be either PolyData or EventData.\n");
	}
	## if 'data' contains both 'PID' and 'EID' columns, assume it is PolyData and ignore 'EID'
	if (is.element("EID", names(data)) && !is.element("PID", names(data))) {
		## replace EID with PID so we can .validatePolyData(), 
		## which is less restrictive than .validateEventData()
		names(data)[which(is.element(names(data), "EID"))] <- "PID";
	}
	data <- .validatePolyData(data);
	if (is.character(data))
		stop(paste("Invalid PolyData 'data'.\n", data, sep=""));

	if (missing(xColumns) || missing(yColumns))
		stop("Must specify 'xColumns' and 'yColumns' vectors.\n");
	if ((length(xColumns) != length(yColumns)) || (length(xColumns)==0))
		stop("Length of 'xColumns' and 'yColumns' must be the same and greater than 0.\n");
	if (!all(is.element(xColumns, names(data))) || !all(is.element(yColumns, names(data))))
		stop(paste("One or more column specified in 'xColumns' or 'yColumns' does not",
			"exist in 'data'.\n", sep="\n"));

	result <- matrix(nrow=(nrow(data) * length(xColumns)), ncol=2);
	for (i in 1:length(xColumns)) {
		result[seq(i, by=length(xColumns), length=nrow(data)), 1] <- data[, xColumns[i]];
		result[seq(i, by=length(yColumns), length=nrow(data)), 2] <- data[, yColumns[i]];
	}
	result <- cbind(rep(1:length(xColumns), nrow(data)), result);
	if (is.element("SID", names(data)))
		result <- cbind(rep(data$SID, each=length(xColumns)), result);
	result <- cbind(rep(data$PID, each=length(xColumns)), result);

	result <- as.data.frame(result);
	if (is.element("SID", names(data))) {
		names(result) <- c("PID", "SID", "POS", "X", "Y");
	} else {
		names(result) <- c("PID", "POS", "X", "Y");
	}

	if (!is.null(result)) {
		## copy the 'projection' and 'zone' attribute from input
		attr(result, "projection") <- attr(data, "projection");
		attr(result, "zone") <- attr(data, "zone");

		if (!is.PolySet(result, fullValidation=FALSE))
			class(result) <- c("PolySet", class(result));
	}

	invisible(result);
}

##==============================================================================
convLP <- function(polyA, polyB, reverse=TRUE)
{
	polyA <- .validatePolySet(polyA);
	if (is.character(polyA))
		stop(paste("Invalid PolySet 'polyA'.\n", polyA, sep=""));
	polyB <- .validatePolySet(polyB);
	if (is.character(polyB))
		stop(paste("Invalid PolySet 'polyB'.\n", polyB, sep=""));

	## test for obvious user errors
	if ((length(unique(polyA$PID)) > 1 || length(unique(polyB$PID)) > 1) ||
		(is.element("SID", names(polyA)) && length(unique(polyA$SID)) > 1) ||
		(is.element("SID", names(polyB)) && length(unique(polyB$SID)) > 1)) {
		warning("Each PolySet should contain only one polyline.\n");
	}

	result <- polyA[, c("X", "Y")];
	if (reverse) {
		result <- rbind(result, polyB[nrow(polyB):1, c("X", "Y")]);
	} else {
		result <- rbind(result, polyB[, c("X", "Y")]);
	}
	result <- cbind(1:nrow(result), result);
	result <- cbind(rep(polyA[1, "PID"], nrow(result)), result);
	names(result) <- c("PID", "POS", "X", "Y");

	if (!is.null(result)) {
		## add 'projection' and 'zone' attributes
		if (is.null(attr(polyA, "projection")) && !is.null(attr(polyB, "projection"))) {
			attr(result, "projection") <- attr(polyB, "projection");
		} else if ((!is.null(attr(polyA, "projection")) && 
			is.null(attr(polyB, "projection"))) || ((!is.null(attr(polyA, "projection")) && 
			!is.null(attr(polyB, "projection"))) && (attr(polyA, "projection") == attr(polyB, "projection")))) {
			attr(result, "projection") <- attr(polyA, "projection");
		}

		if (is.null(attr(polyA, "zone")) && !is.null(attr(polyB, "zone"))) {
			attr(result, "zone") <- attr(polyB, "zone");
		} else if ((!is.null(attr(polyA, "zone")) && is.null(attr(polyB, "zone"))) ||
			((!is.null(attr(polyA, "zone")) && !is.null(attr(polyB, "zone"))) && 
			(attr(polyA, "zone") == attr(polyB, "zone")))) {
			attr(result, "zone") <- attr(polyA, "zone");
		}

		if (!is.PolySet(result, fullValidation=FALSE))
			class(result) <- c("PolySet", class(result));
	}

	invisible (result);
}

##==============================================================================
convUL <- function(xydata, km=TRUE, southern=NULL)
{
	xydata <- .validateXYData(xydata);
	if (is.character(xydata))
		stop(paste("Invalid X/Y data 'xydata'.\n", xydata, sep=""));
	## check for valid projection/zone attributes
	if (!is.element("projection", names(attributes(xydata))) || (!is.element(attr(xydata, "projection"), c("LL", "UTM")))) {
		stop("Missing or invalid projection attribute.\n");
	}

	## automagically set the zone attribute, if possible SIMILAR code to calcArea()
	if ((attr(xydata, "projection") == "LL") && (is.null(attr(xydata, "zone")) || is.na(attr(xydata, "zone")))) {
		m <- mean(xydata$X);
		if ((m <= -180) || (m > 180)) {
			stop(paste("Attempted to automatically calculate the missing 'zone' attribute, but",
				"that failed because the mean longitude falls outside the range",
				"-180 < x <= 180.  Please manually set the attribute and then try again.\n", sep="\n"));
		}
		attr(xydata, "zone") <- ceiling((m + 180) / 6);
		message("convUL: For the UTM conversion, automatically detected zone ", attr(xydata, "zone"), ".");
	}
	else if (!is.element("zone", names(attributes(xydata))) || (attr(xydata, "zone") < 1) || (attr(xydata, "zone") > 60)) {
		stop("Invalid or missing zone attribute; possibly out of valid range.\n");
	}

## TO HERE (formatting code layout only)

	## determine whether in the northern or southern hemisphere
	if (is.null(southern)) {
		if (attr(xydata, "projection") == "LL")
			southern <- (mean(range(xydata$Y)) < 0)
		else if (attr(xydata, "projection") == "UTM")
			southern <- FALSE
		else
			stop("Projection attribute must be UTM or LL.")

		message("convUL: Converting coordinates within the ",
						ifelse(southern, "southern", "northern"), " hemisphere.")
	}

	inXY <- c(xydata$X, xydata$Y);
	outCapacity <- inVerts <- nrow(xydata);
	inProj <- attr(xydata, "projection");
	inZone <- attr(xydata, "zone");
	if (inProj == "UTM") {
		## C function expects data to be in metres...
		if (km) {
			inXY <- inXY * 1000;
		}
		toUTM <- FALSE;
	} else {
		toUTM <- TRUE;
	}

	## call the C function
	results <- .C("convUL",
								inXY=as.double(inXY),
								inVerts=as.integer(inVerts),
								toUTM=as.integer(toUTM),
								zone=as.integer(inZone),
								southern=as.integer(southern),
								outXY=double(2 * outCapacity),
								outVerts=as.integer(outCapacity),
								outStatus=integer(1),
								PACKAGE="PBSmapping");
	## note: outVerts is set to how much space is allocated -- the C function
	##			 should consider this

	if (results$outStatus == 1) {
		stop(
"Insufficient physical memory for processing.\n");
	}
	if (results$outStatus == 2) {
		stop(paste(
"Insufficient memory allocated for output.	Please upgrade to the latest",
"version of the software, and if that does not fix this problem, please",
"file a bug report.\n",
							 sep="\n"));
	}

	## determine the number of rows in the result
	outRows <- as.vector(results$outVerts);

	## extract the data from the C function results
	if (outRows > 0) {
		## C function returns data in metres...
		if (inProj == "LL" && km)
			results$outXY <- results$outXY / 1000;

		xydata$X <- results$outXY[1:outRows];
		xydata$Y <- results$outXY[(outCapacity+1):(outCapacity+outRows)];

		if (inProj == "UTM") {
			attr(xydata, "projection") <- "LL";
		} else {
			attr(xydata, "projection") <- "UTM";
		}

		## no need to set zone attribute, same as input

		## no need to set class, same as input

		return(xydata);
	} else {
		return(NULL);
	}
}

##==============================================================================
dividePolys <- function(polys)
{
	polys <- .validatePolySet(polys)
	if (is.character(polys))
		stop(paste("Invalid PolySet 'polys'.\n", polys, sep=""))

	## save the attributes for the data frame (.validatePolySet returns a data
	## frame)
	attrNames <- setdiff(names(attributes(polys)),
											 c("names", "row.names", "class"));
	attrValues <- attributes(polys)[attrNames];

	IDs <- .createIDs(polys, cols=c("PID", "SID"))
	startsIdx <- which(!duplicated(IDs))

	POSstart <- polys[startsIdx, "POS"]
	## the calculation for POSend works even with one polygon
	POSend <- polys[c(startsIdx[-1] - 1, nrow(polys)), "POS"]
	isHole <- POSstart > POSend
	nVerts <- diff(c(startsIdx, length(IDs) + 1))

	## create new PID/SID columns
	newPID <- rep(1:sum(!isHole),
								times=diff(c(which(!(isHole)), length(isHole)+1)))
	newSID <- (1:length(newPID))
	subSID <- rep(newSID[!duplicated(newPID)] - 1,
								times=diff(c(which(!duplicated(newPID)), length(newPID)+1)))
	newSID <- newSID - subSID

	## insert new PID/SID columns
	polys$PID <- rep(newPID, times=nVerts)
	polys$SID <- rep(newSID, times=nVerts)

	## reorder columns
	columns <- c("PID", "SID", "POS", "X", "Y")
	polys <- polys[, c(columns, setdiff(names(polys), columns))]

	## restore the attributes
	attributes(polys) <- c(attributes(polys), attrValues);

	return (polys)
}

##==============================================================================
extractPolyData <- function(polys)
{
	polys <- .validatePolySet(polys);
	if (is.character(polys))
		stop(paste("Invalid PolySet 'polys'.\n", polys, sep=""));

	## filter out POS, oldPOS, X, Y columns
	polys <- polys[, !is.element(names(polys), c("POS", "oldPOS", "X", "Y"))];
	## check to see if it was reduced to a vector... (or if only PID/SID remained)
	if (length(setdiff(names(polys), c("PID", "SID"))) == 0) {
		stop(
"No extra fields to include in a PolyData object.\n");
	}

	## before building an index, ensure one doesn't already exist
	tempIDX <- date();
	if (is.element(tempIDX, names(polys))) {
		stop(paste(
"Failed attempt to create temporary index column; column already",
"exists.\n",
							 sep="\n"));
	}

	polys[[tempIDX]] <- .createIDs(polys, cols=c("PID", "SID"));
	pdata <- polys[!duplicated(polys[[tempIDX]]),
								 intersect(names(polys), c("PID", "SID", tempIDX))];

	## a function used to process each list element below
	processElement <- function(a) {
		a <- unique(as.vector(a));

		if (length(a) > 1) {
			warning("Added non-unique values to PolyData as 'NA'.\n");
			return (NA);
		} else {
			return (a);
		}
	}

	## process each column individually
	for (column in setdiff(names(polys), c(tempIDX, "PID", "SID"))) {
		temp <- split(polys[[column]], polys[[tempIDX]]);
		temp <- lapply(temp, processElement);
		temp <- data.frame(names(temp), unlist(temp));
		names(temp) <- c(tempIDX, column);
		pdata <- merge(pdata, temp, by=tempIDX, all.x=TRUE);
	}

	## remove index column
	pdata[[tempIDX]] <- NULL;

	## order the results
	if (is.element("SID", names(pdata))) {
		pdata <- pdata[order(pdata$PID, pdata$SID), ];
	} else {
		pdata <- pdata[order(pdata$PID), ];
	}

	## add a class, if necessary
	if (!is.null(pdata) &&
			!is.PolyData(pdata, fullValidation=FALSE))
		class(pdata) <- c("PolyData", class(pdata));

	return (pdata);
}

#findCells------------------------------2014-12-15
# Find events in grid cells
#--------------------------------------------NB|RH
findCells <- function (events, polys, includeBdry=NULL)
{
	events <- .validateEventData(events);
	if (is.character(events))
		stop(paste("Invalid EventData 'events'.\n", events, sep=""));
	polys <- .validatePolySet(polys);
	if (is.character(polys))
		stop(paste("Invalid PolySet 'polys'.\n", polys, sep=""));

	## get parameters to makeGrid: we'll use them to adjust the PID/SID
	## values in the output
	gridPars <- .getGridPars(polys, fullValidation=TRUE);
	if (is.null(gridPars))
		stop("Invalid 'polys'; appears altered since call to 'makeGrid'.\n");

	## prepare data for the C function
	nEvents <- nrow(events);
	brksX <- sort(unique(polys$X));
	brksY <- sort(unique(polys$Y));

	## call the C function for the X's
	resultsX <- .C("findCells",
		inPt=as.double(events$X),
		inPts=as.integer(nEvents),
		inBrk=as.double(brksX),
		inBrks=as.integer(length(brksX)),
		outCell=integer(nEvents),
		outBdry=integer(nEvents),
		outStatus=integer(1),
		PACKAGE="PBSmapping");
	## call the C function for the Y's
	resultsY <- .C("findCells",
		inPt=as.double(events$Y),
		inPts=as.integer(nEvents),
		inBrk=as.double(brksY),
		inBrks=as.integer(length(brksY)),
		outCell=integer(nEvents),
		outBdry=integer(nEvents),
		outStatus=integer(1),
		PACKAGE="PBSmapping");

	## build the output data structures
	d <- data.frame(EID=events$EID,
		PID=resultsX$outCell,
		PIDbdry=resultsX$outBdry,
		SID=resultsY$outCell,
		SIDbdry=resultsY$outBdry);

	bdryX <- d[d$PIDbdry == 1 & d$PID > 1 & d$PID < length(brksX), ];
	bdryX$PID <- bdryX$PID - 1;
	d <- rbind(d, bdryX);
	bdryY <- d[d$SIDbdry == 1 & d$SID > 1 & d$SID < length(brksY), ];
	bdryY$SID <- bdryY$SID - 1;
	d <- rbind(d, bdryY);

	d$Bdry <- as.integer(d$PIDbdry | d$SIDbdry);

	## filter out the bad stuff
	d <- d[d$PID != -1 & d$SID != -1, c("EID", "PID", "SID", "Bdry")];
	if (nrow(d)==0) return(NULL)
	d[d$PID == length(brksX), "PID"] <- length(brksX) - 1;
	d[d$SID == length(brksY), "SID"] <- length(brksY) - 1;

	## at this point:
	## - d uses PID/SID as if called with byrow=T, addSID=T
	## - transform to correct settings...
	if (gridPars$addSID && !gridPars$byrow) {
		## swap
		tmp <- d$PID;
		d$PID <- d$SID;
		d$SID <- tmp;
	} else if (!gridPars$addSID && gridPars$byrow) {
		d$PID <- (d$SID - 1) * (length(brksX) - 1) + d$PID;
		d$SID <- NULL;
	} else if (!gridPars$addSID && !gridPars$byrow) {
		d$PID <- (d$PID - 1) * (length(brksY) - 1) + d$SID;
		d$SID <- NULL;
	}

	## Check new argument `includeBdry' (2014-12-15)
	if (!is.null(includeBdry) && any(d$Bdry==1)) {
		if (includeBdry==0)
			d=d[!is.element(d$Bdry,1),]
		else {
			if (includeBdry==1)
				dd=d[order(d$PID,d$SID),]
			else
				dd=d[rev(order(d$PID,d$SID)),]
			dd=dd[!duplicated(dd$EID),]
			d =dd[order(dd$EID),]
		}
		if (nrow(d)==0) return(NULL)
	}

	if (!is.LocationSet(d, fullValidation=FALSE))
		class(d) <- c("LocationSet", class(d));
	return (d);
}
##~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~findCells


#findPolys------------------------------2014-12-15
# Find events in polygons
#--------------------------------------------NB|RH
findPolys <- function(events, polys, maxRows=1e+05, includeBdry=NULL)
{
	events <- .validateEventData(events);
	if (is.character(events))
		stop(paste("Invalid EventData 'events'.\n", events, sep=""));
	polys <- .validatePolySet(polys);
	if (is.character(polys))
		stop(paste("Invalid PolySet 'polys'.\n", polys, sep=""));

	## create the data structure that the C function expects
	inEvents <- nrow(events);
	inEventsID <- events$EID;
	inEventsXY <- c(events$X, events$Y);

	inPolys <- nrow(polys);
	if (!is.element("SID", names(polys))) {
		inPolysID <- c(polys$PID, integer(length=inPolys), polys$POS);
	} else {
		inPolysID <- c(polys$PID, polys$SID, polys$POS);
	}
	inPolysXY <- c(polys$X, polys$Y);

	## the maximum number of rows in the results occurs when each event
	## lies in all of the polygons; we can easily calculate this case
	## but it often results in too much allocated memory

	## instead, simply use an argument that the user can tune (this makes
	## the behavior more consistent with joinPolys, too)
	outCapacity <- maxRows;

	## call the C function
	results <- .C("findPolys",
		inEventsID=as.integer(inEventsID),
		inEventsXY=as.double(inEventsXY),
		inEvents=as.integer(inEvents),
		inPolysID=as.integer(inPolysID),
		inPolysXY=as.double(inPolysXY),
		inPolys=as.integer(inPolys),
		outID=integer(4 * outCapacity),
		outRows=as.integer(outCapacity),
		outStatus=integer(1),
		PACKAGE="PBSmapping");
	## note: outRows is set to how much space is allocated -- the C function should consider this

	if (results$outStatus == 1) {
		stop(paste(
"Insufficient physical memory for processing.",
"Try reducing this function's 'maxRows' argument to reduce its memory",
"requirements.\n", sep="\n"));
	}
	if (results$outStatus == 2) {
		stop(paste(
"Insufficient memory allocated for output.",
"Try increasing this function's 'maxRows' argument to increase the memory",
"that it allocates for the output.\n", sep="\n"));
	}

	## determine the number of rows in the result
	outRows <- as.vector(results$outRows);
	if (outRows == 0) return(NULL)

	## extract the data from the C function results
	d <- data.frame(EID=results$outID[1:outRows],
		PID=results$outID[(outCapacity+1):(outCapacity+outRows)],
		SID=results$outID[(2*outCapacity+1):(2*outCapacity+outRows)],
		Bdry=results$outID[(3*outCapacity+1):(3*outCapacity+outRows)]);

	## Check new argument `includeBdry' (2014-12-15)
	if (!is.null(includeBdry) && any(d$Bdry==1)) {
		if (includeBdry==0)
			d=d[!is.element(d$Bdry,1),]
		else {
			d0=d[is.element(d$Bdry,0),]
			d1=d[is.element(d$Bdry,1),]
			if (includeBdry==1)
				dd=d1[order(d1$PID,d1$SID),]
			else
				dd=d1[rev(order(d1$PID,d1$SID)),]
			dd=dd[!duplicated(dd$EID),]
			dd=rbind(d0,dd)
			d =dd[order(dd$EID,dd$PID,dd$SID),]
		}
		if (nrow(d)==0) return(NULL)
	}

	if (!is.element("SID", names(polys)))
		d$SID <- NULL;

	if (!is.LocationSet(d, fullValidation=FALSE))
		class(d) <- c("LocationSet", class(d));

	return(d);
}
##~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~findPolys


fixBound <- function(polys, tol)
{
	polys <- .validatePolySet(polys);
	if (is.character(polys))
		stop(paste("Invalid PolySet 'polys'.\n", polys, sep=""));

	if (length(tol) == 1) {
		tolX <- tolY <- tol;
	} else if (length(tol) == 2) {
		tolX <- tol[1]; tolY <- tol[2];
	} else {
		stop(
"Invalid value for 'tol'.	It must be a vector of length 1 or 2.\n");
	}

	xlim <- range(polys$X);
	rangeX <- abs(diff(xlim));
	xDiff <- c(-(tolX * rangeX), (tolX * rangeX));
	leftBnd	<- c(xDiff + xlim[1]);
	rightBnd <- c(xDiff + xlim[2]);

	ylim <- range(polys$Y);
	rangeY <- abs(diff(ylim));
	yDiff <- c(-(tolY * rangeY), (tolY * rangeY));
	lowerBnd <- c(yDiff + ylim[1]);
	upperBnd <- c(yDiff + ylim[2]);

	polys$X[leftBnd[1]	< polys$X & polys$X < leftBnd[2]] <- xlim[1];
	polys$X[rightBnd[1] < polys$X & polys$X < rightBnd[2]] <- xlim[2];
	polys$Y[lowerBnd[1] < polys$Y & polys$Y < lowerBnd[2]] <- ylim[1];
	polys$Y[upperBnd[1] < polys$Y & polys$Y < upperBnd[2]] <- ylim[2];

	## retains attributes including	'projection' and 'zone'

	if (!is.null(polys) &&
			!is.PolySet(polys, fullValidation=FALSE))
		class(polys) <- c("PolySet", class(polys));

	return(polys);
}

##==============================================================================
fixPOS <- function(polys, exteriorCCW=NA)
{
	polys <- .validatePolySet(polys);
	if (is.character(polys))
		stop(paste("Invalid PolySet 'polys'.\n", polys, sep=""));

	if (is.na(exteriorCCW)) {
		exteriorCCW <- -1;
	} else {
		exteriorCCW <- as.logical(exteriorCCW);
		if (is.na(exteriorCCW)) {
			stop("Invalid value for 'exteriorCCW'.	Assign it a Boolean value or NA.\n");
		}
	}

	polys <- .rollupPolys(polys, rollupMode=3, exteriorCCW=exteriorCCW,
												closedPolys=-1, addRetrace=FALSE);

	## .rollupPolys() result retains attributes of 'polys'

	if (!is.null(polys) &&
			!is.PolySet(polys, fullValidation=FALSE))
		class(polys) <- c("PolySet", class(polys));

	return (polys);
}


#=============================================================================
importEvents <- function(EventData, projection=NULL, zone=NULL)
{
	return(as.EventData(read.table(EventData, header=TRUE),
											projection=projection, zone=zone))
}


## importGSHHS--------------------------2023-06-01
## Import a GSHHG binary file provided by Paul Wessel.
## The C call on gshhs using x-limits (0,360) is the only
## limit that ties the polygons near the Greenwhich meridian together.
## The meshing at the International date line is not yet ideal.
## ------------------------------------------NB|RH
importGSHHS <- function(gshhsDB, xlim, ylim, maxLevel=4, n=0, useWest=FALSE)
{
	## Transform all X-coordinates to lie between 0 and 360
	normAngle = function(x){(x + 360.)%%360.}

	## Determine if the specified limits lie in the western hemisphere
	isitWest  = function(lim){
		world = 1:360
		if (diff(lim)>0) span = lim[1]:lim[2]
		else span = setdiff(world,lim[2]:lim[1])
		sum(span>180) > sum(span<=180)
	}
	## Get gshhsDB filename (with and without path information)
	if (missing(gshhsDB))
		gshhsDB <- paste(system.file(package="PBSmapping"), "gshhs_f.b", sep="/")
	else
		gshhsDB <- path.expand(gshhsDB)
	if (!file.exists(gshhsDB))
		stop("Unable to find gshhsDB \"", gshhsDB, "\"")
	gshhsDB.name = basename(substitute(gshhsDB))
	isBorder = grepl("borders",gshhsDB.name)

	##-----FORMAT-BORDERS-CHANGES-------------------
	## Sometimes border longitudes lie between (-180,180) or (0,360).
	## Easiest solution is to import all and use refocusWorld.
	##----------------------------------------------
	if (isBorder) {
		.checkClipLimits(c(xlim, ylim))
		#DBdir  = dirname(substitute(gshhsDB))
		#DBinfo = readLines(paste0(DBdir,"/README.TXT"))
		#DBver  = strsplit(sub("^[[:space:]]+","",DBinfo[grep("Version",DBinfo)[1]]),split=" ")[[1]][2]
		xres <- .Call("importGSHHS", as.character(gshhsDB), as.numeric(c(0,360,-90,90)), 
			as.integer(maxLevel), as.integer(n), PACKAGE = "PBSmapping")
		xPS  = refocusWorld(as.PolySet(as.data.frame(xres),projection="LL"),xlim=xlim,ylim=ylim)
		if (is.null(xPS) || !length(xPS$PID)) return(NULL)
		attr(xPS, "PolyData") <- attr(xres, "PolyData")
		if (useWest) xlim <- xlim - 360.
		if (!all(xlim==attributes(xPS)$rf.xlim)) {  ## refocusWorld adjusts xlims
			xadj = (xlim-attributes(xPS)$rf.xlim)[1]
			xPS$X = xPS$X + xadj
		}
		return(xPS)
	}
	## ----------------------------------------------
	## Headers in GSHHS database range from 0 to 360 while the underlying data
	## range from -180 to 180; if xlim spans the globe use c(0,360) and 
	## refocusWorld the world
	## ----------------------------------------------
	xlim.orig <- NULL
	if (abs(diff(xlim)) >= 360) {
		xlim.orig <- xlim
		xres <- .Call("importGSHHS", as.character(gshhsDB), as.numeric(c(0,360,ylim)), 
			as.integer(maxLevel), as.integer(n), PACKAGE = "PBSmapping")
		xPS  = refocusWorld(as.PolySet(as.data.frame(xres),projection="LL"),xlim=xlim,ylim=ylim)
		if (is.null(xPS) || !length(xPS$PID)) return(NULL)
		PolyData    <- attr(xres, "PolyData")
		clipAsPolys <- attr(xres, "clipAsPolys")
		if (useWest) xlim <- xlim - 360.
		if (!all(xlim==attributes(xPS)$rf.xlim)) {  ## refocusWorld adjusts xlims
			xadj = (xlim-attributes(xPS)$rf.xlim)[1]
			xPS$X = xPS$X + xadj
		}
	}
	else {
		## ----------------------------------------------
		## Try to subset the workd using the C code (tricky)
		## Set up limits: for lon, (-20, 360)
		## ----------------------------------------------
		limits <- list(c(xlim, ylim))
		.checkClipLimits(limits[[1]])
		overlap = ifelse(grepl("gshhs",gshhsDB.name),-20.,0.)
		nxlim   = normAngle(xlim)

		## Check if xlim spans the Greenwich meridian
		## ----------------------------------------------
		green0 = (nxlim[1] >= nxlim[2]) | nxlim[2] > 340.
		green1 = nxlim[1]==0  && grepl("gshhs",gshhsDB.name)
		green2 = nxlim[2]==0  && grepl("gshhs",gshhsDB.name)
		green  = green0 | green1 | green2
		if (green) {
			if (green1) {
				limits[[1]] = c(0.,normAngle(xlim[2]),ylim)
				#limits[[2]] = c(340.,360.,ylim) # maybe not necessary
				bighalf = "R"
			} else if (green2) {
				limits[[1]] = c(normAngle(xlim[1]),360.,ylim)
				limits[[2]] = c(-20.,0.,ylim)
				bighalf = "L"
			} else {
				if ( (360-nxlim[1]) >= nxlim[2] ) {
					limits[[1]] = c(nxlim[1],360.,ylim)
					limits[[2]] = c(overlap,nxlim[2],ylim)
					bighalf = "L"
				}
				else {
					limits[[1]] = c(overlap,xlim[2],ylim)
					limits[[2]] = c(nxlim[1],360.,ylim)
					bighalf = "R"
				}
			}
		}
		isWest = isitWest(limits[[1]][1:2]) # defined solely by user's xlim

		## Initilialize
		x = xbit = list()
		pidlen = 0
		clipAsPolys = logical()
		PolyData = data.frame()

		## Go through 1 or 2 limits, depending on whether `xlim' includes 0 degrees
		## ----------------------------------------------
		for (i in 1:length(limits)) {
			.checkClipLimits(limits[[i]])
			## return an R object -- C call unfortunately converts western hemisphere to negative
			#xres <- .Call("importGSHHS", as.character(gshhsDB), as.numeric(c(0,360,0,90)), as.integer(maxLevel), as.integer(n), PACKAGE = "PBSmapping") # debug only
			xres <- .Call("importGSHHS", as.character(gshhsDB), as.numeric(limits[[i]]),
				as.integer(maxLevel), as.integer(n), PACKAGE = "PBSmapping")
			if (is.null(xres) || !length(xres$PID)) next

			## Documentation states that levels are ignored for lines
			## Check that user-specified levels are honoured (RH 230601)
			if (any(attributes(xres)$PolyData$Level > maxLevel)) {
				pdata = attributes(xres)$PolyData
				keep  = pdata$Level <= maxLevel
				kpid  = pdata$PID[keep]
				pdata = lapply(pdata,function(x){x[keep]})  ## because pdata is a list at this point
				xpid  = is.element(xres[["PID"]], pdata[["PID"]])
				xres  = lapply(xres,function(x){x[xpid]})  ## because xres is a list at this point
				attr(xres, "PolyData") = pdata
			}

			## Attempt to reverse the mess created by the C-code.
			## If xlim spans the Internation Date Line but not Greenwich
			## ----------------------------------------------
			dateline = (normAngle(xlim[1]) >= 45. & normAngle(xlim[2]) <= 315.) & !green
			#if ( any(is.element(c(180,-180),floor(limits[[i]][1]):ceiling(limits[[i]][2]))) & !green ) {
			if (dateline) {
				xneg  = sapply(split(xres[["X"]],xres[["PID"]]),function(x){all(x<0)})
				pneg  = names(xneg[xneg])  ## PIDs with all-negative coordinates
				xmove = is.element(xres[["PID"]],pneg)
				xres[["X"]][xmove] = normAngle(xres[["X"]][xmove])
			}
			#isWestMess = median(xres[["X"]]) < 0  ## `isWestMess' can be different than `isWest'

			## Collect attributes also from each C call
			pdata <- as.data.frame(attr(xres, "PolyData"))
			clipAsPolys <- c(clipAsPolys, attr(xres, "clipAsPolys"))

			## Renumber the PIDs to accommodate a possible second set
			pidnam = unique(xres$PID)
			pid    = (pidlen+1):(pidlen+length(pidnam))
			names(pid) = pidnam
			xres[["oldPID"]] = xres[["PID"]]
			xres[["PID"]] = as.vector(pid[as.character(xres[["PID"]])])
			pdata$oldPID = pdata$PID
			pdata$PID = as.vector(pid[as.character(pdata$PID)])
			pidlen = rev(pid)[1]
			for (j in names(xres)){
				if (length(x)==0) x[[j]] = xres[[j]]
				else              x[[j]] = c(x[[j]],xres[[j]]) # add on results from second limit
			}
			PolyData = rbind(PolyData,pdata)
			xbit[[i]] = xres  ## collect the separate results for debugging only
		} ## end limits loop
#browser();return()

		## Code to activate when debugging
		#if ( "PBSmodelling" %in% rownames(installed.packages()) ) {
		#	if ( exists("XBIT",envir=.PBSmapEnv) )
		#		PBSmodelling::tget(XBIT, tenv=.PBSmapEnv) else XBIT = list()
		#	XBIT[[gshhsDB.name]] = xbit
		#	PBSmodelling::tput(XBIT, tenv=.PBSmapEnv)
		#	PBSmodelling::tput(isWest, tenv=.PBSmapEnv)
		#}
		if (length(x)==0 || !length(x$PID))
			return(NULL)

		## Headers in the GSHHS database range from 0 to 360 while the underlying
		## data ranges from -180 to 180; if our xlim > 180, shift it
		xoff = 0.
		#else if (max(xlim) > 180) { ## too harsh for places like NZ
		if (median(xlim) > 180) {
			xlim.orig <- xlim
			if (useWest) {
				xlim <- xlim - 360.
				xoff = 360.
				if (median(x[["X"]])>0) x[["X"]] = x[["X"]] - 360.
			} else {
				if (median(x[["X"]])<0) x[["X"]] = x[["X"]] + 360.  ## adjust for .Call("importGSHHS")
			}
		}
		XLIM = range(x[["X"]]) # imported data from C-function

		## Convert the list to a PolySet
		xPS <- as.PolySet(as.data.frame(x), projection = "LL")
	} ##----- end else 

	## Clip the PolySet
	## ----------------------------------------------
	if (any(clipAsPolys))
		xPS <- clipPolys(xPS, xlim=xlim, ylim=ylim)
	else
		xPS <- clipLines(xPS, xlim=xlim, ylim=ylim)
	if (is.null(xPS)) return(NULL)
	xPS$oldPOS <- NULL

	## Sort so that holes immediately follow their parents
	xPS <- xPS[order(xPS$PID, xPS$SID), ]

	## Display a message if the output longitudes differ from xlim
	if (!is.null(xlim.orig))
		message("importGSHHS: input xlim was (", paste(xlim.orig, collapse=", "),
					") and the longitude range of the extracted data is (",
					paste(range(xPS$X), collapse=", "), ").")
	attr(xPS, "PolyData") <- PolyData
	return (xPS)
}
##~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~importGSHHS


#=============================================================================
importLocs <- function(LocationSet)
{
	return(as.LocationSet(read.table(LocationSet, header=TRUE)))
}
#=============================================================================

#=============================================================================
importPolys <- function(PolySet, PolyData=NULL, projection=NULL, zone=NULL)
{
	x <- as.PolySet(read.table(PolySet, header=TRUE), projection=projection,
									zone=zone)
	if (!is.null(PolyData)) {
		xDat <- as.PolyData(read.table(PolyData, header=TRUE))
		attr(x, "PolyData") <- xDat
	}
	return(x)
}


### importShapefile----------------------2018-09-06
### This function has several slow parts:
### 1) conversion of the matrix (verts) to X and Y columns in a data frame
### 2) 'lapply's to create POS columns
### Changes 2008-07-15:
###	 Nick's loop to extract data from 'shapeList' has been replaced
###		 by Rowan's series of 'sapply' calls.
###	 Rowan added check for polygons with 0 vertices.
### 2012-04-04: Rowan created function 'placeHoles' 
###	 to place holes under correct solids.
### 2018-09-06: RH modified 'placeHoles' to better deal with orphans,
###	 but can be slow if imported file has many solids|holes|vertices.
### ------------------------------------------NB|RH
#importShapefile <- function (fn, readDBF=TRUE, projection=NULL, zone=NULL, 
#	 minverts=3, placeholes=FALSE, show.progress=FALSE)
#{
#	## initialization
#	.checkRDeps("importShapefile", c("maptools", "foreign"))
#	## call to normalizePath added to perform ~ expansion; otherwise,
#	## pathnames beginning with a ~ fail in the later call to
#	## Rshapeget
#	fn <- normalizePath(fn, mustWork=FALSE)
#	fn <- .getBasename(fn, "shp")
#
#	## test for the required '.shx' file
#	shxFile <- paste(fn, ".shx", sep="")
#	if (!file.exists(shxFile))
#		stop(paste(
#		"Cannot find the index file (\"", shxFile, "\") required to import\n",
#		"the shapefile.\n", sep=""))
#
#	## read shapefile
#	eval(parse(text="shapeList <- .Call(\"Rshapeget\",as.character(fn),as.logical(FALSE),PACKAGE=\"maptools\")"))
#	if (length(shapeList) < 1)
#		stop("The shapefile is empty or an error occurred while importing.\n")
#	shpType=unique(sapply(shapeList,function(x){x$shp.type}))
#	if (length(shpType) != 1)
#		stop ("Supports only a single shape type per shapefile.\n")
#	nVerts=sapply(shapeList,function(x){x$nVerts})
#	v0=is.element(nVerts,0) # any shapefiles with 0 vertices?
#	if (any(v0==TRUE)) {
#		nVerts=nVerts[!v0]; shapeList=shapeList[!v0] }
#	shpID=sapply(shapeList,function(x){x$shpID})
#	nParts=sapply(shapeList,function(x){x$nParts})
#	pStarts=sapply(shapeList,function(x){x$Pstart},simplify=FALSE)
#	if (length(pStarts)!=length(nParts) && all((nParts==sapply(pStarts,length))!=TRUE))
#		stop ("Mismatch in 'nParts' and 'pStarts'.\n")
#	pStarts=unlist(pStarts)
#	v1=unlist(sapply(shapeList,function(x){x$verts[,1]},simplify=FALSE))
#	v2=unlist(sapply(shapeList,function(x){x$verts[,2]},simplify=FALSE))
#	verts=cbind(v1,v2)
#
#	## Keep track of parents and children
#	PC=pStarts
#	zP=is.element(PC,0); PC[zP]=1; PC[!zP]=0
#
#	## reformat results
#	#if (shpType == 3 || shpType == 5) {	## PolySet
#	if (shpType %in% c(3,13,23, 5,15,25)) {	## PolyLine, PolyLineZ, PolyLineM, Polygon, PolygonZ, PolygonM
#		## create preliminary PID/SID columns
#		PID <- rep(1:(length(unique(shpID))), times=nParts)
#		SID <- unlist(lapply(split(nParts, 1:(length(nParts))), "seq"))
#
#		## to determine the number of vertices in each part, we divide the problem
#		## into two cases:
#		## 1) last component/hole of each polygon: the total vertices in the polygon
#		##		less the starting POS of that last component/hole
#		## 2) otherwise: use a "diff" on the starting POS's of each part
#		lastComp <- rev(!duplicated(rev(PID)))
#		nv <- vector()
#		nv[lastComp] <- rep(nVerts, times=nParts)[lastComp] - pStarts[lastComp]
#		nv[!lastComp] <- diff(pStarts)[diff(pStarts) > 0]
#
#		## create PID/SID columns
#		PID <- rep(PID, times=nv)
#		SID <- rep(SID, times=nv)
#		## create POS column; we'll fix the ordering for holes later
#		POS <- unlist(lapply(split(nv, 1:(length(nv))), "seq"))
#		## build the data frame
#		df <- data.frame(PID=PID, SID=SID, POS=POS, X=verts[, 1], Y=verts[, 2])
#
#		#if (shpType == 5) {
#		if (shpType %in% c(5,15,25)) {
#			## PolySet: polygons: reorder vertices for holes
#			or <- .calcOrientation (df)
#			## where "orientation" == -1, we need to reverse the POS
#			or$solid <- is.element(or$orientation,1); or$hole <- !or$solid
#			if (any(or$hole)) {
#				or$nv <- nv
#				toFix <- rep(or$hole, times=or$nv)
#				o <- or$nv[or$hole]
#				newPOS <- unlist(lapply(lapply(split(o, 1:length(o)), "seq"), "rev"))
#				df[toFix, "POS"] <- newPOS	}
#			
#			if (placeholes) {
#				## Fix to the problem where ArcPew does not put solid shapes before holes
#				class(df) <- c("PolySet", setdiff(class(df),"PolySet"))
#				df=placeHoles(df, minVerts=minverts, orient=TRUE, show.progress=show.progress)
#			}
#		}
#		class(df) <- c("PolySet", class(df))
#	#} else if (shpType == 1) {	## EventData
#	} else if (shpType %in% c(1,11,21)) {	## Point, PointZ, PointM
#		EID <- 1:(length(unique(shpID)))
#		df <- data.frame(EID=EID, X=verts[, 1], Y=verts[, 2])
#		class(df) <- c("EventData", class(df))
#	} else {
#		stop ("Shape type not supported.\n");
#	}
#
#	## "cbind" the DBF for EventData or attach as an attribute for PolySets:
#	## According to the "ESRI Shapefile Technical Description", any set of fields
#	## may be present in the DBF file.
#	## The (relevant) requirements are:
#	##	 (1) one record per shape feature (i.e., per PID or EID), and
#	##	 (2) same order as in shape (*.shp) file.
#	dbfFile <- paste(fn, ".dbf", sep="")
#	if (readDBF && !file.exists(dbfFile)) {
#		warning(paste(
#		"The argument 'readDBF' is true but the attribute database\n",
#		"(\"", dbfFile, "\") does not exist.\n", sep=""))
#	} else if (readDBF) {
#		dbf <- foreign::read.dbf (dbfFile)
#		if (shpType == 1) {	## EventData
#			if (nrow(df) != nrow(dbf)) {
#				warning(paste(
#				"The shapefile and its associated DBF do not contain the",
#				"same number of records. DBF ignored.\n", sep="\n"))
#				return (df)
#			}
#			df.class=class(df)
#			df <- cbind(df, dbf)
#			class(df) <- df.class
#		} else if (shpType == 3 || shpType == 5) {
#			## add index to result
#			dbf <- cbind(1:nrow(dbf), dbf)
#			names(dbf)[1] <- "PID"
#			class(dbf) <- c("PolyData", class(dbf))
#			attr(df, "PolyData") <- dbf
#		}
#	## At this point, shpTypes != 1, 3, 5 caused the "stop" above; we do not
#	## need an "else" to check here
#	}
#	attr(df,"parent.child")=PC
#	attr(df,"shpType")=shpType
#	prjFile <- paste(fn, ".prj", sep="")
#	if (file.exists(prjFile)) {
#		prj=scan(prjFile, what="character", quiet=TRUE, skipNul=TRUE)
##browser();return()
#		prj=prj[!is.element(prj,"")][1]
#		if (length(prj)==0 || is.na(prj) || is.null(prj)) prj="Unknown" }
#	else prj="Unknown"
#	attr(df,"prj")=prj
#	xmlFile <- paste(fn, ".shp.xml", sep="")
#	if (file.exists(xmlFile)) {
#		xml=readLines(xmlFile); attr(df,"xml")=xml }
#
#	if (regexpr("GEO",prj)>0 | regexpr("Degree",prj)>0) proj="LL"
#	else if (regexpr("PROJ",prj)>0 && regexpr("UTM",prj)>0) proj="UTM"
#	else proj=1
#	attr(df,"projection")=proj
#
#	if (proj=="UTM" && any(df$X>500))
#		{df$X=df$X/1000; df$Y=df$Y/1000}
#	if (!is.null(zone))
#		attr(df, "zone") <- zone
#	if (!is.null(projection))
#		attr(df,"projection")=projection
##browser();return()
#	return (df) }
###~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~importShapefile


##==============================================================================
is.EventData <- function(x, fullValidation=TRUE)
{
	if (fullValidation) {
		msg <- .validateEventData(x)
		if (is.character(msg))
			return (FALSE)
	}

	return (inherits(x, "EventData", which=TRUE) == 1);
}

##==============================================================================
is.LocationSet <- function(x, fullValidation=TRUE)
{
	if (fullValidation) {
		msg <- .validateLocationSet(x)
		if (is.character(msg))
			return (FALSE)
	}

	return (inherits(x, "LocationSet", which=TRUE) == 1);
}

##==============================================================================
is.PolyData <- function(x, fullValidation=TRUE)
{
	if (fullValidation) {
		msg <- .validatePolyData(x)
		if (is.character(msg))
			return (FALSE)
	}

	return (inherits(x, "PolyData", which=TRUE) == 1);
}

##==============================================================================
is.PolySet <- function(x, fullValidation=TRUE)
{
	if (fullValidation) {
		msg <- .validatePolySet(x)
		if (is.character(msg))
			return (FALSE)
	}

	return (inherits(x, "PolySet", which=TRUE) == 1);
}

##==============================================================================
isConvex <- function(polys)
{
	polys <- .validatePolySet(polys);
	if (is.character(polys))
		stop(paste("Invalid PolySet 'polys'.\n", polys, sep=""));

	## create the data structure that the C function expects
	inPolys <- nrow(polys);
	if (!is.element("SID", names(polys))) {
		inPolysID <- c(polys$PID, vector(length=inPolys));
	} else {
		inPolysID <- c(polys$PID, polys$SID);
	}
	inPolysXY <- c(polys$X, polys$Y);

	## if each point appears in more than two polygons, this number will be
	## too low
	outCapacity <- inPolys;

	## call the C function
	results <- .C("isConvex",
								inID=as.integer(inPolysID),
								inXY=as.double(inPolysXY),
								inVerts=as.integer(inPolys),
								outID=integer(2 * outCapacity),
								outResult=integer(outCapacity),
								outRows=as.integer(outCapacity),
								outStatus=integer(1),
								PACKAGE="PBSmapping");
	## note: outRows is set to how much space is allocated -- the C function
	##			 should consider this

	if (results$outStatus == 1) {
		stop(
"Insufficient physical memory for processing.\n");
	}
	if (results$outStatus == 2) {
		stop(paste(
"Insufficient memory allocated for output.	Please upgrade to the latest",
"version of the software, and if that does not fix this problem, please",
"file a bug report.\n",
							 sep="\n"));
	}

	## determine the number of rows in the result
	outRows <- as.vector(results$outRows);

	## extract the data from the C function results
	if (outRows > 0) {
		d <- data.frame(PID=results$outID[1:outRows],
										SID=results$outID[(outCapacity+1):(outCapacity+outRows)],
										convex=results$outResult[1:outRows]);

		## when I had this conversion in the about data.frame() call, the 'convex'
		## column became factors...
		d$convex <- as.logical(d$convex);

		if (!is.element("SID", names(polys)))
			d$SID <- NULL;

		if (!is.PolyData(d, fullValidation=FALSE))
			class(d) <- c("PolyData", class(d));

		return(d);
	} else {
		return(NULL);
	}
}

##==============================================================================
isIntersecting <- function(polys, numericResult=FALSE)
{
	polys <- .validatePolySet(polys);
	if (is.character(polys))
		stop(paste("Invalid PolySet 'polys'.\n", polys, sep=""));

	## create the data structure that the C function expects
	inPolys <- nrow(polys);
	if (!is.element("SID", names(polys))) {
		inPolysID <- c(polys$PID, integer(length=inPolys), polys$POS);
	} else {
		inPolysID <- c(polys$PID, polys$SID, polys$POS);
	}
	inPolysXY <- c(polys$X, polys$Y);

	## if ecah point appears in more than two polygons, this number will be
	## too low
	outCapacity <- inPolys;

	## call the C function
	results <- .C("isIntersecting",
								inPolysID=as.integer(inPolysID),
								inPolysXY=as.double(inPolysXY),
								inPolys=as.integer(inPolys),
								numericResult=as.integer(numericResult),
								outID=integer(2 * outCapacity),
								outResult=integer(outCapacity),
								outRows=as.integer(outCapacity),
								outStatus=integer(1),
								PACKAGE="PBSmapping");
	## note: outRows is set to how much space is allocated -- the C function
	##			 should consider this

	if (results$outStatus == 1) {
		stop(
"Insufficient physical memory for processing.\n");
	}
	if (results$outStatus == 2) {
		stop(paste(
"Insufficient memory allocated for output.	Please upgrade to the latest",
"version of the software, and if that does not fix this problem, please",
"file a bug report.\n",
							 sep="\n"));
	}

	## determine the number of rows in the result
	outRows <- as.vector(results$outRows);

	## extract the data from the C function results
	if (outRows > 0) {
		d <- data.frame(PID=results$outID[1:outRows],
										SID=results$outID[(outCapacity+1):(outCapacity+outRows)],
										intersecting=results$outResult[1:outRows]);

		if (!numericResult)
			d$intersecting <- as.logical(d$intersecting);

		if (!is.element("SID", names(polys)))
			d$SID <- NULL;

		if (!is.PolyData(d, fullValidation=FALSE))
			class(d) <- c("PolyData", class(d));

		return(d);
	} else {
		return(NULL);
	}
}


#joinPolys------------------------------2014-03-10
# Join one or two PolySets using a logic operation.
#-----------------------------------------------NB
joinPolys <- function(polysA, polysB=NULL, operation="INT")
{
	polysA <- .validatePolySet(polysA)
	if (is.character(polysA))
		stop(paste("Invalid PolySet 'polysA'.\n", polysA, sep=""))

	if (!is.null(polysB)) {
		polysB <- .validatePolySet(polysB)
		if (is.character(polysB))
			stop(paste("Invalid PolySet 'polysB'.\n", polysB, sep=""))
	}

	validOps <- c("INT", "UNION", "DIFF", "XOR")
	if (!is.element(operation, validOps)) {
		stop(paste(
"Invalid \"operation.\"	Must be one of: ", paste(validOps, collapse=", "),
"\n"))
	}
	op <- which(is.element(validOps, operation)) - 1

	inputHasSID <- is.element("SID", names(polysA))
	if (!is.element("SID", names(polysA))) {
		polysA$SID <- 1
	}
	if (!is.element("SID", names(polysB))) {
		polysB$SID <- 1
	}

	## call the C function
	results <- .Call("joinPolys",
									 operation=as.integer(op),
									 sPID=as.integer(polysA$PID),
									 sSID=as.integer(polysA$SID),
									 sPOS=as.integer(polysA$POS),
									 sX	=as.numeric(polysA$X),
									 sY	=as.numeric(polysA$Y),
									 cPID=as.integer(polysB$PID),
									 cSID=as.integer(polysB$SID),
									 cPOS=as.integer(polysB$POS),
									 cX	=as.numeric(polysB$X),
									 cY	=as.numeric(polysB$Y),
									 PACKAGE="PBSmapping")

	if (is.null(results)){
		return(NULL)
	} else if (nrow(results) > 0) {
		## the C routine might add an SID column when one previously didn't exist...
		if (!inputHasSID && all(results$SID == 1)) {
			results$SID <- NULL
		}

		## copy the 'projection' and 'zone' attribute from the input
		attr(results, "projection") <- attr(polysA, "projection")
		attr(results, "zone") <- attr(polysA, "zone")

		if (!is.PolySet(results, fullValidation=FALSE))
			class(results) <- c("PolySet", class(results))

		return(results)
	} else {
		return(NULL)
	}
}
#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~joinPolys

locateEvents <- function(EID, n=512, type="p", ...)
{
	if (!missing(EID)) {
		if (any(!is.numeric(EID))) {
			stop(paste(
"When the 'EID' argument exists, it must contain a vector of numeric",
"values.\n",
								 sep="\n"));
		}
		n <- length(EID);
	}

	events <- locator(n, type, ...);
	## when 0 pts., S-PLUS returns an object of class "missing", whereas
	## R returns NULL; when > 0 pts., both return an object of class "list"

	if (!missing(EID) || (data.class(events) == "list")) {
		if (missing(EID)) {
			EventData <- data.frame(EID=1:length(events$x),
															X=events$x,
															Y=events$y);
		} else {
			EventData <- data.frame(EID=EID,
															X=c(events$x, rep(NA, n - length(events$x))),
															Y=c(events$y, rep(NA, n - length(events$y))));
		}

		if (!is.null(EventData)) {
			## set the projection attribute
			attr(EventData, "projection") <- options()$map.projection;
			## does not set "zone" attribute as per documentation

			if (!is.EventData(EventData, fullValidation=FALSE))
				class(EventData) <- c("EventData", class(EventData));
		}

		return (EventData);
	} else {
		return (NULL);
	}
}

##==============================================================================
locatePolys <- function(pdata, n=512, type="o", ...)
{
	if (!missing(pdata)) {
		pdata <- .validatePolyData(pdata);
		if (is.character(pdata))
			stop(paste("Invalid PolyData 'pdata'.\n", pdata, sep=""));
		polys <- nrow(pdata);
	} else {
		polys <- 1;
	}

	PID <- 1;
	SID <- NULL;
	output <- NULL;
	for (i in 1:polys) {
		## replace as necessary from 'pdata'
		if (!missing(pdata)) {
			PID <- pdata[i, "PID"];
			if (is.element("SID", names(pdata)))
				SID <- pdata[i, "SID"];
			if (is.element("n", names(pdata)))
				n <- pdata[i, "n"];
			if (is.element("type", names(pdata)))
				type <- pdata[i, "type"]
		}

		pts <- locator(n=abs(n), type=type, ...);
		if (length(pts$x) == 0) {
			stop(
"You must locate at least one point.\n");
		}
		if (type == "l" || type == "o")
			lines(pts$x[c(1, length(pts$x))], pts$y[c(1, length(pts$y))], ...);

		if (missing(pdata) || (!missing(pdata) && !is.element("n", names(pdata))))
			## keep the sign in case we are generating a hole with no parent
			n <- sign(n) * length(pts$x);

		## create vectors for data frame
		PID <- rep(PID, abs(n));
		if (!is.null(SID))
			SID <- rep(SID, abs(n));
		if (n > 0) POS <- 1:n else POS <- abs(n):1;
		if (!missing(pdata) && is.element("n", names(pdata))) {
			pts$x <- c(pts$x, rep(NA, abs(n) - length(pts$x)));
			pts$y <- c(pts$y, rep(NA, abs(n) - length(pts$y)));
		}

		if (!is.null(SID)) {
			output <- rbind(output, data.frame(PID=PID, SID=SID, POS=POS,
																				 X=pts$x, Y=pts$y));
		} else {
			output <- rbind(output, data.frame(PID=PID, POS=POS,
																				 X=pts$x, Y=pts$y));
		}
	}

	if (!is.null(output)) {
		## set the projection attribute
		attr(output, "projection") <- options()$map.projection;
		## does not set "zone" attribute as per documentation

		if (!is.PolySet(output, fullValidation=FALSE))
			class(output) <- c("PolySet", class(output));
	}

	return(output)
}


## makeGrid-----------------------------2019-01-04
## Make a regular tesselation (default squares)
## The code checks the Xs and Ys for errors, 
## and then generates the Xs and Ys.
## Hexagonal tesselation added (RH 181120)
## Changed hexagon orientation to match rectangles when byrow=T (RH 190104)
## ------------------------------------------NB|RH
makeGrid <- function(x, y, byrow=TRUE, addSID=TRUE, 
	 projection=NULL, zone=NULL, type="rectangle")
{
	## check Xs and Ys for errors
	dX <- diff(x);
	dY <- diff(y);
	if (any(c(dX, dY) < 0)) {
		stop("Either the 'x' vector or the 'y' vector was not sequential.\n");
	}
	if (length(x) < 2 || length(y) < 2) {
		stop("Both x and y must have a length > 1.\n");
	}
	numX <- length(x);
	numY <- length(y);

	if (type=="rectangle") {
		## generate the Xs and Ys
		nV <- 4; nX=numX - 1; nY=numY -1
		x1 <- rep(x[1:(numX - 1)], times=numY - 1);
		y1 <- rep(y[1:(numY - 1)], each=numX - 1);
		x2 <- x1 + rep(dX, times=numY - 1);
		y2 <- y1 + rep(dY, each=numX - 1);

		xP <- as.vector(rbind(x1, x2, x2, x1))
		yP <- as.vector(rbind(y1, y1, y2, y2))

		## set up for (byrow && addSID)
		pid <- rep(rep(1:nX, each=nV), times=nY);
		sid <- rep(1:nY, each=nV * nX);
		pos <- rep(1:nV, times=(nX * nY));

	} else if (grepl("^hexagon$",type)) { ## added by RH 181120

		nV	 =6; nX=numX; nY=numY
		rX	 =dX[1]/2
		rY	 =dY[1]/2
		ff	 =0.5 ## intrusion/extrusion factor for hexagon vertices

		## Change hexagon behaviour byrow to be consistent with rectangles -- byrow increments PID by column (RH 190104)
		if (byrow) {
			## hexagons (flat-topped) by contiguous columns -- odd-even column centroids need to be staggered
			## Start vertices from left-most
			oY	 =rY * c(0, 1, 1, 0, -1, -1)
			yodd =as.vector(sapply(y,function(a,b){a+b},b=oY))
			yeven=as.vector(sapply(y[-1],function(a,b){a+b},b=oY)) - rY
			yP	 =rep(c(yodd,yeven),floor(nX/2))
			if (!(nX%%2)==0)
				yP=c(yP, yodd)

			## Get PID, SID, POS before xP
			oddY=!(1:nX)%%2==0
			nny =rep(c(nY,nY-1),length(oddY))[1:length(oddY)]
			pid =rep(1:length(nny), nny*nV)
			sid =as.vector(unlist(lapply(table(pid)/nV,function(a,n){rep(1:a,each=n)},n=nV)))
			pos =rep(1:nV, times=sum(table(pid)/nV))

			## Deal with X based on PID
			oX	 =rX * c(-(1+ff), -ff, ff, (1+ff), ff, -ff)
			xlst=lapply(1:nX, function(i,xx,nn,oo) {
				rep(xx[i] + oo, nn[i])
			}, xx=x, nn=as.vector(table(pid)/nV), oo=oX)
			xP=as.vector(unlist(xlst))
		} else {
			## hexagons (pointy-topped) by contiguous rows -- odd-even row centroids need to be staggered
			## Start vertices from bottom-most
			oX	 =rX * c(0,-1,-1,0,1,1)
			#oX	 =rX * sin(arads)
			xodd =as.vector(sapply(x,function(a,b){a+b},b=oX))
			xeven=as.vector(sapply(x[-1],function(a,b){a+b},b=oX)) - rX
			xP	 =rep(c(xodd,xeven),floor(nY/2))
			if (!(nY%%2)==0)
				xP=c(xP, xodd)

			## Get PID, SID, POS before yP
			oddX=!(1:nY)%%2==0
			nnx =rep(c(nX,nX-1),length(oddX))[1:length(oddX)]
			pid =rep(1:length(nnx), nnx*nV)
			sid =as.vector(unlist(lapply(table(pid)/nV,function(a,n){rep(1:a,each=n)},n=nV)))
			pos =rep(1:nV, times=sum(table(pid)/nV))

			## Deal with Y based on PID
			oY	 =rY * c(-(1+ff), -ff, ff, (1+ff), ff, -ff)
			#oY	 =rY * cos(arads) * aspect
			ylst=lapply(1:nY, function(i,yy,nn,oo) {
				rep(yy[i] + oo, nn[i])
			}, yy=y, nn=as.vector(table(pid)/nV), oo=oY)
			yP=as.vector(unlist(ylst))
		}
	}

	## Make data frame
	pSet <- data.frame(PID=pid, SID=sid, POS=pos, X=xP, Y=yP);
	if (!addSID || (!byrow && type=="rectangle")) {
		## adjust for alternative indexing
		pSet <- .createGridIDs(pSet, addSID, byrow);
	}

	if (!is.null(pSet)) {
		if (!is.null(projection))
			attr(pSet, "projection") <- projection;
		if (!is.null(zone))
			attr(pSet, "zone") <- zone;

		if (!is.PolySet(pSet, fullValidation=FALSE))
			class(pSet) <- c("PolySet", class(pSet));
	}

	return(pSet);
}
##~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~makeGrid


##==============================================================================
# the default value for 'propVals' is less than optimal if 'breaks' is a vector
# of length 1...
makeProps <- function(pdata, breaks, propName="col",
											propVals=1:(length(breaks) - 1))
{
	pdata <- .validatePolyData(pdata);
	if (is.character(pdata))
		stop(paste("Invalid PolyData 'pdata'.\n", pdata, sep=""));

	if (!is.element("Z", names(pdata)))
		stop(
"'Z' column is missing in 'pdata'.\n");

	if ((length(breaks) == 1 && breaks < 2) || (length(breaks) == 0))
		stop(
"You must have at least two breaks.\n");

	if (is.null(propVals))
		stop(
"'propVals' cannot be NULL.\n");

	if (length(breaks) == 1) {
		numBreaks <- breaks;
	} else {
		numBreaks <- length(breaks) - 1;
	}

	propVals <- rep(propVals, length.out=numBreaks);
	propIdx <- cut(pdata$Z, breaks, include.lowest=TRUE);
	pdata[[propName]] <- propVals[unclass(propIdx)];

	if (!is.null(pdata) &&
			!is.PolyData(pdata, fullValidation=FALSE))
		class(pdata) <- c("PolyData", class(pdata));

	return(pdata);
}

##==============================================================================
plotLines <- function(polys, xlim=NULL, ylim=NULL, projection=FALSE,
											plt=c(0.11, 0.98, 0.12, 0.88), polyProps=NULL,
											lty=NULL, col=NULL, bg=0, axes=TRUE,
											tckLab=TRUE, tck=0.014, tckMinor=0.5 * tck, ...)
	## The only layout graphics parameter accepted by 'plotLines()' is 'plt'.
	## Since the function changes the plot region to satisfy the aspect ratio,
	## it is assigned a default value to improve consistency among plots.
{
	r <- .plotMaps(polys=polys, xlim=xlim, ylim=ylim, projection=projection,
								 plt=plt, polyProps=polyProps, border=NULL, lty=lty, col=col,
								 density=NULL, angle=NULL, bg=bg, axes=axes, tckLab=tckLab,
								 tck=tck, tckMinor=tckMinor, isType="lines", ...);

	if (!is.null(r) &&
			!is.PolyData(r, fullValidation=FALSE))
		class(r) <- c("PolyData", class(r));

	invisible(r);
}

##==============================================================================
plotMap <- function(polys, xlim=NULL, ylim=NULL, projection=TRUE,
										plt=c(0.11, 0.98, 0.12, 0.88), polyProps=NULL,
										border=NULL, lty=NULL, col=NULL, colHoles=NULL,
										density=NA, angle=NULL, bg=0, axes=TRUE,
										tckLab=TRUE, tck=0.014, tckMinor=0.5 * tck, ...)
	## The only layout graphics parameter accepted by 'plotMap()' is 'plt'.
	## Since the function changes the plot region to satisfy the aspect ratio,
	## it is assigned a default value to improve consistency among plots.
{
	r <- .plotMaps(polys=polys, xlim=xlim, ylim=ylim, projection=projection,
								 plt=plt, polyProps=polyProps, border=border, lty=lty, col=col,
								 colHoles=colHoles, density=density, angle=angle, bg=bg,
								 axes=axes, tckLab=tckLab, tck=tck, tckMinor=tckMinor,
								 isType="polygons", ...);

	if (!is.null(r) &&
			!is.PolyData(r, fullValidation=FALSE))
		class(r) <- c("PolyData", class(r));

	invisible(r);
}

##==============================================================================
plotPoints <- function(data, xlim=NULL, ylim=NULL, projection=FALSE,
											 plt=c(0.11, 0.98, 0.12, 0.88), polyProps=NULL,
											 cex=NULL, col=NULL, pch=NULL, axes=TRUE,
											 tckLab=TRUE, tck=0.014, tckMinor=0.5 * tck, ...)
	## The only layout graphics parameter accepted by 'plotPoints()' is 'plt'.
	## Since the function changes it to satisfy the aspect ratio, it's assigned
	## a default value.
{
	r <- .plotMaps(polys=data, xlim=xlim, ylim=ylim, projection=projection,
								 plt=plt, polyProps=polyProps, border=NULL, lty=NULL, col=col,
								 density=NULL, angle=NULL, bg=NULL,
								 cex=cex, pch=pch, # both become part of '...'
								 axes=axes, tckLab=tckLab, tck=tck, tckMinor=tckMinor,
								 isType="points", ...);

	if (!is.null(r) &&
			!is.PolyData(r, fullValidation=FALSE))
		class(r) <- c("PolyData", class(r));

	invisible(r)
}

##==============================================================================
plotPolys <- function(polys, xlim=NULL, ylim=NULL, projection=FALSE,
											plt=c(0.11, 0.98, 0.12, 0.88), polyProps=NULL,
											border=NULL, lty=NULL, col=NULL, colHoles=NULL,
											density=NA, angle=NULL, bg=0, axes=TRUE,
											tckLab=TRUE, tck=0.014, tckMinor=0.5 * tck, ...)
	## The only layout graphics parameter accepted by 'plotPolys()' is 'plt'.
	## Since the function changes it to satisfy the aspect ratio, it's assigned
	## a default value.
{
	r <- .plotMaps(polys=polys, xlim=xlim, ylim=ylim, projection=projection,
								 plt=plt, polyProps=polyProps, border=border, lty=lty, col=col,
								 colHoles=colHoles, density=density, angle=angle, bg=bg,
								 axes=axes, tckLab=tckLab, tck=tck, tckMinor=tckMinor,
								 isType="polygons", ...);

	if (!is.null(r) &&
			!is.PolyData(r, fullValidation=FALSE))
		class(r) <- c("PolyData", class(r));

	invisible(r)
}

##==============================================================================
print.EventData <- function(x, ...)
{
	if (exists("PBSprint") && (PBSprint || .PBSmapEnv$PBSprint)) {
		print(summary.EventData(x), ...);
	} else {
		print(data.frame(unclass(x)));
	}
	invisible(x);
}

##==============================================================================
print.LocationSet <- function(x, ...)
{
	if (exists("PBSprint") && (PBSprint || .PBSmapEnv$PBSprint)) {
		print(summary.LocationSet(x), ...);
	} else {
		print(data.frame(unclass(x)));
	}
	invisible(x);
}

##==============================================================================
print.PolyData <- function(x, ...)
{
	if (exists("PBSprint") && (PBSprint || .PBSmapEnv$PBSprint)) {
		print(summary.PolyData(x), ...);
	} else {
		print(data.frame(unclass(x)));
	}
	invisible(x);
}

##==============================================================================
print.PolySet <- function(x, ...)
{
	if (exists("PBSprint") && (PBSprint || .PBSmapEnv$PBSprint)) {
		print(summary.PolySet(x), ...);
	} else {
		print(data.frame(unclass(x)));
	}
	invisible(x);
}

##==============================================================================
print.summary.PBS <- function(x, ...)
{
	str <-
		c(x$type, "",
			"", "",
			"Records				 : ", ifelse(is.null(x$records),"NULL", x$records),
			"	Contours			: ", ifelse(x$type == "EventData",
																	 "NA",
																	 ifelse(is.null(x$contours), "NULL",
																					x$contours)),
			"		Holes			 : ", ifelse(x$type != "PolySet",
																	 "NA",
																	 ifelse(is.null(x$holes), "NULL", x$holes)),
			"	Events				: ", ifelse(x$type == "EventData",
																	 ifelse(is.null(x$events), "NULL", x$events),
																	 "NA"),
			"		On boundary : ", ifelse(x$type != "LocationSet",
																	 "NA",
																	 ifelse(is.null(x$boundary), "NULL",
																					x$boundary)),
			"", "",
			"Ranges", "",
			"	X						 : ", ifelse(is.null(x$range.x), "NULL",
																	 paste("[",
																				 paste(range(x$range.x),collapse=", "),
																				 "]", sep="")),
			"	Y						 : ", ifelse(is.null(x$range.y), "NULL",
																	 paste("[",
																				 paste(range(x$range.y),collapse=", "),
																				 "]", sep="")),
			"", "",
			"Attributes", "",
			"	Projection		: ", ifelse(is.null(x$attr.projection), "NULL",
																	 x$attr.projection),
			"	Zone					: ", ifelse(is.null(x$attr.zone), "NULL",
																	 x$attr.zone),
			"", "",
			"Extra columns	 : ", ifelse(is.null(x$columns), "NULL", x$columns));

	for (i in seq(1, by=2, length.out=(length(str)/2))) {
		cat (str[i], str[i+1], "\n", sep="");
	}
	invisible(x);
}


## refocusWorld-------------------------2019-03-14
##	Refocus the 'worldLL'/'worldLLhigh' data sets.
## ------------------------------------------NB|RH
refocusWorld <- function(polys, xlim=NULL, ylim=NULL, clip.AN=TRUE)
{
	## Define a local function to assist with shifting
	.shiftRegion <- function(polys, shift=0) {
		npolys			 <- polys;
		npolys$xmin	<- polys$xmin + (360 * shift);
		npolys$xmax	<- polys$xmax + (360 * shift);
		npolys$shift <- polys$shift + shift;
		return (npolys);
	}
	## Validate input
	polys <- .validatePolySet(polys);
	if (is.character(polys))
		stop(paste("Invalid PolySet 'polys'.\n", polys, sep=""));

	## Default values
	if (is.null(xlim))
		xlim <- range(polys$X)
	if (is.null(ylim))
		ylim <- range(polys$Y)

	## Validate/adjust xlim
	if (abs(diff(xlim)) > 360)
		stop("Invalid 'xlim' range: cannot exceed 360 degrees.\n");
	warn <- FALSE;
	while (min(xlim) < -360) {
		xlim <- xlim + 360;
		warn <- TRUE;
	}
	while (max(xlim) > 360) {
		xlim <- xlim - 360;
		warn <- TRUE;
	}
	if (warn)
		warning(paste("Shifted 'xlim' to ", paste(xlim, collapse=","), ": use new limits in future plots.\n", sep=""));

	## Save the attributes of the data frame
	attrNames	<- setdiff(names(attributes(polys)), c("names", "row.names", "class"));
	attrValues <- attributes(polys)[attrNames];
	
	## Build indices in case both PID *and* SID exist
	polys$IDs <- .createIDs(polys, cols=c("PID", "SID"));
	
	## Build a data frame with the ranges
	xrng <- lapply(split(polys[, c("X")], polys$IDs), range)
	yrng <- lapply(split(polys[, c("Y")], polys$IDs), range)
	ppid <- lapply(split(polys[, c("PID")], polys$IDs), unique)
	npol <- lapply(split(polys[, c("PID")], polys$IDs), length)
	rng	<- data.frame(IDs=names(xrng),
		PID =unlist(ppid),
		xmin=unlist(xrng)[c(TRUE,FALSE)],
		xmax=unlist(xrng)[c(FALSE,TRUE)],
		ymin=unlist(yrng)[c(TRUE,FALSE)],
		ymax=unlist(yrng)[c(FALSE,TRUE)], shift=0);

	## 'rng' currently contains one row per polygon; to wrap the region,
	## let's make the range [-360,360] contain *every* polygon that exists
	## within that range -- doing so will require replicating polygons;
	## for "shifted" polygons, we'll add a "shift" field to the data frame
	nrng <- rng;	# 'nrng' will be the new range data frame

	## Go right until every polygon's X is greater than the range
	srng <- .shiftRegion (rng, 1);
	while (any(srng$xmin < 360)) {
		# typically one iteration
		nrng <- rbind (nrng, srng);
		srng <- .shiftRegion (srng, 1);
	}
	## Go left until every polygon's X is less than the range
	srng <- .shiftRegion (rng, -1);
	while (any(srng$xmax > -360)) {
		# typically one iteration
		nrng <- rbind(nrng, srng);
		srng <- .shiftRegion (srng, -1);
	}

	## We've added polygons that are completely outside of the range; remove those polygons
	nrng <- nrng[nrng$xmax > -360 & nrng$xmin < 360, ]

	## Identify all of the polygons within the region of interest (x/ylim)
	nrng <- nrng[!(nrng$xmax < xlim[1] | nrng$xmin > xlim[2] |
		nrng$ymax < ylim[1] | nrng$ymin > ylim[2]), ];

	## Eliminate duplicated polys for now, i.e., Antarctica (dealt with below if necessary)
	dups <- duplicated(nrng$IDs);
	if (any(dups)) {
		warning(paste("Removed duplicates of following polygons (Antarctica?): ",
			paste(unique(nrng$PID[dups]), collapse=", "), "\n", sep=""));
		nrng <- nrng[!dups, ];
	}
	## Given relevant PIDs and their shift factors, let's create a new data
	## frame containing only those polygons, where we've shifted the X values
	idx <- split(1:nrow(polys), polys$IDs)		## to extract relevant polygons
	len <- lapply(idx, length)								## num vertices in each polygon

	## Indices for relevant polys
	idx <- unlist(idx[as.character(nrng$IDs)])
	## Number vertices in relevant polys
	len <- unlist(len[as.character(nrng$IDs)])
	
	## Build the return value
	npolys	 <- polys[idx, setdiff (names (polys), "IDs")];
	npolys$X <- npolys$X + rep(nrng$shift * 360, len);

	## Deal with Antarctica (RH 181026)
	is.Ant	=sapply(yrng,function(x){all(x< -60)}) & sapply(npol,function(x){x>100})
	if (any(is.Ant)) { ## borders don't usually extend into Antarctica (RH 190314)
		pid.Ant =ppid[[names(is.Ant)[is.Ant]]]
		if (pid.Ant %in% nrng$PID) {
			z.Ant	 =is.element(npolys$PID,pid.Ant) & npolys$X>=0 & npolys$X<=360 & !is.na(npolys$X) & npolys$Y> -90 & !is.na(npolys$Y)
			p0.Ant	=npolys[z.Ant,]
			## Reverse direction of polygon from west to east (if necessary)
			if (p0.Ant$X[1] > rev(p0.Ant$X)[1])
				p0.Ant=p0.Ant[nrow(p0.Ant):1,]
			## Get the Antarctic base polygon and add duplicates to West and East
			pW.Ant	=pE.Ant=p0.Ant
			pW.Ant$X=pW.Ant$X - 360
			pE.Ant$X=pE.Ant$X + 360
			vW.Ant	=pW.Ant[1,,drop=FALSE]
			vE.Ant	=pE.Ant[nrow(pE.Ant),,drop=FALSE]
			vW.Ant$Y=vE.Ant$Y=-90
			## Stitch together lower W vertex, three Antarcticas (West, Central, East), and lower E vertex
			big.Ant =rbind(vW.Ant, pW.Ant, p0.Ant, pE.Ant, vE.Ant)
			big.Ant$POS=1:nrow(big.Ant)
			## Replace Antarctica
			npolys	=npolys[!is.element(npolys$PID,pid.Ant),]
			if (clip.AN) {
				clip.Ant=clipPolys(big.Ant, xlim=range(npolys$X), ylim=ylim)	## Antarctica clipped
				npolys	=rbind(npolys, clip.Ant[,colnames(npolys)])
			} else {
				npolys	=rbind(npolys, big.Ant[,colnames(npolys)])
			}
		}
	}
	row.names(npolys) <- NULL;
	## Restore the attributes
	attributes(npolys) <- c(attributes(npolys), attrValues);
	attr(npolys,"rf.xlim")=xlim # RH added for use in importGSHHS
	invisible(npolys);
}
##~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~refocusWorld


##==============================================================================
summary.EventData <- function(object, ...)
{
	z <- list();
	z$type <- "EventData";
	z$records <- nrow(object);
	z$events <- length(unique(object$EID));
	z$columns <- paste(setdiff(names(object), c("EID", "X", "Y")), collapse=", ");
	z$range.x <- range(object$X);
	z$range.y <- range(object$Y);
	z$attr.projection <- attr(object, "projection");
	z$attr.zone <- attr(object, "zone");
	class(z) <- "summary.PBS";

	z;
}

##==============================================================================
summary.LocationSet <- function(object, ...)
{
	## create 'IDX' vector used later for 'unique'/'duplicated'
	IDX <- .createIDs(object, cols=c("PID", "SID"));

	z <- list();
	z$type <- "LocationSet";
	z$events <- length(unique(object$EID));
	z$contours <- length(unique(IDX));
	z$records <- nrow(object);
	z$boundary <- sum(object$Bdry);
	z$attr.projection <- attr(object, "projection");
	z$attr.zone <- attr(object, "zone");
	z$columns <- paste(setdiff(names(object), c("EID", "PID", "SID", "Bdry")),
										 collapse=", ");
	class(z) <- "summary.PBS";

	z;
}

##==============================================================================
summary.PolyData <- function(object, ...)
{
	## create 'IDX' vector used later for 'unique'/'duplicated'
	IDX <- .createIDs(object, cols=c("PID", "SID"));

	z <- list();
	z$type <- "PolyData";
	z$records <- nrow(object);
	z$contours <- length(unique(IDX));
	z$columns <- paste(setdiff(names(object), c("PID", "SID")), collapse=", ");
	z$attr.projection <- attr(object, "projection");
	z$attr.zone <- attr(object, "zone");
	class(z) <- "summary.PBS";

	z;
}

##==============================================================================
summary.PolySet <- function(object, ...)
{
	## create 'IDX' vector used later for 'unique'/'duplicated'
	## create 'holes' (Boolean) vector to count how many holes exist
	IDX <- .createIDs(object, cols=c("PID", "SID"));
	if (is.element("SID", names(object))) {
		idxFirst <- which(!duplicated(IDX));
		idxLast <- c((idxFirst-1)[-1], length(IDX));
		holes <- object$POS[idxFirst] > object$POS[idxLast];
	} else {
		holes <- NULL;
	}

	z <- list();
	z$type <- "PolySet";
	z$records <- nrow(object);
	z$contours <- length(unique(IDX));
	z$holes <- sum(holes);
	z$range.x <- range(object$X);
	z$range.y <- range(object$Y);
	z$attr.projection <- attr(object, "projection");
	z$attr.zone <- attr(object, "zone");
	z$columns <- paste(setdiff(names(object), c("PID", "SID", "POS", "X", "Y")),
										 collapse=", ");
	class(z) <- "summary.PBS";

	z;
}

##==============================================================================
thickenPolys <- function(polys, tol=1, filter=3, keepOrig=TRUE,
												 close=TRUE)
{
	polys <- .validatePolySet(polys);
	if (is.character(polys))
		stop(paste("Invalid PolySet 'polys'.\n", polys, sep=""));

	## PART 1 --------------------------------------------------------------------
	## - estimate 'outCapacity'
	## - this should use the same distance formula as the underlying C code;
	##	 otherwise, the estimation may be wrong

	## create an 'idx' vector for working with 'polys'
	idx <- .createIDs(polys, cols=c("PID", "SID"));

	## if closing polygons, ensure they all close before we calculate
	## 'outCapacity'; close them all, even if already closed, as a
	## precaution
	if (close) {
		## find start and end of each polygon
		starts <- which(!duplicated(idx));
		## special case of a single polygon
		if (length(starts) == 1) {
			ends <- length(idx);
		} else {
			ends <- c((starts[-1]) - 1, length(idx));
		}

		## create a vector of [start #1, end #1, start #2, end #2, ...]
		toDupe <- vector();
		toDupe[seq(1, by=2, length=length(starts))] <- starts;
		toDupe[seq(2, by=2, length=length(ends))] <- ends;

		## determine next avail. PID
		newPID <- max(polys$PID) + 1;

		## create a new PolySet to use in these calculations
		newRows <- polys[toDupe, ];
		newRows$PID <- rep((newPID:(newPID + length(starts) - 1)), each=2)
		tempPolys <- rbind(polys, newRows);

		## update 'idx'; it won't matter that the new 'idx' values
		## are a different format since they're using new PIDs (i.e., they
		## don't need to match anything); they just need to be relative to
		## earlier 'idx' value, and they will be since we're using a 'paste'
		if (is.element("SID", names(newRows))) {
			idx <- c(idx, paste(newRows$PID, newRows$SID));
		} else {
			idx <- c(idx, newRows$PID);
		}
	} else {
		tempPolys <- polys;
	}

	## compute a distance vector
	projection <- attr(polys, "projection");
	if (!is.null(projection) && !is.na(projection)
			&& ((projection == "UTM") || (projection == "LL")
					|| (projection == 1))) {
		dist <- .calcDist(tempPolys);
	}
	else {
		stop(paste(
"Invalid projection attribute.	Supported projections include \"LL\",",
"\"UTM\", and 1.\n"));
	}
	## filter the distances between polygons (make them not count...)
	dist[((which(!duplicated(idx)) - 1)[-1])] <- 0;

	if (keepOrig == TRUE) {
		## since original points are kept, calculate the capacity based on how
		## many points will be added between the original points
		outCapacity <- sum(ceiling(dist / tol));
	} else {
		## since original points are not kept, calculate the length of each
		## polygon, and then determine how many points will be required to
		## represent that polygon
		dist <- unlist(lapply(split(dist, idx), sum), use.names=FALSE);
		outCapacity <- sum(ceiling(dist / tol));
	}

	## add a fudge factor
	outCapacity <- outCapacity + 1000;

	## PART 2 --------------------------------------------------------------------
	## - thicken the PolySet

	pLL <- 0;
	pUTM <- 1;

	projection <- attributes(polys)$projection;
	if (!is.null(projection) && !is.na(projection) && projection == "LL") {
		proj <- pLL;
	} else {
		proj <- pUTM;
	}

	## save the attributes of the data frame (.validatePolySet returns a data
	## frame); in this case, SAVE ALL ATTRIBUTES since we are "modifying" an
	## existing data set
	attrNames <- setdiff(names(attributes(polys)),
											 c("names", "row.names", "class"));
	attrValues <- attributes(polys)[attrNames];

	inRows <- nrow(polys);

	## create the data structures that the C function expects
	if (!is.element("SID", names(polys))) {
		inID <- c(polys$PID, integer(length=inRows), polys$POS);
	} else {
		inID <- c(polys$PID, polys$SID, polys$POS);
	}
	inXY <- c(polys$X, polys$Y);

	## call the C function
	results <- .C("thickenPolys",
								inID=as.integer(inID),
								inXY=as.double(inXY),
								inVerts=as.integer(inRows),
								tolerance=as.double(tol),
								filter=as.integer(filter),
								units=as.integer(proj),
								keepOrig=as.integer(keepOrig),
								close=as.integer(close),
								outID=integer(3 * outCapacity),
								outXY=double(2 * outCapacity),
								outRows=as.integer(outCapacity),
								outStatus=integer(1),
								PACKAGE="PBSmapping");
	## note: outRows is set to how much space is allocated -- the C function
	##			 should take this into consideration

	if (results$outStatus == 1) {
		stop(
"Insufficient physical memory for processing.\n");
	}
	if (results$outStatus == 2) {
		stop(paste(
"Insufficient memory allocated for output.	Please upgrade to the latest",
"version of the software, and if that does not fix this problem, please",
"file a bug report.\n",
							 sep="\n"));
	}

	## determine the number of rows in the result
	outRows <- as.vector(results$outRows);

	## extract the data from the C function results
	if (outRows > 0) {
		d <- data.frame(PID=results$outID[1:outRows],
										SID=results$outID[(outCapacity+1):(outCapacity+outRows)],
										POS=results$outID[(2*outCapacity+1):(2*outCapacity+outRows)],
										X=results$outXY[1:outRows],
										Y=results$outXY[(outCapacity+1):(outCapacity+outRows)]);

		if (!is.element("SID", names(polys)))
			d$SID <- NULL;

		## restore the attributes (including projection and zone)
		attributes(d) <- c(attributes(d), attrValues);

		if (!is.PolySet(d, fullValidation=FALSE))
			class(d) <- c("PolySet", class(d));

		return(d);
	} else {
		return(NULL);
	}
}

##==============================================================================
thinPolys <- function(polys, tol=1, filter=3)
{
	polys <- .validatePolySet(polys);
	if (is.character(polys))
		stop(paste("Invalid PolySet 'polys'.\n", polys, sep=""));

	pLL <- 0;
	pUTM <- 1;

	projection <- attributes(polys)$projection;
	if (!is.null(projection) && !is.na(projection) && projection == "LL") {
		proj <- pLL;
	} else {
		proj <- pUTM;
	}

	## save the attributes of the data frame (.validatePolySet returns a data
	## frame); in this case, SAVE ALL ATTRIBUTES since we are "modifying" an
	## existing data set
	attrNames <- setdiff(names(attributes(polys)),
											 c("names", "row.names", "class"));
	attrValues <- attributes(polys)[attrNames];

	inRows <- nrow(polys);
	outCapacity <- as.integer(inRows);

	## create the data structures that the C function expects
	if (!is.element("SID", names(polys))) {
		inID <- c(polys$PID, integer(length=inRows), polys$POS);
	} else {
		inID <- c(polys$PID, polys$SID, polys$POS);
	}

	## adjust units to what the C function expects
	if (proj == pLL) {
		## degrees >> micro-degrees
		inXY <- c(polys$X, polys$Y) * 10^6;
	} else {
		## kilometers >> meters
		inXY <- c(polys$X, polys$Y) * 10^3;
	}

	## call the C function
	results <- .C("thinPolys",
								inID=as.integer(inID),
								inXY=as.integer(inXY),
								inVerts=as.integer(inRows),
								tolerance=as.double(tol),
								filter=as.integer(filter),
								units=as.integer(proj),
								outID=integer(3 * outCapacity),
								outXY=integer(2 * outCapacity),
								outRows=as.integer(outCapacity),
								outStatus=integer(1),
								PACKAGE="PBSmapping");
	## note: outRows is set to how much space is allocated -- the C function
	##			 should take this into consideration

	if (results$outStatus == 1) {
		stop(
"Insufficient physical memory for processing.\n");
	}
	if (results$outStatus == 2) {
		stop(paste(
"Insufficient memory allocated for output.	Please upgrade to the latest",
"version of the software, and if that does not fix this problem, please",
"file a bug report.\n",
							 sep="\n"));
	}

	## determine the number of rows in the result
	outRows <- as.vector(results$outRows);

	## extract the data from the C function results
	if (outRows > 0) {
		if (proj == pLL) {
			results$outXY <- results$outXY / 10^6;
		} else {
			results$outXY <- results$outXY / 10^3;
		}

		d <- data.frame(PID=results$outID[1:outRows],
										SID=results$outID[(outCapacity+1):(outCapacity+outRows)],
										POS=results$outID[(2*outCapacity+1):(2*outCapacity+outRows)],
										X=results$outXY[1:outRows],
										Y=results$outXY[(outCapacity+1):(outCapacity+outRows)]);

		if (!is.element("SID", names(polys)))
			d$SID <- NULL;

		## restore the attributes (including projection and zone)
		attributes(d) <- c(attributes(d), attrValues);

		if (!is.PolySet(d, fullValidation=FALSE))
			class(d) <- c("PolySet", class(d));

		return(d);
	} else {
		return(NULL);
	}
}