| Title: | Mapping Fisheries Data and Spatial Analysis Tools |
|---|---|
| Description: | This software has evolved from fisheries research conducted at the Pacific Biological Station (PBS) in 'Nanaimo', British Columbia, Canada. It extends the R language to include two-dimensional plotting features similar to those commonly available in a Geographic Information System (GIS). Embedded C code speeds algorithms from computational geometry, such as finding polygons that contain specified point events or converting between longitude-latitude and Universal Transverse Mercator (UTM) coordinates. Additionally, we include 'C++' code developed by Angus Johnson for the 'Clipper' library, data for a global shoreline, and other data sets in the public domain. Under the user's R library directory '.libPaths()', specifically in './PBSmapping/doc', a complete user's guide is offered and should be consulted to use package functions effectively. |
| Authors: | Jon T. Schnute [aut], Nicholas Boers [aut], Rowan Haigh [aut, cre], Alex Couture-Beil [ctb], Denis Chabot [ctb], Chris Grandin [ctb], Alan Murta [ctb], Angus Johnson [ctb], Paul Wessel [ctb], Franklin Antonio [ctb], Nicholas J. Lewin-Koh [ctb], Roger Bivand [ctb], Sean Anderson [ctb] |
| Maintainer: | Rowan Haigh <[email protected]> |
| License: | GPL (>= 2) |
| Version: | 2.74.1 |
| Built: | 2024-10-25 06:18:06 UTC |
| Source: | https://github.com/pbs-software/pbs-mapping |
Add bubbles proportional to some EventData's Z column
(e.g., catch or effort) to an existing plot, where each unique
EID describes a bubble.
addBubbles(events, type=c("perceptual","surface","volume"), z.max=NULL, min.size=0, max.size=0.8, symbol.zero="+", symbol.fg=rgb(0,0,0,0.60), symbol.bg=rgb(0,0,0,0.30), legend.pos="bottomleft", legend.breaks=NULL, show.actual=FALSE, legend.type=c("nested","horiz","vert"), legend.title="Abundance", legend.cex=.8, neg.col="RYB", ...)addBubbles(events, type=c("perceptual","surface","volume"), z.max=NULL, min.size=0, max.size=0.8, symbol.zero="+", symbol.fg=rgb(0,0,0,0.60), symbol.bg=rgb(0,0,0,0.30), legend.pos="bottomleft", legend.breaks=NULL, show.actual=FALSE, legend.type=c("nested","horiz","vert"), legend.title="Abundance", legend.cex=.8, neg.col="RYB", ...)
events |
EventData to use (required). |
type |
|
z.max |
|
min.size |
|
max.size |
|
symbol.zero |
|
symbol.fg |
|
symbol.bg |
|
legend.pos |
|
legend.breaks |
|
show.actual |
|
legend.type |
|
legend.title |
|
legend.cex |
|
neg.col |
|
... |
|
Modified from (and for the legend, strongly inspired by) Tanimura et al. (2006) by Denis Chabot to work with PBSmapping.
Furthermore, Chabot's modifications make it possible to draw
several maps with bubbles that all have the same scale
(instead of each bubble plot having a scale that depends on
the maximum z-value for that plot).
This is done by making z.max equal to the largest z-value
from all maps that will be plotted.
The user can also add a legend in one of four corners
(see legend) or at a specific c(X,Y) position.
If legend.pos is NULL, no legend is drawn.
Denis Chabot, Research Scientist
Maurice-Lamontagne Institute, Fisheries & Oceans Canada (DFO), Mont-Joli QC
Maintainer: Rowan Haigh, Program Head – Offshore Rockfish
Pacific Biological Station (PBS), Fisheries & Oceans Canada (DFO), Nanaimo BC
locus opus: Offsite, Vancouver BC
Last modified Rd: 2024-09-25
Tanimura, S., Kuroiwa, C., and Mizota, T. (2006) Proportional symbol mapping in R. Journal of Statistical Software 15(5).
In package PBSmapping:addPolys,
EventData,
RGB2RYB,
surveyData
local(envir=.PBSmapEnv,expr={ oldpar = par(no.readonly=TRUE) #--- common code for both examples below data(nepacLL,surveyData,envir=.PBSmapEnv) surveyData$Z <- surveyData$catch #--- plot a version that only varies the size plotMap(nepacLL, xlim=c(-131.8,-127.2), ylim=c(50.5,52.7), col="gainsboro",plt=c(.08,.99,.08,.99), cex.axis=1.2, cex.lab=1.5) addBubbles(surveyData, symbol.bg=rgb(.9,.5,0,.6), legend.type="nested", symbol.zero="+", col="grey") #--- plot a version that uses different symbol colours plotMap(nepacLL, xlim=c(-131.8,-127.2), ylim=c(50.5,52.7), col="gainsboro",plt=c(.08,.99,.08,.99), cex.axis=1.2, cex.lab=1.5) subset <- surveyData[surveyData$Z <= 1000, ] addBubbles(subset, symbol.bg=c("red", "yellow", "green"), legend.type="horiz", legend.breaks=pretty(range(subset$Z), n=11), symbol.zero=FALSE, col="grey", min.size=0.1, max.size=0.4) par(oldpar) })local(envir=.PBSmapEnv,expr={ oldpar = par(no.readonly=TRUE) #--- common code for both examples below data(nepacLL,surveyData,envir=.PBSmapEnv) surveyData$Z <- surveyData$catch #--- plot a version that only varies the size plotMap(nepacLL, xlim=c(-131.8,-127.2), ylim=c(50.5,52.7), col="gainsboro",plt=c(.08,.99,.08,.99), cex.axis=1.2, cex.lab=1.5) addBubbles(surveyData, symbol.bg=rgb(.9,.5,0,.6), legend.type="nested", symbol.zero="+", col="grey") #--- plot a version that uses different symbol colours plotMap(nepacLL, xlim=c(-131.8,-127.2), ylim=c(50.5,52.7), col="gainsboro",plt=c(.08,.99,.08,.99), cex.axis=1.2, cex.lab=1.5) subset <- surveyData[surveyData$Z <= 1000, ] addBubbles(subset, symbol.bg=c("red", "yellow", "green"), legend.type="horiz", legend.breaks=pretty(range(subset$Z), n=11), symbol.zero=FALSE, col="grey", min.size=0.1, max.size=0.4) par(oldpar) })
Add a compass rose to an existing map, similar to those found on nautical charts showing both true north and magnetic north.
addCompass(X, Y, rot="magN", useWest=TRUE, year, cex=1, col.compass=c("gainsboro","blue","yellow","black"), ...)addCompass(X, Y, rot="magN", useWest=TRUE, year, cex=1, col.compass=c("gainsboro","blue","yellow","black"), ...)
X |
|
Y |
|
rot |
|
useWest |
|
year |
|
cex |
|
col.compass |
|
... |
|
The basic idea comes from Jim Lemon (see References), but is modified here to reflect a compass rose used on BC nautical charts.
The default rotation ("magN") is a calculation of the initial
bearing of a great-circle arc from the compass position to the north
geomagnetic rot using the function calcGCdist.
The default year is the current year, but the user can choose years from
1900 to 2025 for approximate rot locations using model output from NOAA's IGRF-13
(International Geomagnetic Reference Field).
The user can also specify a fixed rotation (e.g. rot=-30) or no
rotation (either rot=0 or rot="trueN").
No value returned.
Rowan Haigh, Program Head – Offshore Rockfish
Pacific Biological Station (PBS), Fisheries & Oceans Canada (DFO), Nanaimo BC
locus opus: Regional Headquarters, Vancouver BC
Last modified Rd: 2022-07-05
[R-sig-Geo] How to display a compass rose on a map
Magnetic North, Geomagnetic and Magnetic Poles
addBubbles,
addLabels,
addPoints,
addStipples,
calcGCdist
local(envir=.PBSmapEnv,expr={ data(nepacLL,envir=.PBSmapEnv) par(mfrow=c(1,1),mar=c(3,4,0.5,0.5)) plotMap(nepacLL, xlim=c(-134.5,-124.5), ylim=c(48,55), plt=NULL, col="lightyellow", cex.axis=1.2, cex.lab=1.5) addCompass(-132, 49.5, rot=-12, cex=1.5) })local(envir=.PBSmapEnv,expr={ data(nepacLL,envir=.PBSmapEnv) par(mfrow=c(1,1),mar=c(3,4,0.5,0.5)) plotMap(nepacLL, xlim=c(-134.5,-124.5), ylim=c(48,55), plt=NULL, col="lightyellow", cex.axis=1.2, cex.lab=1.5) addCompass(-132, 49.5, rot=-12, cex=1.5) })
Add the label column of data to the existing plot.
addLabels (data, xlim = NULL, ylim = NULL, polyProps = NULL, placement = "DATA", polys = NULL, rollup = 3, cex = NULL, col = NULL, font = NULL, ...)addLabels (data, xlim = NULL, ylim = NULL, polyProps = NULL, placement = "DATA", polys = NULL, rollup = 3, cex = NULL, col = NULL, font = NULL, ...)
data |
|
xlim |
range of X-coordinates. |
ylim |
range of Y-coordinates. |
polyProps |
PolyData specifying which labels to plot and their
properties. |
placement |
one of |
polys |
PolySet to use for calculating label placement. |
rollup |
level of detail at which to process |
cex |
vector describing character expansion factors (cycled by
|
col |
vector describing colours (cycled by |
font |
vector describing fonts (cycled by |
... |
If data is EventData, it must minimally contain the columns
EID, X, Y, and label. Since the
EID column does not match a column in polys, set
placement = "DATA". The function plots each label at
its corresponding X/Y coordinate.
If data is PolyData, it must minimally contain the columns
PID and label. If it also contains X and
Y columns, set placement = "DATA" to plot labels at
those coordinates. Otherwise, set placement to one of
"CENTROID", "MEAN_RANGE", or "MEAN_XY". When
placement != "DATA", supply a PolySet polys. Using this
PolySet, the function calculates a centroid, mean range, or mean X/Y
coordinate for each polygon, and then links those PolyData with
data by PID/SID to determine label
coordinates.
If data contains both PID and EID columns, the
function assumes it is PolyData and ignores the EID
column.
For additional help on the arguments cex, col, and
font, please see par.
EventData or PolyData with X and Y columns
that can subsequently reproduce the labels on the plot. Modify this
data frame to tweak label positions.
Nicholas M. Boers, Staff Software Engineer
Jobber, Edmonton AB
Last modified Rd: 2013-04-10
addPoints,
calcCentroid,
calcMidRange,
calcSummary,
EventData,
plotPoints,
PolyData.
local(envir=.PBSmapEnv,expr={ oldpar = par(no.readonly=TRUE) #--- create sample PolyData to label Vancouver Island labelData <- data.frame(PID=33, label="Vancouver Island"); #--- load data if (!is.null(version$language) && (version$language == "R")) data(nepacLL,envir=.PBSmapEnv) #--- plot the map plotMap(nepacLL,xlim=c(-129,-122.6),ylim=c(48,51.1),col="lemonchiffon") #--- add the labels addLabels(labelData,placement="CENTROID",polys=nepacLL,cex=1.2,col=2,font=2) par(oldpar) })local(envir=.PBSmapEnv,expr={ oldpar = par(no.readonly=TRUE) #--- create sample PolyData to label Vancouver Island labelData <- data.frame(PID=33, label="Vancouver Island"); #--- load data if (!is.null(version$language) && (version$language == "R")) data(nepacLL,envir=.PBSmapEnv) #--- plot the map plotMap(nepacLL,xlim=c(-129,-122.6),ylim=c(48,51.1),col="lemonchiffon") #--- add the labels addLabels(labelData,placement="CENTROID",polys=nepacLL,cex=1.2,col=2,font=2) par(oldpar) })
Add a PolySet to an existing plot, where each unique (PID,
SID) describes a polyline.
addLines (polys, xlim = NULL, ylim = NULL, polyProps = NULL, lty = NULL, col = NULL, arrows = FALSE, ...)addLines (polys, xlim = NULL, ylim = NULL, polyProps = NULL, lty = NULL, col = NULL, arrows = FALSE, ...)
polys |
PolySet to add (required). |
xlim |
range of X-coordinates. |
ylim |
range of Y-coordinates. |
polyProps |
PolyData specifying which polylines to plot and their
properties. |
lty |
vector of line types (cycled by |
col |
vector of colours (cycled by |
arrows |
Boolean value; if |
... |
The plotting routine does not connect the last vertex of each discrete
polyline to the first vertex of that polyline. It clips polys
to xlim and ylim before plotting.
For additional help on the arguments lty and col, please
see par.
PolyData consisting of the PolyProps used to create the plot.
Nicholas M. Boers, Staff Software Engineer
Jobber, Edmonton AB
Last modified Rd: 2013-04-10
calcLength,
clipLines,
closePolys,
convLP,
fixBound,
fixPOS,
locatePolys,
plotLines,
thinPolys,
thickenPolys.
local(envir=.PBSmapEnv,expr={ oldpar = par(no.readonly=TRUE) #--- create a PolySet to plot polys <- data.frame(PID=rep(1,4),POS=1:4,X=c(0,1,1,0),Y=c(0,0,1,1)) polys <- as.PolySet(polys, projection=1) #--- plot the PolySet plotLines(polys, xlim=c(-.5,1.5), ylim=c(-.5,1.5), projection=1) #--- add the PolySet to the plot (in a different style) addLines(polys, lwd=5, col=3) par(oldpar) })local(envir=.PBSmapEnv,expr={ oldpar = par(no.readonly=TRUE) #--- create a PolySet to plot polys <- data.frame(PID=rep(1,4),POS=1:4,X=c(0,1,1,0),Y=c(0,0,1,1)) polys <- as.PolySet(polys, projection=1) #--- plot the PolySet plotLines(polys, xlim=c(-.5,1.5), ylim=c(-.5,1.5), projection=1) #--- add the PolySet to the plot (in a different style) addLines(polys, lwd=5, col=3) par(oldpar) })
Add EventData/PolyData to an existing plot, where each
unique EID describes a point.
addPoints (data, xlim = NULL, ylim = NULL, polyProps = NULL, cex = NULL, col = NULL, pch = NULL, ...)addPoints (data, xlim = NULL, ylim = NULL, polyProps = NULL, cex = NULL, col = NULL, pch = NULL, ...)
data |
|
xlim |
range of X-coordinates. |
ylim |
range of Y-coordinates. |
polyProps |
PolyData specifying which points to plot and their
properties. |
cex |
vector describing character expansion factors (cycled by
|
col |
vector describing colours (cycled by |
pch |
vector describing plotting characters (cycled by |
... |
This function clips data to xlim and ylim before
plotting. It only adds PolyData containing X and
Y columns.
For additional help on the arguments cex, col, and
pch, please see par.
PolyData consisting of the PolyProps used to create the plot.
Nicholas M. Boers, Staff Software Engineer
Jobber, Edmonton AB
Last modified Rd: 2013-04-10
combineEvents,
convDP,
findPolys,
locateEvents,
plotPoints.
local(envir=.PBSmapEnv,expr={ oldpar = par(no.readonly=TRUE) #--- load the data (if using R) if (!is.null(version$language) && (version$language=="R")) data(nepacLL,surveyData,envir=.PBSmapEnv) #--- plot a map plotMap(nepacLL, xlim=c(-136, -125), ylim=c(48, 57)) #--- add events addPoints(surveyData, col=1:7) par(oldpar) })local(envir=.PBSmapEnv,expr={ oldpar = par(no.readonly=TRUE) #--- load the data (if using R) if (!is.null(version$language) && (version$language=="R")) data(nepacLL,surveyData,envir=.PBSmapEnv) #--- plot a map plotMap(nepacLL, xlim=c(-136, -125), ylim=c(48, 57)) #--- add events addPoints(surveyData, col=1:7) par(oldpar) })
Add a PolySet to an existing plot, where each unique (PID,
SID) describes a polygon.
addPolys(polys, xlim = NULL, ylim = NULL, polyProps = NULL, border = NULL, lty = NULL, col = NULL, colHoles = NULL, density = NA, angle = NULL, ...) .insertNAs(polys, idx) .preparePolyProps(polysPID, polysSID, polyProps) .rollupPolys(polys, rollupMode, exteriorCCW, closedPolys, addRetrace)addPolys(polys, xlim = NULL, ylim = NULL, polyProps = NULL, border = NULL, lty = NULL, col = NULL, colHoles = NULL, density = NA, angle = NULL, ...) .insertNAs(polys, idx) .preparePolyProps(polysPID, polysSID, polyProps) .rollupPolys(polys, rollupMode, exteriorCCW, closedPolys, addRetrace)
polys |
|
xlim |
|
ylim |
|
polyProps |
|
border |
|
lty |
|
col |
|
colHoles |
|
density |
|
angle |
|
... |
|
idx |
|
polysPID |
|
polysSID |
|
rollupMode |
|
exteriorCCW |
|
closedPolys |
|
addRetrace |
|
The plotting routine connects the last vertex of each discrete polygon
to the first vertex of that polygon. It supports both
borders (border, lty) and fills (col,
density, angle). It clips polys to xlim
and ylim before plotting.
For additional help on the arguments border, lty,
col, density, and angle, please see
polygon and par.
PolyData consisting of the PolyProps used to create the plot.
Auxiliary dot function '.insertNAs' facilitates (hastens)
the plotting of polygons and polylines. It also reduces the incidence
of retrace lines.
Auxiliary dot function '.preparePolyProps' performs
at least one of the following tasks:
1) creates 'polyProps' if it equals NULL;
2) adds 'SID' column to 'polyProps' if one exists in 'polys';
3) removes from 'polyProps' any PIDs that do not exist in 'polys'.
Returns a polyProps object.
Auxiliary dot function '.rollupPolys' does not validate a PolySet;
returns a rolled-up PolySet or NULL.
Nicholas M. Boers, Software Engineer, Jobber, Edmonton AB
Maintainer: Rowan Haigh, Program Head – Offshore Rockfish
Pacific Biological Station (PBS), Fisheries & Oceans Canada (DFO), Nanaimo BC
locus opus: Remote office, Vancouver BC
Last modified Rd: 2024-09-03
In package PBSmapping:addLabels,
addLines (also uses '.preparePolyProps' and '.insertNAs'),
addPoints,
addStipples'.rollupPolys' is also called by:calcArea,
calcCentroid,
calcLength,
calcSummary,
fixPOS
local(envir=.PBSmapEnv,expr={ oldpar = par(no.readonly=TRUE) #--- create a PolySet to plot polys <- data.frame(PID=rep(1,4),POS=1:4,X=c(0,1,1,0),Y=c(0,0,1,1)) polys <- as.PolySet(polys, projection=1) #--- plot the PolySet plotPolys(polys,xlim=c(-.5,1.5),ylim=c(-.5,1.5),density=0,projection=1) #--- add the PolySet to the plot (in a different style) addPolys(polys,col="green",border="blue",lwd=3) par(oldpar) })local(envir=.PBSmapEnv,expr={ oldpar = par(no.readonly=TRUE) #--- create a PolySet to plot polys <- data.frame(PID=rep(1,4),POS=1:4,X=c(0,1,1,0),Y=c(0,0,1,1)) polys <- as.PolySet(polys, projection=1) #--- plot the PolySet plotPolys(polys,xlim=c(-.5,1.5),ylim=c(-.5,1.5),density=0,projection=1) #--- add the PolySet to the plot (in a different style) addPolys(polys,col="green",border="blue",lwd=3) par(oldpar) })
Add stipples to an existing plot.
addStipples (polys, xlim=NULL, ylim=NULL, polyProps=NULL, side=1, density=1, distance=4, ...)addStipples (polys, xlim=NULL, ylim=NULL, polyProps=NULL, side=1, density=1, distance=4, ...)
polys |
PolySet that provides the stipple boundaries (required). |
xlim |
range of X-coordinates. |
ylim |
range of Y-coordinates. |
polyProps |
PolyData specifying which polygons to stipple
and their properties. |
side |
one of |
density |
density of points, relative to the default. |
distance |
distance to offset points, measured as a percentage of
the absolute difference in |
... |
This function locates stipples based on the PolySet
polys and does not stipple degenerate lines.
PolyData consisting of the PolyProps used to create the plot.
Nicholas M. Boers, Staff Software Engineer
Jobber, Edmonton AB
Last modified Rd: 2013-04-10
addPoints,
addPolys,
plotMap,
plotPoints,
plotPolys,
points,
PolySet.
local(envir=.PBSmapEnv,expr={ oldpar = par(no.readonly=TRUE) #--- load the data (if using R) if (!is.null(version$language) && (version$language=="R")) data(nepacLL,envir=.PBSmapEnv) #--- plot a map plotMap(nepacLL,xlim=c(-128.66,-122.83),ylim=c(48.00,51.16)) #--- add stippling addStipples(nepacLL,col="purple",pch=20,cex=0.25,distance=2) par(oldpar) })local(envir=.PBSmapEnv,expr={ oldpar = par(no.readonly=TRUE) #--- load the data (if using R) if (!is.null(version$language) && (version$language=="R")) data(nepacLL,envir=.PBSmapEnv) #--- plot a map plotMap(nepacLL,xlim=c(-128.66,-122.83),ylim=c(48.00,51.16)) #--- add stippling addStipples(nepacLL,col="purple",pch=20,cex=0.25,distance=2) par(oldpar) })
Append a two-column matrix to a PolySet, assigning PID and
possibly SID values automatically or as specified in its
arguments.
appendPolys (polys, mat, PID = NULL, SID = NULL, isHole = FALSE)appendPolys (polys, mat, PID = NULL, SID = NULL, isHole = FALSE)
polys |
existing PolySet; if |
mat |
two-column matrix to append (required). |
PID |
new polygon's |
SID |
new polygon's |
isHole |
Boolean value; if |
If the PID argument is NULL, the appended polygon's
PID will be one greater than the maximum within polys
(if defined); otherwise, it will be 1.
If polys contains an SID column and the SID
argument equals NULL, this function uses the next available
SID for the corresponding PID.
If polys does not contain an SID column and the
caller passes an SID argument, all existing polygons will
receive an SID of 1. The new polygon's SID will
match the SID argument.
If isHole = TRUE, the polygon's POS values will
appropriately represent a hole (reverse order of POS).
If (PID, SID) already exists in the PolySet, the
function will issue a warning and duplicate those identifiers.
PolySet containing mat appended to polys. The
function retains attributes from polys.
Nicholas M. Boers, Staff Software Engineer
Jobber, Edmonton AB
Last modified Rd: 2013-04-10
addPolys,
clipPolys,
closePolys,
convLP,
fixBound,
fixPOS,
joinPolys,
plotMap,
plotPolys.
local(envir=.PBSmapEnv,expr={ #--- create two simple matrices a <- matrix(data=c(0,0,1,0,1,1,0,1),ncol=2,byrow=TRUE); b <- matrix(data=c(2,2,3,2,3,3,2,3), ncol=2,byrow=TRUE); #--- build a PolySet from them polys <- appendPolys(NULL, a); polys <- appendPolys(polys, b); #--- print the result print (polys); })local(envir=.PBSmapEnv,expr={ #--- create two simple matrices a <- matrix(data=c(0,0,1,0,1,1,0,1),ncol=2,byrow=TRUE); b <- matrix(data=c(2,2,3,2,3,3,2,3), ncol=2,byrow=TRUE); #--- build a PolySet from them polys <- appendPolys(NULL, a); polys <- appendPolys(polys, b); #--- print the result print (polys); })
Bathymetry data spanning British Columbia's coast.
data(bcBathymetry)data(bcBathymetry)
Three-element list: x = vector of horizontal grid line
locations, y = vector of vertical grid line locations, z
= (x by y) matrix containing water depths measured in
meters. Positive values indicate distance below sea level and
negative values above it.
The functions 'graphics::contour' and 'grDevices::contourLines' expect data in this format.
Function convCP converts the output from 'grDevices::contourLines' into a PolySet.
In R, the data must be loaded using the 'utils::data' function.
Nicholas M. Boers, Staff Software Engineer
Jobber, Edmonton AB
Maintainer: Rowan Haigh, Program Head – Offshore Rockfish
Pacific Biological Station (PBS), Fisheries & Oceans Canada (DFO), Nanaimo BC
locus opus: Regional Headquarters, Vancouver BC
Last modified Rd: 2022-07-07
Bathymetry data acquired from the Scripps Institution of Oceanography at the University of San Diego.
Using their online form, we requested bathymetry data for the complete
nepacLL region.
At forty megabytes, the data were not suitable for distribution in our mapping package.
Therefore, we reduced the data to the range
and .
Smith, W.H.F. and Sandwell, D.T. (1997)
Global seafloor topography from satellite altimetry and ship depth soundings.
Science 277, 1957-1962.
Website: https://topex.ucsd.edu/WWW_html/mar_topo.html
In package graphics:contour
In package grDevices:contourLines
In package PBSmapping:convCP,
nepacLL,
nepacLLhigh
Calculate the areas of polygons found in a PolySet.
calcArea (polys, rollup = 3)calcArea (polys, rollup = 3)
polys |
PolySet to use. |
rollup |
level of detail in the results; |
If rollup equals 1, the results contain an area for each
unique PID only. When it equals 2, they contain entries
for outer contours only. Finally, setting it to 3 prevents
roll-up, and they contain areas for each unique (PID,
SID).
Outer polygons have positive areas and inner polygons negative areas. When polygons are rolled up, the routine sums the positive and negative areas and consequently accounts for holes.
If the PolySet's projection attribute equals
"LL", the function projects the PolySet in UTM first.
If the PolySet's zone attribute exists, it uses it for
the conversion. Otherwise, it computes the mean longitude and uses
that value to determine the zone. The longitude range of zone
is .
PolyData with columns PID, SID (may be
missing), and area. If the projection equals "LL" or
"UTM", the units of area are square kilometres.
Nicholas M. Boers, Staff Software Engineer
Jobber, Edmonton AB
Last modified Rd: 2022-09-06
calcCentroid,
calcLength,
calcMidRange,
calcSummary,
locatePolys.
local(envir=.PBSmapEnv,expr={ #--- load the data (if using R) if (!is.null(version$language) && (version$language == "R")) data(nepacLL,envir=.PBSmapEnv) #--- convert LL to UTM so calculation makes sense attr(nepacLL, "zone") <- 9 nepacUTM <- convUL(nepacLL) #--- calculate and print the areas print(calcArea(nepacUTM)) })local(envir=.PBSmapEnv,expr={ #--- load the data (if using R) if (!is.null(version$language) && (version$language == "R")) data(nepacLL,envir=.PBSmapEnv) #--- convert LL to UTM so calculation makes sense attr(nepacLL, "zone") <- 9 nepacUTM <- convUL(nepacLL) #--- calculate and print the areas print(calcArea(nepacUTM)) })
Calculate the centroids of polygons found in a PolySet.
calcCentroid (polys, rollup = 3)calcCentroid (polys, rollup = 3)
polys |
PolySet to use. |
rollup |
level of detail in the results; |
If rollup equals 1, the results contain a centroid for
each unique PID only. When it equals 2, they contain
entries for outer contours only. Finally, setting it to 3
prevents roll-up, and they contain a centroid for each unique
(PID, SID).
PolyData with columns PID, SID (may be missing),
X, and Y.
Nicholas M. Boers, Staff Software Engineer
Jobber, Edmonton AB
Last modified Rd: 2013-04-10
calcArea,
calcLength,
calcMidRange,
calcSummary,
locateEvents,
locatePolys.
local(envir=.PBSmapEnv,expr={ #--- load the data (if using R) if (!is.null(version$language) && (version$language=="R")) data(nepacLL,envir=.PBSmapEnv) #--- calculate and print the centroids for several polygons print(calcCentroid(nepacLL[is.element(nepacLL$PID,c(33,39,47)),])) })local(envir=.PBSmapEnv,expr={ #--- load the data (if using R) if (!is.null(version$language) && (version$language=="R")) data(nepacLL,envir=.PBSmapEnv) #--- calculate and print the centroids for several polygons print(calcCentroid(nepacLL[is.element(nepacLL$PID,c(33,39,47)),])) })
Calculate the convex hull for a set of points.
calcConvexHull(xydata, keepExtra=FALSE) .closestPoint(pts, pt)calcConvexHull(xydata, keepExtra=FALSE) .closestPoint(pts, pt)
xydata |
|
keepExtra |
|
pts |
|
pt |
|
Uses the function chull() in the package grDevices.
By default, it ignores all columns other than X and Y;
however, the user can choose to retain additional columns in xydata
by specifying keepExtra=TRUE.
PolySet with columns PID, POS, X, Y,
and additional columns in xydata if keepExtra=TRUE.
Auxiliary dot function '.closestPoint' returns a vector of length
'pts' where TRUE indicates that the point is closest to 'pt'.
Returns several TRUE values when several points are equidistant.
Nicholas M. Boers, Software Engineer, Jobber, Edmonton AB
Maintainer: Rowan Haigh, Program Head – Offshore Rockfish
Pacific Biological Station (PBS), Fisheries & Oceans Canada (DFO), Nanaimo BC
locus opus: Remote office, Vancouver BC
Last modified Rd: 2024-09-03
In package PBSmapping:addPoints,
addPolys,
calcArea,
calcCentroid,
calcMidRange,
calcSummary,
locateEvents,
plotMap,
plotPoints,
plotPolys.
local(envir=.PBSmapEnv,expr={ oldpar = par(no.readonly=TRUE) data(surveyData,envir=.PBSmapEnv) #--- plot the convex hull, and then plot the points plotMap(calcConvexHull(surveyData),col="moccasin") addPoints(surveyData,col="blue",pch=17,cex=.6) par(oldpar) })local(envir=.PBSmapEnv,expr={ oldpar = par(no.readonly=TRUE) data(surveyData,envir=.PBSmapEnv) #--- plot the convex hull, and then plot the points plotMap(calcConvexHull(surveyData),col="moccasin") addPoints(surveyData,col="blue",pch=17,cex=.6) par(oldpar) })
Calculate the great-circle distance between geographic (LL) coordinates. Also calculate the initial bearing of the great-circle arc (at its starting point).
calcGCdist(lon1, lat1, lon2, lat2, R=6371.2)calcGCdist(lon1, lat1, lon2, lat2, R=6371.2)
lon1 |
|
lat1 |
|
lon2 |
|
lat2 |
|
R |
|
The great-circle distance is calculated between two points along a
spherical surface using the shortest distance and disregarding
topography.
Method 1: Haversine Formula
where = latitude (in radians), = longitude (in radians), = radius (km) of the Earth, = square of half the chord length between the points, = angular distance in radians, = great-circle distance (km) between two points.
Method 2: Spherical Law of Cosines
The initial bearing (aka forward azimuth) for the start point can be calculated using:
A list obect containing:a – Haversine = square of half the chord length between the points,c – Haversine = angular distance in radians,d – Haversine = great-circle distance (km) between two points,d2 – Law of Cosines = great-circle distance (km) between two points,theta – Initial bearing (degrees) for the start point.
If one uses the north geomagnetic pole as an end point,
crudely approximates the magnetic declination.
Rowan Haigh, Program Head – Offshore Rockfish
Pacific Biological Station (PBS), Fisheries & Oceans Canada (DFO), Nanaimo BC
locus opus: Remote office, Vancouver BC
Last modified Rd: 2024-03-11
Movable Type Scripts – Calculate distance, bearing and more between Latitude/Longitude points
In package PBSmapping:addCompass,
calcArea,
calcCentroid,
calcLength
local(envir=.PBSmapEnv,expr={ #-- Distance between southern BC waters and north geomagnetic pole print(calcGCdist(-126.5,48.6,-72.7,80.4)) })local(envir=.PBSmapEnv,expr={ #-- Distance between southern BC waters and north geomagnetic pole print(calcGCdist(-126.5,48.6,-72.7,80.4)) })
Calculate the length of polylines found in a PolySet.
calcLength (polys, rollup = 3, close = FALSE)calcLength (polys, rollup = 3, close = FALSE)
polys |
PolySet to use. |
rollup |
level of detail in the results; |
close |
Boolean value; if |
If rollup equals 1, the results contain an entry for
each unique PID only. Setting it to 3 prevents roll-up,
and they contain an entry for each unique (PID, SID).
If the projection attribute equals "LL", this routine uses
Great Circle distances to compute the surface length of each polyline.
In doing so, the algorithm simplifies Earth to a sphere.
If the projection attribute equals "UTM" or 1, this
routine uses Pythagoras' Theorem to calculate lengths.
PolyData with columns PID, SID (may be missing),
and length. If projection equals "UTM" or
"LL", lengths are in kilometres. Otherwise, lengths are in the
same unit as the input PolySet.
Nicholas M. Boers, Staff Software Engineer
Jobber, Edmonton AB
Last modified Rd: 2013-04-10
calcArea,
calcCentroid,
calcMidRange,
calcSummary,
locatePolys.
local(envir=.PBSmapEnv,expr={ #--- load the data (if using R) if (!is.null(version$language) && (version$language=="R")) data(nepacLL,envir=.PBSmapEnv) #--- calculate the perimeter of Vancouver Island print(calcLength(nepacLL[nepacLL$PID==33, ])) })local(envir=.PBSmapEnv,expr={ #--- load the data (if using R) if (!is.null(version$language) && (version$language=="R")) data(nepacLL,envir=.PBSmapEnv) #--- calculate the perimeter of Vancouver Island print(calcLength(nepacLL[nepacLL$PID==33, ])) })
Calculate the midpoint of the X/Y ranges of polygons
found in a PolySet.
calcMidRange (polys, rollup = 3)calcMidRange (polys, rollup = 3)
polys |
PolySet to use. |
rollup |
level of detail in the results; |
If rollup equals 1, the results contain a mean range for
each unique PID only. When it equals 2, they contain
entries for outer contours only. Finally, setting it to 3
prevents roll-up, and they contain a mean range for each unique
(PID, SID).
PolyData with columns PID, SID (may be missing),
X, and Y.
Nicholas M. Boers, Staff Software Engineer
Jobber, Edmonton AB
Last modified Rd: 2013-04-10
calcArea,
calcCentroid,
calcLength,
calcSummary.
local(envir=.PBSmapEnv,expr={ #--- load the data (if using R) if (!is.null(version$language) && (version$language=="R")) data(nepacLL,envir=.PBSmapEnv) #--- calculate and print the centroids for several polygons print(calcMidRange(nepacLL[is.element(nepacLL$PID,c(33,39,47)),])) })local(envir=.PBSmapEnv,expr={ #--- load the data (if using R) if (!is.null(version$language) && (version$language=="R")) data(nepacLL,envir=.PBSmapEnv) #--- calculate and print the centroids for several polygons print(calcMidRange(nepacLL[is.element(nepacLL$PID,c(33,39,47)),])) })
Apply functions to polygons in a PolySet.
calcSummary (polys, rollup = 3, FUN, ...)calcSummary (polys, rollup = 3, FUN, ...)
polys |
PolySet to use. |
rollup |
level of detail in the results; |
FUN |
the function to apply; it must accept a vector and return a vector or scalar. |
... |
optional arguments for |
If rollup equals 1, the results contain an entry for
each unique PID only. When it equals 2, they contain
entries for outer contours only. Finally, setting it to 3
prevents roll-up, and they contain an entry for each unique
(PID, SID).
PolyData with columns PID, SID (may be missing),
X, and Y. If FUN returns a vector of length
greater than 1 (say n), names the columns X1, X2,
..., Xn and Y1, Y2, ..., Yn.
Nicholas M. Boers, Staff Software Engineer
Jobber, Edmonton AB
Last modified Rd: 2013-04-10
calcArea,
calcCentroid,
calcConvexHull,
calcLength,
calcMidRange,
combineEvents,
findPolys,
locateEvents,
locatePolys,
makeGrid,
makeProps.
local(envir=.PBSmapEnv,expr={ #--- load the data (if using R) if (!is.null(version$language) && (version$language=="R")) data(nepacLL,envir=.PBSmapEnv) #--- calculate and print the centroids for several polygons print(calcSummary(nepacLL[is.element(nepacLL$PID,c(33,39,47)),], rollup=3, FUN=mean)) })local(envir=.PBSmapEnv,expr={ #--- load the data (if using R) if (!is.null(version$language) && (version$language=="R")) data(nepacLL,envir=.PBSmapEnv) #--- calculate and print the centroids for several polygons print(calcSummary(nepacLL[is.element(nepacLL$PID,c(33,39,47)),], rollup=3, FUN=mean)) })
Calculate the Voronoi (Dirichlet) tesselation for a set of points.
calcVoronoi(xydata, xlim=NULL, ylim=NULL, eps=1e-09, frac=0.0001) .expandEdges(polys, pts, xlim, ylim)calcVoronoi(xydata, xlim=NULL, ylim=NULL, eps=1e-09, frac=0.0001) .expandEdges(polys, pts, xlim, ylim)
xydata |
|
xlim |
|
ylim |
|
eps |
|
frac |
|
polys |
|
pts |
|
This routine ignores all columns other than 'X' and 'Y'.
If the user leaves 'xlim' and 'ylim' unspecified, the
function defaults to the range of the data with each extent expanded
by ten percent of the range.
This function sets the attribute 'projection' to 1 and the
attribute 'zone' to NULL as it assumes this projection in
its calculations.
'PolySet' with columns 'PID', 'POS', 'X', and 'Y'.
Auxiliary dot function '.expandEdges' returns an expanded 'PolySet'.
Nicholas M. Boers, Software Engineer, Jobber, Edmonton AB
Maintainer: Rowan Haigh, Program Head – Offshore Rockfish
Pacific Biological Station (PBS), Fisheries & Oceans Canada (DFO), Nanaimo BC
locus opus: Remote office, Vancouver BC
Last modified Rd: 2024-09-03
In package PBSmapping:addPoints,
addPolys,
calcArea,
calcCentroid,
calcConvexHull,
calcMidRange,
calcSummary,
locateEvents,
plotMap,
plotPoints,
plotPolys,
PolySet
local(envir=.PBSmapEnv,expr={ oldpar = par(no.readonly=TRUE) #--- create some EventData events <- as.EventData(data.frame( EID=1:200, X=rnorm(200), Y=rnorm(200)), projection=1) #--- calculate the Voronoi tesselation polys <- calcVoronoi(events) #--- create PolyData to color it based on area polyData <- calcArea(polys) names(polyData)[is.element(names(polyData), "area")] <- "Z" colSeq <- seq(0.4, 0.95, length=4) polyData <- makeProps(polyData, breaks=quantile(polyData$Z,c(0,.25,.5,.75,1)), propName="col", propVals=rgb(colSeq,colSeq,colSeq)) #--- plot the tesselation plotMap(polys, polyProps=polyData) #--- plot the points addPoints(events, pch=19) par(oldpar) })local(envir=.PBSmapEnv,expr={ oldpar = par(no.readonly=TRUE) #--- create some EventData events <- as.EventData(data.frame( EID=1:200, X=rnorm(200), Y=rnorm(200)), projection=1) #--- calculate the Voronoi tesselation polys <- calcVoronoi(events) #--- create PolyData to color it based on area polyData <- calcArea(polys) names(polyData)[is.element(names(polyData), "area")] <- "Z" colSeq <- seq(0.4, 0.95, length=4) polyData <- makeProps(polyData, breaks=quantile(polyData$Z,c(0,.25,.5,.75,1)), propName="col", propVals=rgb(colSeq,colSeq,colSeq)) #--- plot the tesselation plotMap(polys, polyProps=polyData) #--- plot the points addPoints(events, pch=19) par(oldpar) })
Clip a 'PolySet', where each unique ('PID', 'SID')
describes a polygon or polyline.
clipPolys (polys, xlim, ylim, keepExtra = FALSE) clipLines (polys, xlim, ylim, keepExtra = FALSE) .clip(polys, xlim, ylim, isPolygons, keepExtra)clipPolys (polys, xlim, ylim, keepExtra = FALSE) clipLines (polys, xlim, ylim, keepExtra = FALSE) .clip(polys, xlim, ylim, isPolygons, keepExtra)
polys |
|
xlim |
|
ylim |
|
keepExtra |
|
isPolygons |
|
For each discrete polygon, the function connects vertices 1 and N
(does not connect vertices 1 and N for discrete polylines).
It recalculates the 'POS' values for each vertex, saving the old
values in a column named 'oldPOS'. For new vertices, it sets
'oldPOS' to NA.
'PolySet' containing the input data, with some points added or
removed. A new column 'oldPOS' records the original 'POS'
value for each vertex.
Nicholas M. Boers, Software Engineer, Jobber, Edmonton AB
Maintainer: Rowan Haigh, Program Head – Offshore Rockfish
Pacific Biological Station (PBS), Fisheries & Oceans Canada (DFO), Nanaimo BC
locus opus: Remote office, Vancouver BC
Last modified Rd: 2024-09-03
In package PBSmapping:addLines,
addPolys,
addStipples,
PolySet
local(envir=.PBSmapEnv,expr={ oldpar = par(no.readonly=TRUE) #--- create a triangle that will be clipped polys <- data.frame(PID=rep(1, 3), POS=1:3, X=c(0,1,.5), Y=c(0,0,1)) #--- clip the triangle in the X direction, and plot the results plotPolys(clipPolys(polys,xlim=c(0,.75),ylim=range(polys[,"Y"])),col=2) par(oldpar) })local(envir=.PBSmapEnv,expr={ oldpar = par(no.readonly=TRUE) #--- create a triangle that will be clipped polys <- data.frame(PID=rep(1, 3), POS=1:3, X=c(0,1,.5), Y=c(0,0,1)) #--- clip the triangle in the X direction, and plot the results plotPolys(clipPolys(polys,xlim=c(0,.75),ylim=range(polys[,"Y"])),col=2) par(oldpar) })
Close a PolySet of polylines to form polygons.
closePolys (polys)closePolys (polys)
polys |
PolySet to close. |
Generally, run fixBound before this function. The ranges of a
PolySet's X and Y columns define the boundary.
For each discrete polygon, this function determines if the first and
last points lie on a boundary. If both points lie on the same
boundary, it adds no points. However, if they lie on different
boundaries, it may add one or two corners to the polygon.
When the boundaries are adjacent, one corner will be added as follows:
top boundary + left boundary implies add top-left corner;
top boundary + right boundary implies add top-right corner;
bottom boundary + left boundary implies add bottom-left corner;
bottom boundary + right boundary implies add bottom-right corner.
When the boundaries are opposite, it first adds the corner closest to a starting or ending polygon vertex. This determines a side (left-right or bottom-top) that connects the opposite boundaries. Then, it adds the other corner of that side to close the polygon.
PolySet identical to polys, except for possible
additional corner points.
Nicholas M. Boers, Staff Software Engineer
Jobber, Edmonton AB
Last modified Rd: 2013-04-10
local(envir=.PBSmapEnv,expr={ oldpar = par(no.readonly=TRUE) #--- 4 corners polys <- data.frame( PID = c(1, 1, 2, 2, 3, 3, 4, 4), POS = c(1, 2, 1, 2, 1, 2, 1, 2), X = c(0, 1, 2, 3, 0, 1, 2, 3), Y = c(1, 0, 0, 1, 2, 3, 3, 2)) plotPolys(closePolys(polys), col=2) #--- 2 corners and 1 opposite polys <- data.frame( PID = c(1, 1, 2, 2, 3, 3, 3), POS = c(1, 2, 1, 2, 1, 2, 3), X = c(0, 1, 0, 1, 5, 6, 1.5), Y = c(1, 0, 2, 3, 0, 1.5, 3)) plotPolys(closePolys(polys), col=2) par(oldpar) })local(envir=.PBSmapEnv,expr={ oldpar = par(no.readonly=TRUE) #--- 4 corners polys <- data.frame( PID = c(1, 1, 2, 2, 3, 3, 4, 4), POS = c(1, 2, 1, 2, 1, 2, 1, 2), X = c(0, 1, 2, 3, 0, 1, 2, 3), Y = c(1, 0, 0, 1, 2, 3, 3, 2)) plotPolys(closePolys(polys), col=2) #--- 2 corners and 1 opposite polys <- data.frame( PID = c(1, 1, 2, 2, 3, 3, 3), POS = c(1, 2, 1, 2, 1, 2, 3), X = c(0, 1, 0, 1, 5, 6, 1.5), Y = c(1, 0, 2, 3, 0, 1.5, 3)) plotPolys(closePolys(polys), col=2) par(oldpar) })
Combine measurements associated with events that occur in the same polygon.
combineEvents (events, locs, FUN, ..., bdryOK = TRUE)combineEvents (events, locs, FUN, ..., bdryOK = TRUE)
events |
EventData with at least four columns ( |
locs |
LocationSet usually resulting from a call to
|
FUN |
a function that produces a scalar from a vector
(e.g., |
... |
optional arguments for |
bdryOK |
Boolean value; if |
This function combines measurements associated with events that occur
in the same polygon. Each event (EID) has a corresponding
measurement Z. The locs data frame (usually output from
findPolys) places events within polygons. Thus, each
polygon (PID, SID) determines a set of events within it,
and a corresponding vector of measurements Zv. The function
returns FUN(Zv), a summary of measurements within each polygon.
PolyData with columns PID, SID (if in
locs), and Z.
Nicholas M. Boers, Staff Software Engineer
Jobber, Edmonton AB
Last modified Rd: 2013-04-10
findCells,
findPolys,
locateEvents,
locatePolys,
makeGrid,
makeProps.
local(envir=.PBSmapEnv,expr={ #--- create an EventData data frame: let each event have Z = 1 events <- data.frame(EID=1:10, X=1:10, Y=1:10, Z=rep(1, 10)) #--- example output from findPolys where 1 event occurred in the first #--- polygon, 3 in the second, and 6 in the third locs <- data.frame(EID=1:10,PID=c(rep(1,1),rep(2,3),rep(3,6)),Bdry=rep(0,10)) #--- sum the Z column of the events in each polygon, and print the result print(combineEvents(events=events, locs=locs, FUN=sum)) })local(envir=.PBSmapEnv,expr={ #--- create an EventData data frame: let each event have Z = 1 events <- data.frame(EID=1:10, X=1:10, Y=1:10, Z=rep(1, 10)) #--- example output from findPolys where 1 event occurred in the first #--- polygon, 3 in the second, and 6 in the third locs <- data.frame(EID=1:10,PID=c(rep(1,1),rep(2,3),rep(3,6)),Bdry=rep(0,10)) #--- sum the Z column of the events in each polygon, and print the result print(combineEvents(events=events, locs=locs, FUN=sum)) })
Combine several polygons into a single polygon by modifying the
PID and SID indices.
combinePolys (polys)combinePolys (polys)
polys |
PolySet with one or more polygons, each with possibly several components/holes. |
This function accepts a PolySet containing one or more polygons
(PIDs), each with one or more components or holes
(SIDs). The SID column need not exist in the input.
The function combines these polygons into a single polygon by simply
renumbering the PID and SID indices. The resulting
PolySet contains a single PID (with the value 1) and uses
the SID value to differentiate between polygons, their
components, and holes.
PolySet, possibly with the addition of an SID column if
it did not already exist. The function may also reorder columns such
that PID, SID, POS, X and Y appear
first, in that order.
Nicholas M. Boers, Staff Software Engineer
Jobber, Edmonton AB
Last modified Rd: 2007-06-06
Convert output from contourLines into a PolySet.
convCP (data, projection = NULL, zone = NULL)convCP (data, projection = NULL, zone = NULL)
data |
contour line data, often from the
|
projection |
optional |
zone |
optional |
data contains a list as described below. The
contourLines function create a list suitable for the
data argument.
A three-element list describes each contour. The named elements in
this list include the scalar level, the vector x, and
the vector y. Vectors x and y must have equal
lengths. A higher-level list (data) contains one or more of
these contours lists.
A list with two named elements PolySet and PolyData.
The PolySet element contains a PolySet representation of the
contour lines. The PolyData element links each contour line
(PID, SID) with a level.
Nicholas M. Boers, Staff Software Engineer
Jobber, Edmonton AB
Last modified Rd: 2013-04-10
contour,
contourLines,
convLP,
makeTopography.
local(envir=.PBSmapEnv,expr={ oldpar = par(no.readonly=TRUE) #--- create sample data for the contourLines() function x <- seq(-0.5, 0.8, length=50); y <- x z <- outer(x, y, FUN = function(x,y) { sin(2*pi*(x^2+y^2))^2; } ) data <- contourLines(x, y, z, levels=c(0.2, 0.8)) #--- pass that sample data into convCP() result <- convCP(data) #--- plot the result plotLines(result$PolySet, projection=1) print(result$PolyData) par(oldpar) })local(envir=.PBSmapEnv,expr={ oldpar = par(no.readonly=TRUE) #--- create sample data for the contourLines() function x <- seq(-0.5, 0.8, length=50); y <- x z <- outer(x, y, FUN = function(x,y) { sin(2*pi*(x^2+y^2))^2; } ) data <- contourLines(x, y, z, levels=c(0.2, 0.8)) #--- pass that sample data into convCP() result <- convCP(data) #--- plot the result plotLines(result$PolySet, projection=1) print(result$PolyData) par(oldpar) })
Convert EventData/PolyData into a PolySet.
convDP (data, xColumns, yColumns)convDP (data, xColumns, yColumns)
data |
|
xColumns |
vector of X-column names. |
yColumns |
vector of Y-column names. |
This function expects data to contain several X- and Y-columns.
For example, consider data with columns x1, y1,
x2, and y2. Suppose xColumns = c("x1", "x2") and
yColumns = c("y1", "y2"). The result will contain
nrow(data) polygons. Each one will have two vertices,
(x1, y1) and (x2, y2) and POS values 1 and 2,
respectively. If data includes an SID column, so will
the result.
If data contains an EID and not a PID column,
the function uses the EIDs as PIDs.
If data contains both PID and EID columns,
the function assumes it is PolyData and ignores the EID
column.
PolySet with the same PIDs as those given in data. If
data has an SID column, the result will include it.
Nicholas M. Boers, Staff Software Engineer
Jobber, Edmonton AB
Last modified Rd: 2013-04-10
local(envir=.PBSmapEnv,expr={ oldpar = par(no.readonly=TRUE) #--- create sample PolyData polyData <- data.frame(PID=c(1, 2, 3), x1=c(1, 3, 5), y1=c(1, 3, 2), x2=c(1, 4, 5), y2=c(2, 4, 1), x3=c(2, 4, 6), y3=c(2, 3, 1)) #--- print PolyData print(polyData) #--- make a PolySet from PolyData polys <- convDP(polyData, xColumns=c("x1", "x2", "x3"), yColumns=c("y1", "y2", "y3")) #--- print and plot the PolySet print(polys) plotLines(polys, xlim=c(0,7), ylim=c(0,5), col=2) par(oldpar) })local(envir=.PBSmapEnv,expr={ oldpar = par(no.readonly=TRUE) #--- create sample PolyData polyData <- data.frame(PID=c(1, 2, 3), x1=c(1, 3, 5), y1=c(1, 3, 2), x2=c(1, 4, 5), y2=c(2, 4, 1), x3=c(2, 4, 6), y3=c(2, 3, 1)) #--- print PolyData print(polyData) #--- make a PolySet from PolyData polys <- convDP(polyData, xColumns=c("x1", "x2", "x3"), yColumns=c("y1", "y2", "y3")) #--- print and plot the PolySet print(polys) plotLines(polys, xlim=c(0,7), ylim=c(0,5), col=2) par(oldpar) })
Convert two polylines into a polygon.
convLP (polyA, polyB, reverse = TRUE)convLP (polyA, polyB, reverse = TRUE)
polyA |
PolySet containing a polyline. |
polyB |
PolySet containing a polyline. |
reverse |
Boolean value; if |
The resulting PolySet contains all the vertices from
polyA in their original order. If reverse = TRUE, this
function appends the vertices from polyB in the reverse order
(nrow(polyB):1). Otherwise, it appends them in their original
order. The PID column equals the PID of polyA.
No SID column appears in the result. The resulting polygon is
an exterior boundary.
PolySet with a single PID that is the same as
polyA. The result contains all the vertices in polyA and
polyB. It has the same projection and zone
attributes as those in the input PolySets. If an input PolySet's
attributes equal NULL, the function uses the other
PolySet's. If the PolySet attributes conflict, the result's attribute
equals NULL.
Nicholas M. Boers, Staff Software Engineer
Jobber, Edmonton AB
Last modified Rd: 2013-04-10
addLines,
appendPolys,
closePolys,
convCP,
joinPolys,
plotLines.
local(envir=.PBSmapEnv,expr={ oldpar = par(no.readonly=TRUE) #--- create two polylines polyline1 <- data.frame(PID=rep(1,2),POS=1:2,X=c(1,4),Y=c(1,4)) polyline2 <- data.frame(PID=rep(1,2),POS=1:2,X=c(2,5),Y=c(1,4)) #--- create two plots to demonstrate the effect of `reverse' par(mfrow=c(2, 1)) plotPolys(convLP(polyline1, polyline2, reverse=TRUE), col=2) plotPolys(convLP(polyline1, polyline2, reverse=FALSE), col=3) par(oldpar) })local(envir=.PBSmapEnv,expr={ oldpar = par(no.readonly=TRUE) #--- create two polylines polyline1 <- data.frame(PID=rep(1,2),POS=1:2,X=c(1,4),Y=c(1,4)) polyline2 <- data.frame(PID=rep(1,2),POS=1:2,X=c(2,5),Y=c(1,4)) #--- create two plots to demonstrate the effect of `reverse' par(mfrow=c(2, 1)) plotPolys(convLP(polyline1, polyline2, reverse=TRUE), col=2) plotPolys(convLP(polyline1, polyline2, reverse=FALSE), col=3) par(oldpar) })
Convert coordinates between UTM and Lon/Lat.
convUL (xydata, km=TRUE, southern=NULL)convUL (xydata, km=TRUE, southern=NULL)
xydata |
|
km |
|
southern |
|
The object xydata must possess a projection attribute that
identifies the current projection. If the data frame contains UTM
coordinates, it must also have a zone attribute equal to a
number between 1 and 60 (inclusive). If it contains
geographic (longitude/latitude) coordinates and the zone attribute is
missing, the function computes the mean longitude and uses that value
to determine the zone. The longitude range of zone is
.
This function converts the X and Y columns of
xydata from "LL" to "UTM" or vice-versa. If the
data span more than one zone to the right or left of the intended
central zone, the underlying algorithm may produce erroneous
results. This limitation means that the user should use
the most central zone of the mapped region, or allow the function to determine
the central zone when converting from geographic to UTM coordinates.
After the conversion, this routine adjusts the data frame's attributes accordingly.
A data frame identical to xydata, except that the X and
Y columns contain the results of the conversion, and the
projection attribute matches the new projection.
Nicholas M. Boers, Software Engineer, Jobber, Edmonton AB
Maintainer: Rowan Haigh, Program Head – Offshore Rockfish
Pacific Biological Station (PBS), Fisheries & Oceans Canada (DFO), Nanaimo BC
locus opus: Regional Headquarters, Vancouver BC
Last modified Rd: 2022-09-06
Ordnance Survey. (2020) A guide to coordinate systems in Great Britain. Copyright Ordnance Survey 2018 (v3.6). Southampton, UK.
local(envir=.PBSmapEnv,expr={ oldpar = par(no.readonly=TRUE) #--- load the data data(nepacLL,envir=.PBSmapEnv) #--- set the zone attribute #--- use a zone that is most central to the mapped region attr(nepacLL, "zone") <- 6 #--- convert and plot the result nepacUTM <- convUL(nepacLL) plotMap(nepacUTM) par(oldpar) })local(envir=.PBSmapEnv,expr={ oldpar = par(no.readonly=TRUE) #--- load the data data(nepacLL,envir=.PBSmapEnv) #--- set the zone attribute #--- use a zone that is most central to the mapped region attr(nepacLL, "zone") <- 6 #--- convert and plot the result nepacUTM <- convUL(nepacLL) plotMap(nepacUTM) par(oldpar) })
Divide a single polygon (with several outer-contour components) into
several polygons, a polygon for each outer contour, by modifying the
PID and SID indices.
dividePolys (polys)dividePolys (polys)
polys |
PolySet with one or more polygons, each with possibly several components/holes. |
Given the input PolySet, this function renumbers the PID
and SID indices so that each outer contour has a unique PID and
is followed by all of its holes, identifying them with SIDs
greater than one.
PolySet, possibly with the addition of an SID column if
it did not already exist. The function may also reorder columns such
that PID, SID, POS, X and Y appear
first, in that order.
Nicholas M. Boers, Staff Software Engineer
Jobber, Edmonton AB
Last modified Rd: 2007-06-06
Routines to add various items to existing plots.
.addAxis(xlim, ylim, tckLab, tck, tckMinor, ...) .addAxis2(side=1:2, xlim, ylim, tckLab, tck, tckMinor, ...) .addBubblesLegend(radii.leg, usr.xdiff, usr.ydiff, symbol.zero, symbol.fg, symbol.bg, legend.pos, legend.breaks, legend.type, legend.title, legend.cex, ...) .addCorners(polys, ptSummary) .addFeature(feature, data, polyProps, isEventData, cex=NULL, col=NULL, font=NULL, pch=NULL, ...) .addLabels(projection=NULL, ...) .addProps(type, polyProps, ...).addAxis(xlim, ylim, tckLab, tck, tckMinor, ...) .addAxis2(side=1:2, xlim, ylim, tckLab, tck, tckMinor, ...) .addBubblesLegend(radii.leg, usr.xdiff, usr.ydiff, symbol.zero, symbol.fg, symbol.bg, legend.pos, legend.breaks, legend.type, legend.title, legend.cex, ...) .addCorners(polys, ptSummary) .addFeature(feature, data, polyProps, isEventData, cex=NULL, col=NULL, font=NULL, pch=NULL, ...) .addLabels(projection=NULL, ...) .addProps(type, polyProps, ...)
xlim |
|
ylim |
|
tckLab |
|
tck |
|
tckMinor |
|
... |
|
side |
|
radii.leg |
|
usr.xdiff |
|
usr.ydiff |
|
symbol.zero |
|
symbol.fg |
|
symbol.bg |
|
legend.pos |
|
legend.breaks |
|
legend.type |
|
legend.title |
|
legend.cex |
|
polys |
|
ptSummary |
|
feature |
|
data |
|
polyProps |
|
isEventData |
|
cex |
|
col |
|
font |
|
pch |
|
projection |
|
type |
|
Internal functions add features to an existing map plot.
Nothing in particular.
Nicholas M. Boers, Software Engineer, Jobber, Edmonton AB
Maintainer: Rowan Haigh, Program Head – Offshore Rockfish
Pacific Biological Station (PBS), Fisheries & Oceans Canada (DFO), Nanaimo BC
locus opus: Remote office, Vancouver BC
Last modified Rd: 2024-09-03
In package PBSmapping:plotMap,
addBubbles,
addPolys,
addLabels,
addStipples
Calculate distance and/or orientation.
.calcDist(polys) .calcOrientation(polys).calcDist(polys) .calcOrientation(polys)
polys |
|
.calcdist : Equatorial radius 6,378.14 km; Polar radius 6,356.78 km; Mean radius 6,371.3 km
.calcOrientation : Calls C code 'calcOrientation' using .C().
.calcdist : distance vector (distances between each point)
.calcOrientation : data frame with 'orientation' column
(-1 when counter-clockwise; 0 when N/A; +1 when clockwise) or
NULL if no rows in output
Nicholas M. Boers, Software Engineer, Jobber, Edmonton AB
Maintainer: Rowan Haigh, Program Head – Offshore Rockfish
Pacific Biological Station (PBS), Fisheries & Oceans Canada (DFO), Nanaimo BC
locus opus: Remote office, Vancouver BC
Last modified Rd: 2024-09-03
In package PBSmapping:calcLength,
placeHoles,
thickenPolys
Auxiliary (dot) functions to check various attributes of PBmapping objects.
.checkClipLimits(limits) .checkProjection(projectionPlot, projectionPoly) .checkRDeps(caller="unspecified", requires=NULL).checkClipLimits(limits) .checkProjection(projectionPlot, projectionPoly) .checkRDeps(caller="unspecified", requires=NULL)
limits |
|
projectionPlot |
|
projectionPoly |
|
caller |
|
requires |
|
Auxiliary dot functions to facilitate the machinations of PBSmapping.
.checkClipLimits : nada; just a limit checker.checkProjection : nada; just a projection checker.checkRDeps : nada; just a package checker
Nicholas M. Boers, Software Engineer, Jobber, Edmonton AB
Maintainer: Rowan Haigh, Program Head – Offshore Rockfish
Pacific Biological Station (PBS), Fisheries & Oceans Canada (DFO), Nanaimo BC
locus opus: Remote office, Vancouver BC
Last modified Rd: 2024-09-03
In package PBSmapping:.checkClipLimits used by:importGSHHS.checkProjection used by:addLabels,
addLines,
addPoints,
addPolys,
addStipples.checkRDeps used by:calcVoronoi,
importShapefile (temporarily unavailable)
Create indices for indexing map data structures (PolySets).
.createIDs(x, cols, fastIDdig=NULL) .createFastIDdig(polysA, polysB=NULL, cols) .createGridIDs(d, addSID, byrow).createIDs(x, cols, fastIDdig=NULL) .createFastIDdig(polysA, polysB=NULL, cols) .createGridIDs(d, addSID, byrow)
x |
|
cols |
|
fastIDdig |
|
polysA |
|
polysB |
|
d |
|
addSID |
|
byrow |
|
.createIDs : create IDs (or IDX) column from its input.
.createFastIDdig : determine the maximum number of digits in
the second column of a data frame.
If given two data frames ('polysA' and 'polysB'),
determine the maximum between the two data frames.
.createGridIDs : Create IDs for a grid according to the
addSID and byrow arguments.
.createIDs : a vector of integer or real-number indices.
.createFastIDdig : maximum number of digits to use in real-number index.
.createGridIDs : a modified grid PolySet.
Nicholas M. Boers, Software Engineer, Jobber, Edmonton AB
Maintainer: Rowan Haigh, Program Head – Offshore Rockfish
Pacific Biological Station (PBS), Fisheries & Oceans Canada (DFO), Nanaimo BC
locus opus: Remote office, Vancouver BC
Last modified Rd: 2024-09-25
In package PBSmapping:plotMap,
addBubbles,
addPolys,
addLabels,
addStipples
Fix things like GSHHSWorld.
.fixGSHHSWorld(world).fixGSHHSWorld(world)
world |
|
Determine PID of Antarctica, which is used to extract the current Antarctica. The continent is rebuilt from left to right, creating a new, very wide Antarctica. It is then clipped and merged into the existing world (with other polygons).
A revised PolySet of the world with a peeled version of Antarctica.
This function is not called by any PBSmapping functions.
It can be used after 'importGSHHS' uses ylim=c(-90,90).
Nicholas M. Boers, Software Engineer, Jobber, Edmonton AB
Maintainer: Rowan Haigh, Program Head – Offshore Rockfish
Pacific Biological Station (PBS), Fisheries & Oceans Canada (DFO), Nanaimo BC
locus opus: Remote office, Vancouver BC
Last modified Rd: 2024-09-03
In package PBSmapping:importGSHHS,
clipPolys
## Not run: worldLL = importGSHHS("gshhs_l.b",xlim=c(-20,360),ylim=c(-90,90),level=1,n=15,xoff=0) worldLL = .fixGSHHSWorld(worldLL) ## End(Not run)## Not run: worldLL = importGSHHS("gshhs_l.b",xlim=c(-20,360),ylim=c(-90,90),level=1,n=15,xoff=0) worldLL = .fixGSHHSWorld(worldLL) ## End(Not run)
Routines to get various attributes from PBmapping objects.
.getBasename(fn, ext) .getGridPars(polys, fullValidation=TRUE).getBasename(fn, ext) .getGridPars(polys, fullValidation=TRUE)
fn |
|
ext |
|
polys |
|
fullValidation |
|
Auxiliary dot functions to facilitate the machinations of PBSmapping.
.getBasename : base name of a shapefile set (usually with multiple extensions)..getGridPars : PBSmapping grid object attributes.
Nicholas M. Boers, Software Engineer, Jobber, Edmonton AB
Maintainer: Rowan Haigh, Program Head – Offshore Rockfish
Pacific Biological Station (PBS), Fisheries & Oceans Canada (DFO), Nanaimo BC
locus opus: Remote office, Vancouver BC
Last modified Rd: 2024-09-03
In package PBSmapping:importShapefile (temporarily unavailable),
makeGrid,
is.PolySet
Determines which points (EventData) are located inside a polygon (PolySet).
.is.in(events, polys, use.names=TRUE).is.in(events, polys, use.names=TRUE)
events |
|
polys |
|
use.names |
|
Taps into the PBSmapping C code for 'findPolys'.
Reports events inside the polygon, outside the polygon, and whether all events are inside, outside, or on the boundary of the polygon.
Rowan Haigh, Program Head – Offshore Rockfish
Pacific Biological Station (PBS), Fisheries & Oceans Canada (DFO), Nanaimo BC
locus opus: Remote office, Vancouver BC
Last modified Rd: 2024-03-11
In PBSmapping:findPolys,
as.EventData,
is.PolySet,
findPolys
Convert PBS matrices to data frames, preserving certain attributes.
.mat2df(data).mat2df(data)
data |
|
Converts a PBS matrix object to a PBS data frame.
Preserves PBS class objects 'EventData', 'LocationSet',
'PolyData', 'PolySet', if any of them exist in the matrix.
Transformed data object (from matrix to data frame).
Nicholas M. Boers, Software Engineer, Jobber, Edmonton AB
Maintainer: Rowan Haigh, Program Head – Offshore Rockfish
Pacific Biological Station (PBS), Fisheries & Oceans Canada (DFO), Nanaimo BC
locus opus: Remote office, Vancouver BC
Last modified Rd: 2024-09-03
In package PBSmapping:addPoints,
.addFeature,
.validateData
Run built-in example figures for PBSmapping.
.PBSfigs(nfigs=1:10, wait=TRUE) .PBSfig01() .PBSfig02() .PBSfig03() .PBSfig04() .PBSfig05() .PBSfig06() .PBSfig07() .PBSfig08() .PBSfig09() .PBSfig10().PBSfigs(nfigs=1:10, wait=TRUE) .PBSfig01() .PBSfig02() .PBSfig03() .PBSfig04() .PBSfig05() .PBSfig06() .PBSfig07() .PBSfig08() .PBSfig09() .PBSfig10()
nfigs |
|
wait |
|
The wrapper function '.PBSfigs' displays pre-defined figures
from 1 to 10, depending on the numbers that the user specifies.
User can also call each figure directly, e.g., '.PBSfig05()'.
NULL
Nicholas M. Boers, Software Engineer, Jobber, Edmonton AB
Maintainer: Rowan Haigh, Program Head – Offshore Rockfish
Pacific Biological Station (PBS), Fisheries & Oceans Canada (DFO), Nanaimo BC
locus opus: Offsite, Vancouver BC
Last modified Rd: 2023-10-30
.PBSfigs(nfigs=1:3, wait=FALSE).PBSfigs(nfigs=1:3, wait=FALSE)
An environment set aside for PBSmapping.
.PBSmapEnv.PBSmapEnv
A new environment with a .GlobalEnv parent.
The environment is created in 'zzz.r' and can be
used by PBSmodelling functions 'lisp', 'tget',
'tput', 'tprint', and 'tcall', if this package has be loaded.
Currently, '.PBSmapEnv' is only used in '.validateData'
(in the namespace), and by a set of print functions specific to PBSmapping.
Generated by a call to the base function new.env().
In PBSmapping:print.EventData,
print.LocationSet,
print.PolyData,
print.PolySet,
print.summary.PBS
In PBSmodelling:lisp,
tget,
packList
Routines to validate data structures (e.g., required field names, attributes) that are passed to package functions.
.validateData(data, className, requiredCols=NULL, requiredAttr=NULL, noFactorCols=NULL, noNACols=NULL, keyCols=NULL, numericCols=NULL) .validateEventData(EventData) .validateLocationSet(LocationSet) .validatePolyData(PolyData) .validatePolyProps(polyProps, parCols=NULL) .validatePolySet(polys) .validateXYData(xyData).validateData(data, className, requiredCols=NULL, requiredAttr=NULL, noFactorCols=NULL, noNACols=NULL, keyCols=NULL, numericCols=NULL) .validateEventData(EventData) .validateLocationSet(LocationSet) .validatePolyData(PolyData) .validatePolyProps(polyProps, parCols=NULL) .validatePolySet(polys) .validateXYData(xyData)
data |
|
className |
|
requiredCols |
|
requiredAttr |
|
noFactorCols |
|
noNACols |
|
keyCols |
|
numericCols |
|
EventData |
|
LocationSet |
|
PolyData |
|
polyProps |
|
parCols |
|
polys |
|
xyData |
|
Internal functions check the validity of data objects used by PBSmapping.
The primary function is '.validateData'; other dot validate functions
are wrappers for the four main data structures:
'PolySet', 'PolyData', 'EventData', and 'LocationSet'.
The data object that is validated.
Nicholas M. Boers, Software Engineer, Jobber, Edmonton AB
Maintainer: Rowan Haigh, Program Head – Offshore Rockfish
Pacific Biological Station (PBS), Fisheries & Oceans Canada (DFO), Nanaimo BC
locus opus: Remote office, Vancouver BC
Last modified Rd: 2024-09-03
In package PBSmapping:PolySet,
PolyData,
EventData,
LocationSet
An EventData object comprises a data frame with at least three fields named
EID, X, and Y; each row specifies an event that occurs
at a specific point.
PBSmapping functions that expect EventData will accept properly formatted data frames in their place (see 'Details').
as.EventData attempts to coerce a data frame to an object with
class EventData.
is.EventData returns TRUE if its argument is of class
EventData.
as.EventData(x, projection = NULL, zone = NULL) is.EventData(x, fullValidation = TRUE)as.EventData(x, projection = NULL, zone = NULL) is.EventData(x, fullValidation = TRUE)
x |
data frame to be coerced or tested. |
projection |
optional |
zone |
optional |
fullValidation |
Boolean value; if |
Conceptually, an EventData object describes events (EID) that take place at specific points (X,Y) in two-dimensional space. Additional fields can specify measurements associated with these events. In a fishery context, EventData could describe fishing events associated with trawl tows, based on the fields:
EID - fishing event (tow) identification number;
X, Y - fishing location;
Duration - length of time for the tow;
Depth - average depth of the tow;
Catch - biomass captured.
Like PolyData, EventData can have attributes projection
and zone, which may be absent. Inserting the string
"EventData" as the class attribute's first element alters the
behaviour of some functions, including print (if
PBSprint is TRUE) and summary.
The as.EventData method returns an object with classes
"EventData" and "data.frame", in that order.
Nicholas M. Boers, Staff Software Engineer
Jobber, Edmonton AB
Last modified Rd: 2015-04-23
PolySet, PolyData, LocationSet
Extract PolyData from a PolySet. Columns for the
PolyData include those other than PID, SID,
POS, oldPOS, X, and Y.
extractPolyData (polys)extractPolyData (polys)
polys |
PolySet to use. |
This function identifies the PolySet's extra columns and
determines if those columns contain unique values for each
(PID, SID). Where they do, the (PID, SID)
will appear in the PolyData output with that unique value.
Where they do not, the extra column will contain NAs for that
(PID, SID).
PolyData with columns PID, SID, and any extra
columns.
Nicholas M. Boers, Staff Software Engineer
Jobber, Edmonton AB
Last modified Rd: 2013-04-10
local(envir=.PBSmapEnv,expr={ #--- create a PolySet with an extra column polys <- data.frame(PID = c(rep(1, 10), rep(2, 10)), POS = c(1:10, 1:10), X = c(rep(1, 10), rep(1, 10)), Y = c(rep(1, 10), rep(1, 10)), colour = (c(rep("green", 10), rep("red", 10)))) #--- extract the PolyData print(extractPolyData(polys)) })local(envir=.PBSmapEnv,expr={ #--- create a PolySet with an extra column polys <- data.frame(PID = c(rep(1, 10), rep(2, 10)), POS = c(1:10, 1:10), X = c(rep(1, 10), rep(1, 10)), Y = c(rep(1, 10), rep(1, 10)), colour = (c(rep("green", 10), rep("red", 10)))) #--- extract the PolyData print(extractPolyData(polys)) })
Find the grid cells in a PolySet that contain events specified in
EventData. Similar to findPolys, except this
function requires a PolySet resulting from
makeGrid. This restriction allows this function to
calculate the result with greater efficiency.
findCells (events, polys, includeBdry=NULL)findCells (events, polys, includeBdry=NULL)
events |
EventData to use. |
polys |
PolySet to use. |
includeBdry |
numeric: determines how points on boundaries are handled: |
The resulting data frame, a LocationSet, contains the columns
EID, PID, SID (if in polys), and
Bdry, where an event (EID) occurs in a polygon
(PID, SID). The Boolean (0,1) variable Bdry indicates
whether an event lies on a polygon's edge. Note that if an event lies
properly outside of all the polygons, then a record with (EID,
PID, SID) does not occur in the output. It may happen,
however, that an event occurs in multiple polygons (i.e., on two or
more boundaries). Thus, the same EID can occur more than once
in the output.
If an event happens to lie at the boundary intersection of four (or two) grid cells
then one EID will be associated with four (or two) grid cells. A user
can choose to manipulate this result by setting the argument includeBdry
to a numeric value that constrains the association of a boundary event to
0 or 1 grid cell (see argument description above).
LocationSet that links events with polygons.
Nicholas M. Boers, Staff Software Engineer
Jobber, Edmonton AB
Last modified Rd: 2014-12-15
findPolys,
makeGrid,
combineEvents,
locateEvents,
locatePolys,
LocationSet.
local(envir=.PBSmapEnv,expr={ oldpar = par(no.readonly=TRUE) #--- create some EventData: points in a diagonal line events <- data.frame(EID=1:11, X=seq(0, 2, length=11), Y=seq(0, 2, length=11)) events <- as.EventData(events, projection=1); #--- create a PolySet (a grid) polys <- makeGrid (x=seq(0, 2, by=0.50), y=seq(0, 2, by=0.50), projection=1) #--- show a picture plotPolys(polys, xlim=range(polys$X)+c(-0.1, 0.1), ylim=range(polys$Y)+c(-0.1, 0.1), projection=1) addPoints(events, col=2) #--- run findCells and print the results fc <- findCells(events, polys) fc <- fc[order(fc$EID, fc$PID, fc$SID), ] fc$label <- paste(fc$PID, fc$SID, sep=", ") print (fc) #--- add labels to the graph addLabels(as.PolyData(fc[!duplicated(paste(fc$PID,fc$SID)), ], projection=1), placement="CENTROID", polys=as.PolySet(polys, projection=1), col=4) par(oldpar) })local(envir=.PBSmapEnv,expr={ oldpar = par(no.readonly=TRUE) #--- create some EventData: points in a diagonal line events <- data.frame(EID=1:11, X=seq(0, 2, length=11), Y=seq(0, 2, length=11)) events <- as.EventData(events, projection=1); #--- create a PolySet (a grid) polys <- makeGrid (x=seq(0, 2, by=0.50), y=seq(0, 2, by=0.50), projection=1) #--- show a picture plotPolys(polys, xlim=range(polys$X)+c(-0.1, 0.1), ylim=range(polys$Y)+c(-0.1, 0.1), projection=1) addPoints(events, col=2) #--- run findCells and print the results fc <- findCells(events, polys) fc <- fc[order(fc$EID, fc$PID, fc$SID), ] fc$label <- paste(fc$PID, fc$SID, sep=", ") print (fc) #--- add labels to the graph addLabels(as.PolyData(fc[!duplicated(paste(fc$PID,fc$SID)), ], projection=1), placement="CENTROID", polys=as.PolySet(polys, projection=1), col=4) par(oldpar) })
Find the polygons in a PolySet that contain events specified in EventData.
findPolys (events, polys, maxRows = 1e+05, includeBdry=NULL)findPolys (events, polys, maxRows = 1e+05, includeBdry=NULL)
events |
EventData to use. |
polys |
PolySet to use. |
maxRows |
estimated maximum number of rows in the output LocationSet. |
includeBdry |
numeric: determines how points on boundaries are handled: |
The resulting data frame, a LocationSet, contains the columns
EID, PID, SID (if in polys), and
Bdry, where an event (EID) occurs in a polygon
(PID, SID) and SID does not correspond to an
inner boundary. The Boolean variable Bdry indicates whether an
event lies on a polygon's edge. Note that if an event lies properly
outside of all the polygons, then a record with (EID,
PID, SID) does not occur in the output. It may happen,
however, that an event occurs in multiple polygons. Thus, the same
EID can occur more than once in the output.
If an event happens to lie at the boundary intersection of two or more polygons
then one EID will be associated with two or more polygons. A user
can choose to manipulate this result by setting the argument includeBdry
to a numeric value that constrains the association of a boundary event to
0 or 1 polygon (see argument description above).
LocationSet that links events with polygons.
Nicholas M. Boers, Staff Software Engineer
Jobber, Edmonton AB
Last modified Rd: 2014-12-15
combineEvents,
findCells,
locateEvents,
locatePolys,
LocationSet,
makeGrid.
local(envir=.PBSmapEnv,expr={ oldpar = par(no.readonly=TRUE) #--- create some EventData: a column of points at X = 0.5 events <- data.frame(EID=1:10, X=.5, Y=seq(0, 2, length=10)) events <- as.EventData(events, projection=1) #--- create a PolySet: two squares with the second above the first polys <- data.frame(PID=c(rep(1, 4), rep(2, 4)), POS=c(1:4, 1:4), X=c(0, 1, 1, 0, 0, 1, 1, 0), Y=c(0, 0, 1, 1, 1, 1, 2, 2)) polys <- as.PolySet(polys, projection=1) #--- show a picture plotPolys(polys, xlim=range(polys$X)+c(-0.1, 0.1), ylim=range(polys$Y)+c(-0.1, 0.1), projection=1); addPoints(events, col=2); #--- run findPolys and print the results print(findPolys(events, polys)) par(oldpar) })local(envir=.PBSmapEnv,expr={ oldpar = par(no.readonly=TRUE) #--- create some EventData: a column of points at X = 0.5 events <- data.frame(EID=1:10, X=.5, Y=seq(0, 2, length=10)) events <- as.EventData(events, projection=1) #--- create a PolySet: two squares with the second above the first polys <- data.frame(PID=c(rep(1, 4), rep(2, 4)), POS=c(1:4, 1:4), X=c(0, 1, 1, 0, 0, 1, 1, 0), Y=c(0, 0, 1, 1, 1, 1, 2, 2)) polys <- as.PolySet(polys, projection=1) #--- show a picture plotPolys(polys, xlim=range(polys$X)+c(-0.1, 0.1), ylim=range(polys$Y)+c(-0.1, 0.1), projection=1); addPoints(events, col=2); #--- run findPolys and print the results print(findPolys(events, polys)) par(oldpar) })
The ranges of a PolySet's X and Y columns define
its boundary. This function fixes a PolySet's vertices by
moving vertices near a boundary to the actual boundary.
fixBound (polys, tol)fixBound (polys, tol)
polys |
PolySet to fix. |
tol |
vector (length 1 or 2) specifying a percentage of
the ranges to use in defining near to a boundary. If
|
When moving vertices to a boundary, the function moves them strictly horizontally or vertically, as appropriate.
PolySet identical to the input, except for possible changes in
the X and Y columns.
Nicholas M. Boers, Staff Software Engineer
Jobber, Edmonton AB
Last modified Rd: 2013-04-10
closePolys,
fixPOS,
isConvex,
isIntersecting,
PolySet.
local(envir=.PBSmapEnv,expr={ oldpar = par(no.readonly=TRUE) #--- set up a long horizontal and long vertical line to extend the plot's #--- limits, and then try fixing the bounds of a line in the top-left #--- corner and a line in the bottom-right corner polys <- data.frame(PID=c(1, 1, 2, 2, 3, 3, 4, 4), POS=c(1, 2, 1, 2, 1, 2, 1, 2), X = c(0, 10, 5, 5, 0.1, 4.9, 5.1, 9.9), Y = c(5, 5, 0, 10, 5.1, 9.9, 0.1, 4.9)) polys <- fixBound(polys, tol=0.0100001) plotLines(polys) par(oldpar) })local(envir=.PBSmapEnv,expr={ oldpar = par(no.readonly=TRUE) #--- set up a long horizontal and long vertical line to extend the plot's #--- limits, and then try fixing the bounds of a line in the top-left #--- corner and a line in the bottom-right corner polys <- data.frame(PID=c(1, 1, 2, 2, 3, 3, 4, 4), POS=c(1, 2, 1, 2, 1, 2, 1, 2), X = c(0, 10, 5, 5, 0.1, 4.9, 5.1, 9.9), Y = c(5, 5, 0, 10, 5.1, 9.9, 0.1, 4.9)) polys <- fixBound(polys, tol=0.0100001) plotLines(polys) par(oldpar) })
Fix the POS column of a PolySet by recalculating it
using sequential integers.
fixPOS (polys, exteriorCCW = NA)fixPOS (polys, exteriorCCW = NA)
polys |
PolySet to fix. |
exteriorCCW |
Boolean value; if |
This function recalculates the POS values of each (PID,
SID) as either 1 to N or N to 1, depending on the order of
POS (ascending or descending) in the input data. POS
values in the input must be properly ordered (ascending or
descending), but they may contain fractional values. For example,
POS = 2.5 might correspond to a point manually added between
POS = 2 and POS = 3. If exteriorCCW = NA, all
other columns remain unchanged. Otherwise, it orders the X and
Y columns according to exteriorCCW.
PolySet with the same columns as the input, except for possible
changes to the POS, X, and Y columns.
Nicholas M. Boers, Staff Software Engineer
Jobber, Edmonton AB
Last modified Rd: 2013-04-10
closePolys,
fixBound,
isConvex,
isIntersecting,
PolySet.
local(envir=.PBSmapEnv,expr={ #--- create a PolySet with broken POS numbering polys <- data.frame(PID = c(rep(1, 10), rep(2, 10)), POS = c(seq(2, 10, length = 10), seq(10, 2, length = 10)), X = c(rep(1, 10), rep(1, 10)), Y = c(rep(1, 10), rep(1, 10))) #--- fix the POS numbering polys <- fixPOS(polys) #--- print the results print(polys) })local(envir=.PBSmapEnv,expr={ #--- create a PolySet with broken POS numbering polys <- data.frame(PID = c(rep(1, 10), rep(2, 10)), POS = c(seq(2, 10, length = 10), seq(10, 2, length = 10)), X = c(rep(1, 10), rep(1, 10)), Y = c(rep(1, 10), rep(1, 10))) #--- fix the POS numbering polys <- fixPOS(polys) #--- print the results print(polys) })
Import a text file and convert into EventData.
importEvents(EventData, projection=NULL, zone=NULL)importEvents(EventData, projection=NULL, zone=NULL)
EventData |
|
projection |
|
zone |
|
An imported EventData.
Nicholas M. Boers, Staff Software Engineer
Jobber, Edmonton AB
Last modified Rd: 2023-10-30
importPolys, importLocs, importGSHHS
Import data from a GSHHS database and convert data into a PolySet
with a PolyData attribute.
The database was originally called ‘Global Self-consistent, Hierarchical, High-resolution Shoreline’
(GSHHS, Wessel and Smith 1996), but ‘Shoreline’ was subsequently expanded to include more ‘Geography’ (GSHHG).
importGSHHS(gshhsDB, xlim, ylim, maxLevel=4, n=0, useWest=FALSE)importGSHHS(gshhsDB, xlim, ylim, maxLevel=4, n=0, useWest=FALSE)
gshhsDB |
|
xlim |
|
ylim |
|
maxLevel |
|
n |
|
useWest |
|
This routine requires a binary GSHHG (Global Self-consistent, Hierarchical,
High-resolution Geography) database file.
The GSHHG database has been released in the public domain.
(www.soest.hawaii.edu/pwessel/gshhg)
At the time of writing, the most recent binary database was the archive
file called gshhg-bin-2.3.7.zip.
The archive contains multiple binary files that contain geographical
coordinates for shorelines (gshhs), rivers (wdb_rivers), and
borders (wdb_borders).
The latter two come from World DataBank II (WDBII).
The five resolutions available are:
full (f), high (h), intermediate (i), low (l), and coarse (c).
This routine returns a PolySet object with an associated
PolyData attribute. The attribute contains four fields: (a) PID,
(b) SID, (c) Level, and (d) Source.
Each record corresponds to a line/polygon in the PolySet. The
Level indicates the line's/polygon's level (1=land, 2=lake,
3=island, 4=pond). The Source identifies the data source
(1=WVS, 0=CIA (WDBII)).
A PolySet with a PolyData attribute.
The function calls a C routine, also called importGSHHS, which returns
a set of map coordinates that is not always predictably laid out.
This issue stems from how the world is divided at the Greenwich meridian and
at the International Date Line.
The unpredictability occurs when user-specified X-limits span either of the
longitudinal meridians – (0, 360) or
(-180, 180).
This version of the R function attempts to stitch together the overlapping
edges of gshhs that run from -20 to
360 (see example map 5 below).
At present, no attempt has been made to deal with the overlap
at the International Date Line where Russia overlaps the Aleutian Islands of
Alaska. To some extent, the C-code can deal with this, but not in all cases.
Therefore, the user will likely experience some limitations when using
importGSHHS.
The solution is to import the whole dataset with this function using
xlim=c(0,360), and then apply the function refocusWorld
with user-desired X-limits.
The Y-limits are generally not problematic unless the user wants to focus on
either pole.
Nicholas M. Boers, Senior Software Engineer
Jobber, Edmonton AB
Maintainer: Rowan Haigh, Program Head – Offshore Rockfish
Pacific Biological Station (PBS), Fisheries & Oceans Canada (DFO), Nanaimo BC
locus opus: Offsite, Vancouver BC
Last modified Rd: 2024-09-25
Wessel, P., and Smith, W.H.F. (1996)
A Global Self-consistent, Hierarchical, High-resolution Shoreline Database.
J. Geophys. Res. 101 8741-8743.
(www.soest.hawaii.edu/pwessel/gshhg/Wessel+Smith_1996_JGR.pdf)
In package PBSmapping:importEvents,
importLocs,
importPolys
## Not run: useWest=FALSE useVers=c("2.2.0","2.2.3","2.3.0","2.3.4") # GSHHG versions mapswitch = 5 for (i in c("land","rivers","borders")) if (exists(i)) eval(parse(text=paste0("rm(",i,")"))) switch( mapswitch, # 1. Canada------------------------------------------------ {vN=4; useWest=T; xlim=c(-150,-50)+360;ylim=c(40,75)}, # 2. NW Canada & America----------------------------------- {vN=4; useWest=T;xlim=c(-136,-100)+360;ylim=c(40,75)}, # 3. Black Sea (user Ivailo)------------------------------- {vN=4; xlim=c(27.5, 34.3); ylim=c(40.9, 46.7)}, # 4. W Europe, NW Africa (user Uli)------------------------ {vN=4; xlim=c(-20,10); ylim=c(20,50)}, # 5. W Europe + Iceland------------------------------------ {vN=4; xlim=c(-25, 20); ylim=c(40, 68)}, # 6. New Zealand------------------------------------------- {vN=4; xlim=c(163, 182); ylim=c(-48,-34)}, # 7. Australia--------------------------------------------- {vN=4; xlim=c(112,155); ylim=c(-44,-10)}, # 8. Japan------------------------------------------------- {vN=4; xlim=c(127,148); ylim=c(30,47)}, # 9. Central America--------------------------------------- {vN=4; useWest=T; xlim=c(-95,-60)+360;ylim=c(-10,25)}, #10. North Pacific----------------------------------------- {vN=4; useWest=T; xlim=c(150,220); ylim=c(45,80)}, #11. Pacific Ocean----------------------------------------- {vN=4; xlim=c(112,240); ylim=c(-48,80)}, #12. North Atlantic (world coordinates)-------------------- {vN=4; xlim=c(285,360); ylim=c(40,68)}, #13. North Atlantic (western hemisphere coordinates)------- {vN=4; xlim=c(-75,0); ylim=c(40,68)}, #14. Atlantic Ocean---------------------------------------- {vN=4; xlim=c(285,380); ylim=c(-50,68)}, #15. Northern hemisphere----------------------------------- {vN=4; xlim=c(-180,180); ylim=c(0,85)}, #16. Asia-------------------------------------------------- {vN=4; xlim=c(0,180); ylim=c(0,80)}, #17. North America----------------------------------------- {vN=4; xlim=c(-180,0); ylim=c(0,80)}, #18. International date line------------------------------- {vN=4; xlim=c(45,315); ylim=c(0,80)}, #19. Indian Ocean------------------------------------------ {vN=4; xlim=c(20,130); ylim=c(-40,40)}, #20. Moose County ("400 miles north of everywhere")-------- {vN=4; xlim=c(272.5,280.5); ylim=c(43,47.5)} ) db=paste0("gshhg-bin-",useVers[vN]) # database version folder gshhg = paste0("C:/Ruser/GSHHG/",db,"/") # directory with binary files land = importGSHHS(paste0(gshhg,"gshhs_i.b"), xlim=xlim,ylim=ylim,maxLevel=4,useWest=useWest) rivers = importGSHHS(paste0(gshhg,"wdb_rivers_i.b"), xlim=xlim,ylim=ylim,useWest=useWest) borders = importGSHHS(paste0(gshhg,"wdb_borders_i.b"), xlim=xlim,ylim=ylim,useWest=useWest,maxLevel=1) if(exists("land")){ plotMap(land,xlim=xlim-ifelse(useWest,360,0),ylim=ylim, col="lemonchiffon",bg="aliceblue") if(!is.null(rivers)) addLines(rivers,col="blue") if(!is.null(borders)) addLines(borders,col="red",lwd=2) } ## End(Not run)## Not run: useWest=FALSE useVers=c("2.2.0","2.2.3","2.3.0","2.3.4") # GSHHG versions mapswitch = 5 for (i in c("land","rivers","borders")) if (exists(i)) eval(parse(text=paste0("rm(",i,")"))) switch( mapswitch, # 1. Canada------------------------------------------------ {vN=4; useWest=T; xlim=c(-150,-50)+360;ylim=c(40,75)}, # 2. NW Canada & America----------------------------------- {vN=4; useWest=T;xlim=c(-136,-100)+360;ylim=c(40,75)}, # 3. Black Sea (user Ivailo)------------------------------- {vN=4; xlim=c(27.5, 34.3); ylim=c(40.9, 46.7)}, # 4. W Europe, NW Africa (user Uli)------------------------ {vN=4; xlim=c(-20,10); ylim=c(20,50)}, # 5. W Europe + Iceland------------------------------------ {vN=4; xlim=c(-25, 20); ylim=c(40, 68)}, # 6. New Zealand------------------------------------------- {vN=4; xlim=c(163, 182); ylim=c(-48,-34)}, # 7. Australia--------------------------------------------- {vN=4; xlim=c(112,155); ylim=c(-44,-10)}, # 8. Japan------------------------------------------------- {vN=4; xlim=c(127,148); ylim=c(30,47)}, # 9. Central America--------------------------------------- {vN=4; useWest=T; xlim=c(-95,-60)+360;ylim=c(-10,25)}, #10. North Pacific----------------------------------------- {vN=4; useWest=T; xlim=c(150,220); ylim=c(45,80)}, #11. Pacific Ocean----------------------------------------- {vN=4; xlim=c(112,240); ylim=c(-48,80)}, #12. North Atlantic (world coordinates)-------------------- {vN=4; xlim=c(285,360); ylim=c(40,68)}, #13. North Atlantic (western hemisphere coordinates)------- {vN=4; xlim=c(-75,0); ylim=c(40,68)}, #14. Atlantic Ocean---------------------------------------- {vN=4; xlim=c(285,380); ylim=c(-50,68)}, #15. Northern hemisphere----------------------------------- {vN=4; xlim=c(-180,180); ylim=c(0,85)}, #16. Asia-------------------------------------------------- {vN=4; xlim=c(0,180); ylim=c(0,80)}, #17. North America----------------------------------------- {vN=4; xlim=c(-180,0); ylim=c(0,80)}, #18. International date line------------------------------- {vN=4; xlim=c(45,315); ylim=c(0,80)}, #19. Indian Ocean------------------------------------------ {vN=4; xlim=c(20,130); ylim=c(-40,40)}, #20. Moose County ("400 miles north of everywhere")-------- {vN=4; xlim=c(272.5,280.5); ylim=c(43,47.5)} ) db=paste0("gshhg-bin-",useVers[vN]) # database version folder gshhg = paste0("C:/Ruser/GSHHG/",db,"/") # directory with binary files land = importGSHHS(paste0(gshhg,"gshhs_i.b"), xlim=xlim,ylim=ylim,maxLevel=4,useWest=useWest) rivers = importGSHHS(paste0(gshhg,"wdb_rivers_i.b"), xlim=xlim,ylim=ylim,useWest=useWest) borders = importGSHHS(paste0(gshhg,"wdb_borders_i.b"), xlim=xlim,ylim=ylim,useWest=useWest,maxLevel=1) if(exists("land")){ plotMap(land,xlim=xlim-ifelse(useWest,360,0),ylim=ylim, col="lemonchiffon",bg="aliceblue") if(!is.null(rivers)) addLines(rivers,col="blue") if(!is.null(borders)) addLines(borders,col="red",lwd=2) } ## End(Not run)
Import a text file and convert it into a LocationSet.
importLocs(LocationSet)importLocs(LocationSet)
LocationSet |
|
An imported LocationSet.
Nicholas M. Boers, Staff Software Engineer
Jobber, Edmonton AB
Last modified Rd: 2023-10-30
importPolys, importEvents, importGSHHS
Import a text file and convert it into a PolySet with optional PolyData attribute.
importPolys(PolySet, PolyData=NULL, projection=NULL, zone=NULL)importPolys(PolySet, PolyData=NULL, projection=NULL, zone=NULL)
PolySet |
|
PolyData |
|
projection |
|
zone |
|
An imported PolySet with optional PolyData attribute.
Nicholas M. Boers, Staff Software Engineer
Jobber, Edmonton AB
Last modified Rd: 2023-10-30
importEvents, importLocs, importGSHHS
Determine whether polygons found in a PolySet are convex.
isConvex (polys)isConvex (polys)
polys |
PolySet to use. |
Convex polygons do not self-intersect. In a convex polygon, only the first and last vertices may share the same coordinates (i.e., the polygons are optionally closed).
The function does not give special consideration to holes. It returns
a value for each unique (PID, SID), regardless of
whether a contour represents a hole.
PolyData with columns PID, SID (may be missing),
and convex. Column convex contains Boolean values.
Nicholas M. Boers, Staff Software Engineer
Jobber, Edmonton AB
Last modified Rd: 2013-04-10
local(envir=.PBSmapEnv,expr={ #--- load the data (if using R) if (!is.null(version$language) && (version$language=="R")) data(nepacLL,envir=.PBSmapEnv) #--- calculate then print the polygons that are convex p <- isConvex(nepacLL); #--- nepacLL actually contains no convex polygons print(p[p$convex,]) })local(envir=.PBSmapEnv,expr={ #--- load the data (if using R) if (!is.null(version$language) && (version$language=="R")) data(nepacLL,envir=.PBSmapEnv) #--- calculate then print the polygons that are convex p <- isConvex(nepacLL); #--- nepacLL actually contains no convex polygons print(p[p$convex,]) })
Determine whether polygons found in a PolySet are self-intersecting.
isIntersecting (polys, numericResult = FALSE)isIntersecting (polys, numericResult = FALSE)
polys |
PolySet to use. |
numericResult |
Boolean value; if |
When numericResult = TRUE, this function counts intersections
as the algorithm processes them. It counts certain types (i.e., those
involving vertices and those where an edge retraces over an edge) more
than once.
The function does not give special consideration to holes. It returns
a value for each unique (PID, SID), regardless of
whether a contour represents a hole.
PolyData with columns PID, SID (may be missing),
and intersecting. If numericResult is TRUE,
intersecting contains the number of intersections. Otherwise,
it contains a Boolean value.
Nicholas M. Boers, Staff Software Engineer
Jobber, Edmonton AB
Last modified Rd: 2013-04-10
local(envir=.PBSmapEnv,expr={ #--- load the data (if using R) if (!is.null(version$language) && (version$language=="R")) data(nepacLL,envir=.PBSmapEnv) #--- calculate then print the polygons that are self-intersecting p <- isIntersecting(nepacLL, numericResult = FALSE) print(p[p$intersecting,]) })local(envir=.PBSmapEnv,expr={ #--- load the data (if using R) if (!is.null(version$language) && (version$language=="R")) data(nepacLL,envir=.PBSmapEnv) #--- calculate then print the polygons that are self-intersecting p <- isIntersecting(nepacLL, numericResult = FALSE) print(p[p$intersecting,]) })
Join one or two PolySets using a logic operation.
joinPolys(polysA, polysB=NULL, operation="INT")joinPolys(polysA, polysB=NULL, operation="INT")
polysA |
PolySet to join. |
polysB |
optional second PolySet with which to join. |
operation |
one of |
This function interfaces with the first version of the Clipper Library,
specifically version 6.2.1 released 2014-10-31, developed by Angus Johnson.
Angus now offers Clipper 2 Library.
Prior to 2013-03-23, 'joinPolys' used the General Polygon Clipper library
by Alan Murta at the University of Manchester.
We keep this historic reference to GPC because 'joinPolys' remains faithful
to Murta's definition of a generic polygon, which we describe below.
Murta (2004) defines a generic polygon (or polygon set)
as zero or more disjoint boundaries of arbitrary configuration. He
relates a boundary to a contour, where each may be convex,
concave or self-intersecting. In a PolySet, the polygons associated
with each unique PID loosely correspond to a generic polygon,
as they can represent both inner and outer boundaries. Our use of the
term generic polygon includes the restrictions imposed by a
PolySet. For example, the polygons for a given PID cannot
be arranged arbitrarily.
If 'polysB' is NULL, this function sequentially applies
the logic 'operation' between the generic polygons in 'polysA'.
For example, suppose 'polysA' contains three generic polygons (A, B, C).
The function outputs the PolySet containing ((A op B) op C).
If 'polysB' is not NULL, this function applies the logic 'operation'
between each generic polygon in 'polysA' and each one in 'polysB'.
For example, suppose 'polysA' contains two generic polygons (A, B)
and 'polysB' contains two generic polygons (C, D).
The function's output is the concatenation of
A op C, B op C, A op D, B op D, with PIDs 1 to 4, respectively.
Generally there are n times m comparisons, where n = number of polygons
in 'polysA' and m = number of polygons in 'polysB'.
If 'polysB' contains only one generic polygon, the function
maintains the PIDs from 'polysA'.
It also maintains them when 'polysA' contains only one generic polygon and
the 'operation' is "DIFF" (difference).
Otherwise, if 'polysA' contains only one generic polygon, it maintains
the PIDs from 'polysB'.
If 'polysB' is NULL, the resulting PolySet contains
a single generic polygon (one PID), possibly with several
components (SIDs).
The function recalculates the PID and SID columns.
If 'polysB' is not NULL, the resulting PolySet contains one or more
generic polygons (PIDs), each with possibly several components (SIDs).
The function recalculates the SID column, and depending on the input,
it may recalculate the PID column.
C code: Angus Johnson, Computer Programmer
Implementation: Nicholas M. Boers, Staff Software Engineer
Jobber, Edmonton AB
Maintainer: Rowan Haigh, Program Head – Offshore Rockfish
Pacific Biological Station (PBS), Fisheries & Oceans Canada (DFO), Nanaimo BC
locus opus: Remote office, Vancouver BC
Last modified Rd: 2024-09-03
Murta, A. (2004) A General Polygon Clipping Library. Accessed: Jul 29, 2004.
Johnson, A. (2014) Clipper – an open source freeware library for clipping and offsetting lines and polygons. Accessed: Oct 31, 2014.
In package PBSmapping:addPolys,
appendPolys,
clipPolys,
closePolys,
fixBound,
fixPOS,
locatePolys,
plotMap,
plotPoints,
thickenPolys,
thinPolys
local(envir=.PBSmapEnv,expr={ oldpar = par(no.readonly=TRUE) #--- load the data (if using R) if (!is.null(version$language) && (version$language=="R")) data(nepacLL,envir=.PBSmapEnv) ### Example 1. Cut a triangle out of Vancouver Island par(mfrow=c(1,1)) #--- create a triangle to use in clipping polysB <- data.frame(PID=rep(1, 3), POS=1:3, X=c(-127.5, -124.5, -125.6), Y = c(49.2, 50.3, 48.6)) #--- intersect nepacLL with the single polygon, and plot the result plotMap(joinPolys(nepacLL, polysB), col="cyan") #--- add nepacLL in a different line type to emphasize the intersection addPolys(nepacLL, border="purple", lty=3, density=0) box() ### Example 2. Cut Texada and Lasqueti Islands out of Boxes xlim = list(box1=c(-124.8,-124),box2=c(-124,-123.9)) ylim = list(box1=c(49.4,49.85), box2=c(49.85,49.9)) Xlim = extendrange(xlim); Ylim=extendrange(ylim) polyA = as.PolySet(data.frame( PID = rep(1:2,each=4), POS = rep(1:4,2), X = as.vector(sapply(xlim,function(x){x[c(1,1,2,2)]})), Y = as.vector(sapply(ylim,function(x){x[c(1,2,2,1)]})) ), projection="LL") data(nepacLLhigh,envir=.PBSmapEnv) polyB = nepacLLhigh[is.element(nepacLLhigh$PID,c(736,1912)),] polyC = joinPolys(polyA, polyB, "DIFF") par(mfrow=c(2,2),cex=1,mgp=c(2,0.5,0)) plotMap(polyA,col="lightblue",xlim=Xlim,ylim=Ylim) addPolys(polyB,col="gold"); text(mean(Xlim)-0.05,Ylim-0.04,"Boxes (A,B) and Isles (C,D)") labs = calcCentroid(polyA) labs[1,c("X","Y")] = labs[2,c("X","Y")]+c(-0.1,-0.05) text(labs[,"X"],labs[,"Y"],c("A","B"),font=2) plotMap(polyC[is.element(polyC$PID,1),],col="pink",xlim=Xlim,ylim=Ylim) text(mean(Xlim)-0.05,Ylim-0.04,"Box A \"DIFF\" Isle C") plotMap(polyC[is.element(polyC$PID,3),],col="green",xlim=Xlim,ylim=Ylim) text(mean(Xlim)-0.05,Ylim-0.04,"Box A \"DIFF\" Isle D") plotMap(polyC[is.element(polyC$PID,c(1,3)),],col="cyan",xlim=Xlim,ylim=Ylim) text(mean(Xlim)-0.05,Ylim-0.04,"Box A \"DIFF\" Isles (C,D)") par(oldpar) })local(envir=.PBSmapEnv,expr={ oldpar = par(no.readonly=TRUE) #--- load the data (if using R) if (!is.null(version$language) && (version$language=="R")) data(nepacLL,envir=.PBSmapEnv) ### Example 1. Cut a triangle out of Vancouver Island par(mfrow=c(1,1)) #--- create a triangle to use in clipping polysB <- data.frame(PID=rep(1, 3), POS=1:3, X=c(-127.5, -124.5, -125.6), Y = c(49.2, 50.3, 48.6)) #--- intersect nepacLL with the single polygon, and plot the result plotMap(joinPolys(nepacLL, polysB), col="cyan") #--- add nepacLL in a different line type to emphasize the intersection addPolys(nepacLL, border="purple", lty=3, density=0) box() ### Example 2. Cut Texada and Lasqueti Islands out of Boxes xlim = list(box1=c(-124.8,-124),box2=c(-124,-123.9)) ylim = list(box1=c(49.4,49.85), box2=c(49.85,49.9)) Xlim = extendrange(xlim); Ylim=extendrange(ylim) polyA = as.PolySet(data.frame( PID = rep(1:2,each=4), POS = rep(1:4,2), X = as.vector(sapply(xlim,function(x){x[c(1,1,2,2)]})), Y = as.vector(sapply(ylim,function(x){x[c(1,2,2,1)]})) ), projection="LL") data(nepacLLhigh,envir=.PBSmapEnv) polyB = nepacLLhigh[is.element(nepacLLhigh$PID,c(736,1912)),] polyC = joinPolys(polyA, polyB, "DIFF") par(mfrow=c(2,2),cex=1,mgp=c(2,0.5,0)) plotMap(polyA,col="lightblue",xlim=Xlim,ylim=Ylim) addPolys(polyB,col="gold"); text(mean(Xlim)-0.05,Ylim-0.04,"Boxes (A,B) and Isles (C,D)") labs = calcCentroid(polyA) labs[1,c("X","Y")] = labs[2,c("X","Y")]+c(-0.1,-0.05) text(labs[,"X"],labs[,"Y"],c("A","B"),font=2) plotMap(polyC[is.element(polyC$PID,1),],col="pink",xlim=Xlim,ylim=Ylim) text(mean(Xlim)-0.05,Ylim-0.04,"Box A \"DIFF\" Isle C") plotMap(polyC[is.element(polyC$PID,3),],col="green",xlim=Xlim,ylim=Ylim) text(mean(Xlim)-0.05,Ylim-0.04,"Box A \"DIFF\" Isle D") plotMap(polyC[is.element(polyC$PID,c(1,3)),],col="cyan",xlim=Xlim,ylim=Ylim) text(mean(Xlim)-0.05,Ylim-0.04,"Box A \"DIFF\" Isles (C,D)") par(oldpar) })
Locate events on the current plot (using the locator
function).
locateEvents (EID, n = 512, type = "p", ...)locateEvents (EID, n = 512, type = "p", ...)
EID |
vector of event IDs (optional). |
n |
maximum number of events to locate. |
type |
one of |
... |
This function allows its user to define events with mouse clicks on
the current plot via the locator function. The
arguments n and type are the usual parameters of the
locator function. If EID is not missing, then
n = length(EID).
On exit from locator, suppose the user defined m
events. If EID was missing, then the output data frame will
contain m events. However, if EID exists, then the
output data frame will contain length(EID) events, and both
X and Y will be NA for events
EID[(m+1):n]. The na.omit function
can remove rows with NAs.
EventData with columns EID, X, and Y, and
projection attribute equal to the map's projection. The
function does not set the zone attribute.
Nicholas M. Boers, Staff Software Engineer
Jobber, Edmonton AB
Last modified Rd: 2007-06-06
addPoints,
combineEvents,
convDP,
EventData,
findCells,
findPolys,
plotPoints.
#--- define five events on the current plot, numbering them 10 to 14 ## Not run: events <- locateEvents(EID = 10:14)#--- define five events on the current plot, numbering them 10 to 14 ## Not run: events <- locateEvents(EID = 10:14)
Locate polygons on the current plot (using the locator
function).
locatePolys (pdata, n = 512, type = "o", ...)locatePolys (pdata, n = 512, type = "o", ...)
pdata |
PolyData (optional) with
columns |
n |
maximum number of points to locate. |
type |
one of |
... |
This function allows its user to define polygons with mouse clicks on
the current plot via the locator function. The
arguments n and type are the usual parameters for the
locator function, but the user can specify them for each
individual (PID, SID) in a pdata object.
If a pdata object exists, the function ignores columns other
than PID, SID, n, and type. If pdata
includes n, then an outer boundary has n > 0 and an
inner boundary has n < 0.
On exit from locator, suppose the user defined m
vertices for a given polygon. For that polygon, the X and
Y columns will contain NAs where POS =
(m+1):n for outer-boundaries and POS =
(|n|-m):1 for inner-boundaries. The
na.omit function can remove rows with NAs.
If a pdata object does not exist, the output contains only one
polygon with a PID equal to 1. One inner-boundary polygon
(POS goes from n to 1) can be generated by
supplying a negative n.
If type = "o" or type = "l", the function draws a line
connecting the last and first vertices.
PolySet with projection attribute equal to the map's
projection. The function does not set the zone attribute.
Nicholas M. Boers, Staff Software Engineer
Jobber, Edmonton AB
Last modified Rd: 2007-06-06
addPolys,
appendPolys,
clipPolys,
closePolys,
findCells,
findPolys,
fixPOS,
joinPolys,
plotMap,
plotPolys,
thickenPolys,
thinPolys.
#--- define one polygon with up to 5 vertices on the current plot ## Not run: polys <- locatePolys(n = 5)#--- define one polygon with up to 5 vertices on the current plot ## Not run: polys <- locatePolys(n = 5)
A LocationSet comprises a data frame that summarises which EventData points
(EID) lie in which PolySet polygons (PID) or
(PID, SID).
Events not located in target polygons are not reported. If an event lies on a
polygon boundary, an additional LocationSet field called Bdry is set
to TRUE. One event can also occur in multiple polygons.
PBSmapping functions that expect LocationSet's will accept properly formatted data frames in their place (see 'Details').
as.LocationSet attempts to coerce a data frame to an object with
class LocationSet.
is.LocationSet returns TRUE if its argument is of class
LocationSet.
as.LocationSet(x) is.LocationSet(x, fullValidation = TRUE)as.LocationSet(x) is.LocationSet(x, fullValidation = TRUE)
x |
data frame to be coerced or tested. |
fullValidation |
Boolean value; if |
A PolySet can define regional boundaries for drawing a map, and
EventData can give event points on the map. Which events occur in
which regions? Our function findPolys resolves this
problem. The output lies in a LocationSet, a data frame with three or
four columns (EID, PID, SID, Bdry), where
SID may be missing. One row in a LocationSet means that the event
EID occurs in the polygon (PID, SID). The boundary
(Bdry) field specifies whether (Bdry=T) or not
(Bdry=F) the event lies on the polygon boundary. If SID
refers to an inner polygon boundary, then EID occurs in
(PID, SID) only if Bdry=T. An event may occur in
multiple polygons. Thus, the same EID can occur in multiple
records. If an EID does not fall in any (PID, SID),
or if it falls within a hole, it does not occur in the output
LocationSet. Inserting the string "LocationSet" as the first
element of a LocationSet's class attribute alters the behaviour
of some functions, including print (if
PBSprint is TRUE) and summary.
The as.LocationSet method returns an object with classes
"LocationSet" and "data.frame", in that order.
Nicholas M. Boers, Staff Software Engineer
Jobber, Edmonton AB
Last modified Rd: 2015-04-23
Make a grid of polygons, using PIDs and SIDs according
to the input arguments.
makeGrid(x, y, byrow=TRUE, addSID=TRUE, projection=NULL, zone = NULL, type="rectangle")makeGrid(x, y, byrow=TRUE, addSID=TRUE, projection=NULL, zone = NULL, type="rectangle")
x |
|
y |
|
byrow |
|
addSID |
|
projection |
|
zone |
|
type |
|
This function makes a grid of polygons, labeling them according to
byrow and addSID.
For rectangular tesselations (grid cells), the variables and
indicate column and row numbers, respectively, where the lower-left cell
of the grid is (1, 1):
byrow TRUE and addSID
FALSE implies PID
byrow FALSE and addSID
FALSE implies PID
byrow TRUE and addSID
TRUE implies PID , SID
byrow FALSE and addSID
TRUE implies PID , SID
For hexagonal tesselations (grid cells), indicates columns for flat-topped
hexagons and rows for pointy-topped hexagons. The reverse is true for .
Stemming from their six-sided nature, hexagons will adjoin along a long-edge by row when
their orientation is such that one vertex is higher than all the others.
Hexagons will adjoin along a long-edge by column when their orientation shows two
uppermost vertices.
PolySet with columns PID, SID
(if addSID=TRUE), POS, X, and Y.
The PolySet is a set of rectangular grid cells when type='rectangle', with
vertices:.
The PolySet is a set of hexagonal grid cells when type='hexagon'.
Nicholas M. Boers, Staff Software Engineer
Jobber, Edmonton AB
Rowan Haigh, Program Head – Offshore Rockfish
Pacific Biological Station (PBS), Fisheries & Oceans Canada (DFO), Nanaimo BC
locus opus: Institute of Ocean Sciences (IOS), Sidney BC
Last modified Rd: 2019-01-04
addPolys,
clipPolys,
combineEvents,
findCells,
findPolys,
thickenPolys.
local(envir=.PBSmapEnv,expr={ oldpar = par(no.readonly=TRUE) ##--- make a 10 x 10 grid polyGrid <- makeGrid(x=0:10, y=0:10) ##--- plot the grid plotPolys(polyGrid, density=0, projection=1) par(oldpar) })local(envir=.PBSmapEnv,expr={ oldpar = par(no.readonly=TRUE) ##--- make a 10 x 10 grid polyGrid <- makeGrid(x=0:10, y=0:10) ##--- plot the grid plotPolys(polyGrid, density=0, projection=1) par(oldpar) })
Append a column for a polygon property (e.g., border or
lty) to PolyData based on measurements in the
PolyData's Z column.
makeProps(pdata,breaks,propName="col",propVals=1:(length(breaks)-1))makeProps(pdata,breaks,propName="col",propVals=1:(length(breaks)-1))
pdata |
PolyData with a |
breaks |
either a vector of cut points or a scalar denoting the
number of intervals that |
propName |
name of the new column to append to |
propVals |
vector of values to associate with |
This function acts like the cut function to produce
PolyData suitable for the polyProps plotting argument
(see addLabels, addLines,
addPoints, addPolys,
addStipples, plotLines,
plotMap,plotPoints, and
plotPolys). The Z column of pdata is
equivalent to the data vector x of the cut
function.
PolyData with the same columns as pdata plus an
additional column propName.
Nicholas M. Boers, Staff Software Engineer
Jobber, Edmonton AB
Last modified Rd: 2013-04-10
addLabels,
addLines,
addPoints,
addPolys,
addStipples,
plotLines,
plotMap,
plotPoints,
plotPolys,
PolyData,
PolySet.
local(envir=.PBSmapEnv,expr={ #--- create a PolyData object pd <- data.frame(PID=1:10, Z=1:10) #--- using 3 intervals, create a column named `col' and populate it with #--- the supplied values makeProps(pdata=pd, breaks=3, propName="col", propVals=c(1:3)) })local(envir=.PBSmapEnv,expr={ #--- create a PolyData object pd <- data.frame(PID=1:10, Z=1:10) #--- using 3 intervals, create a column named `col' and populate it with #--- the supplied values makeProps(pdata=pd, breaks=3, propName="col", propVals=c(1:3)) })
Make topography data suitable for the 'graphics::contour' and
'grDevices::contourLines' functions using freely available global
seafloor topography data.
makeTopography (dat, digits=2, func=NULL)makeTopography (dat, digits=2, func=NULL)
dat |
|
digits |
|
func |
|
Suitable data can be obtained through the Topex acquisition form.
The function 'utils::read.table' will import dowloaded ASCII files into R,
creating objects suitable for the argument 'dat' in 'makeTopography'.
When creating data for regions with longitude values spanning
-180 to 0, consider
subtracting 360 from the result's longitude coordinates (x).
When creating bathymetry data, consider negating the result's
elevations (z) to give depths positive values.
Combinations of (x,y) do not need to be complete (z[x,y]=NA) or
unique (z[x,y] = func(z[x,y])).
List with elements x, y, and z. Elements x and
y are vectors, while z is a matrix with rownames x
and colnames y.
The functions 'graphics::contour' and 'grDevices::contourLines'
expect data conforming to this list format.
Nicholas M. Boers, Staff Software Engineer
Jobber, Edmonton AB
Maintainer: Rowan Haigh, Program Head – Offshore Rockfish
Pacific Biological Station (PBS), Fisheries & Oceans Canada (DFO), Nanaimo BC
locus opus: Institute of Ocean Sciences (IOS), Sidney BC
Last modified Rd: 2021-01-11
In package graphics:contour
In package grDevices:contourLines
In package PBSmapping:convCP
local(envir=.PBSmapEnv,expr={ oldpar = par(no.readonly=TRUE) #--- Example 1: Sample data frame and conversion. file <- data.frame(X=c(1,1,2,2),Y=c(3,4,3,4),Z=c(5,6,7,8)) print(makeTopography(file)) #--- Example 2: Aleutian Islands bathymetry isob <- c(100,500,1000,2500,5000) icol <- rgb(0,0,seq(255,100,len=length(isob)),max=255) afile <- paste(system.file(package="PBSmapping"), "/Extra/aleutian.txt",sep="") aleutian <- read.table(afile, header=FALSE, col.names=c("x","y","z")) aleutian$x <- aleutian$x - 360 aleutian$z <- -aleutian$z alBathy <- makeTopography(aleutian) alCL <- contourLines(alBathy,levels=isob) alCP <- convCP(alCL) alPoly <- alCP$PolySet attr(alPoly,"projection") <- "LL" plotMap(alPoly, type="n", cex.axis=1.2, cex.lab=1.5) addLines(alPoly,col=icol) data(nepacLL,envir=.PBSmapEnv) addPolys(nepacLL,col="gold") legend(x="topleft",bty="n",col=icol,lwd=2,legend=as.character(isob)) par(oldpar) })local(envir=.PBSmapEnv,expr={ oldpar = par(no.readonly=TRUE) #--- Example 1: Sample data frame and conversion. file <- data.frame(X=c(1,1,2,2),Y=c(3,4,3,4),Z=c(5,6,7,8)) print(makeTopography(file)) #--- Example 2: Aleutian Islands bathymetry isob <- c(100,500,1000,2500,5000) icol <- rgb(0,0,seq(255,100,len=length(isob)),max=255) afile <- paste(system.file(package="PBSmapping"), "/Extra/aleutian.txt",sep="") aleutian <- read.table(afile, header=FALSE, col.names=c("x","y","z")) aleutian$x <- aleutian$x - 360 aleutian$z <- -aleutian$z alBathy <- makeTopography(aleutian) alCL <- contourLines(alBathy,levels=isob) alCP <- convCP(alCL) alPoly <- alCP$PolySet attr(alPoly,"projection") <- "LL" plotMap(alPoly, type="n", cex.axis=1.2, cex.lab=1.5) addLines(alPoly,col=icol) data(nepacLL,envir=.PBSmapEnv) addPolys(nepacLL,col="gold") legend(x="topleft",bty="n",col=icol,lwd=2,legend=as.character(isob)) par(oldpar) })
PolySet of polygons for the shorelines of the northeast Pacific Ocean and of the world, both in normal and high resolution.
data(nepacLL) data(nepacLLhigh) data(worldLL) data(worldLLhigh)data(nepacLL) data(nepacLLhigh) data(worldLL) data(worldLLhigh)
Data frame consisting of 4 columns: PID = primary polygon ID,
POS = position of each vertex within a given polygon, X
= longitude coordinate, and Y = latitude coordinate. Attributes:
projection = "LL".
In R, the data must be loaded using the data function.
Nicholas M. Boers, Staff Software Engineer
Jobber, Edmonton AB
Last modified Rd: 2023-10-30
Polygon data from the GSHHG (Global Self-consistent, Hierarchical,
High-resolution Geography) Database.
Download the native binary files of shoreline polygons, rivers, and borders
contained in the latest zip archive (version 2.3.4 ) at
(www.soest.hawaii.edu/pwessel/gshhg).
nepacLL <- importGSHHS("gshhs_h.b", xlim=c(-190,-110), ylim=c(34,72),
level=1, n=15, xoff=-360)
nepacLLhigh <- importGSHHS("gshhs_f.b", xlim=c(-190,-110),
ylim=c(34,72), level=1, n=0, xoff=-360)
nepacLLhigh <- thinPolys(nepacLLhigh, tol=0.1, filter=3)
worldLL <- importGSHHS("gshhs_l.b", xlim=c(-20,360), ylim=c(-90,90),
level=1, n=15, xoff=0)
worldLL <- PBSmapping:::.fixGSHHSWorld(worldLL)
worldLLhigh <- importGSHHS("gshhs_i.b", xlim=c(-20,360),
ylim=c(-90,90), level=1, n=15, xoff=0)
worldLLhigh <- PBSmapping:::.fixGSHHSWorld(worldLLhigh)
Wessel, P. and Smith, W.H.F. (1996) A global, self-consistent,
hierarchical, high-resolution shoreline database. Journal of
Geophysical Research 101, 8741–8743.
(www.soest.hawaii.edu/pwessel/gshhg/Wessel+Smith_1996_JGR.pdf)
Data:bcBathymetry,
surveyData,
towData
Functions:importGSHHS,
plotMap,
plotPolys,
addPolys,
clipPolys,
refocusWorld,
thickenPolys,
thinPolys
This software has evolved from fisheries research conducted at the Pacific Biological Station (PBS) in Nanaimo, British Columbia, Canada. It extends the R language to include two-dimensional plotting features similar to those commonly available in a Geographic Information System (GIS). Embedded C code speeds algorithms from computational geometry, such as finding polygons that contain specified point events or converting between longitude-latitude and Universal Transverse Mercator (UTM) coordinates. It includes data for a global shoreline and other data sets in the public domain.
For a complete user's guide, see the file PBSmapping-UG.pdf
in the R directory .../library/PBSmapping/doc.
PBSmapping includes 10 demos that appear as figures in
the User's Guide. To see them, run the function .PBSfigs().
More generally, a user can view all demos available from locally
installed packages with the function runDemos() in our
related (and recommended) package PBSmodelling.
Nicholas M. Boers, Staff Software Engineer
Jobber, Edmonton AB
Last modified Rd: 2008-09-03
Specify whether PBS Mapping should print object summaries or not. If not, data objects are displayed as normal.
PBSprintPBSprint
If PBSprint = TRUE, the mapping software will print summaries
rather than the data frames for EventData, LocationSet, PolyData, and
PolySet objects. If PBSprint = FALSE, it will print the data
frames.
This variable's default value is FALSE.
TRUE or FALSE, depending on the user's preference.
Nicholas M. Boers, Staff Software Engineer
Jobber, Edmonton AB
Last modified Rd: 2007-06-06
Place secondary polygons identified as holes (inner contours) under primary polygons identified as solids (outer contours) if the vertices of the holes lie completely within the vertices of the solids. This operation is performed for each primary polygon until all holes have been assigned.
placeHoles(polyset, minVerts=3, orient=FALSE, show.progress=FALSE, unique.pid=FALSE)placeHoles(polyset, minVerts=3, orient=FALSE, show.progress=FALSE, unique.pid=FALSE)
polyset |
a valid PBSmapping PolySet. |
minVerts |
|
orient |
|
show.progress |
|
unique.pid |
|
The algorithm identifies outer contours (solids) and inner contours (holes)
using either the default PBSmapping method (solids = increasing "POS",
holes = decreasing "POS") or the rotational direction of the polygons
(solids = clockwise, holes = counter-clockwise).
It then systematically starts matching holes with solids based on their vertices being completely within the boundaries of the solid. If a hole happens to match a current solid completley (all vertices on the boundary), then the hole is not matched to this solid because it is a hole in another solid that creates space for the current solid.
To facilitate computation, the algorithm assumes that once a hole is located in a solid, it will not occur in any other solid. This means that for each successive solid, the number of candidate holes will either decrease or stay the same.
This function makes use of the PBSmapping hidden function ".is.in"
which uses the C code "findPolys". The latter only returns events found
in a polygon (or on the boundary) but .is.in evaluates all events and
returns a list containing:"e.in" – events within the polygon,"e.out" – events outside the polygon,"all.in" – logical value of whether all events are in the polygon,"all.out" – logical value of whether all events are outside the polygon,"all.bdry" – logical value of whether all events occur on the boundary of the polygon
Returns the input PolySet where holes have been arranged
beneath appropriate solids for each PID (original or redefined).
Rowan Haigh, Program Head – Offshore Rockfish
Pacific Biological Station (PBS), Fisheries & Oceans Canada (DFO), Nanaimo BC
locus opus: Offsite, Vancouver BC
Last modified Rd: 2023-10-30
In package PBSmapping:findPolys,
is.PolySet,
dot-is.in
Plot a PolySet as polylines.
plotLines (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, ...)plotLines (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, ...)
polys |
PolySet to plot (required). |
xlim |
range of X-coordinates. |
ylim |
range of Y-coordinates. |
projection |
desired projection when PolySet lacks a
|
plt |
four element numeric vector |
polyProps |
PolyData specifying which polylines to plot and their
properties. |
lty |
vector describing line types (cycled by |
col |
vector describing colours (cycled by |
bg |
background colour of the plot. |
axes |
Boolean value; if |
tckLab |
Boolean vector (length 1 or 2); if |
tck |
numeric vector (length 1 or 2) describing the length
of tick marks as a fraction of the smallest dimension. If
|
tckMinor |
numeric vector (length 1 or 2) describing the length of tick marks as a fraction of the smallest dimension. These tick marks can not be automatically labelled. If given a two-element vector, the first element describes the tick marks on the x-axis and the second element describes those on the y-axis. |
... |
additional |
This function plots a PolySet, where each unique (PID,
SID) describes a polyline. It does not connect each polyline's
last vertex to its first. Unlike plotMap, the function
ignores the aspect ratio. It clips polys to xlim and
ylim before plotting.
The function creates a blank plot when polys equals
NULL. In this case, the user must supply both xlim and
ylim arguments. Alternatively, it accepts the argument
type = "n" as part of ..., which is equivalent to specifying
polys = NULL, but requires a PolySet. In both cases,
the function's behaviour changes slightly. To resemble the
plot function, it plots the border, labels, and other
parts according to par parameters such as col.
For additional help on the arguments lty and col, please
see par.
PolyData consisting of the PolyProps used to create the plot.
To satisfy the aspect ratio, this plotting routine resizes the plot
region. Consequently, par parameters such as
plt, mai, and mar will change. When the function
terminates, these changes persist to allow for additions to the plot.
Nicholas M. Boers, Staff Software Engineer
Jobber, Edmonton AB
Last modified Rd: 2013-04-10
addLines, calcLength, clipLines,
closePolys, convLP, fixBound,
fixPOS,
locatePolys, thinPolys, thickenPolys.
local(envir=.PBSmapEnv,expr={ oldpar = par(no.readonly=TRUE) #--- create a PolySet to plot polys <- data.frame(PID=rep(1,4),POS=1:4,X=c(0,1,1,0),Y=c(0,0,1,1)) #--- plot the PolySet plotLines(polys, xlim=c(-.5,1.5), ylim=c(-.5,1.5)) par(oldpar) })local(envir=.PBSmapEnv,expr={ oldpar = par(no.readonly=TRUE) #--- create a PolySet to plot polys <- data.frame(PID=rep(1,4),POS=1:4,X=c(0,1,1,0),Y=c(0,0,1,1)) #--- plot the PolySet plotLines(polys, xlim=c(-.5,1.5), ylim=c(-.5,1.5)) par(oldpar) })
Plot a PolySet as a map, using the correct aspect ratio.
plotMap(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, ...) .plotMaps(polys, xlim, ylim, projection, plt, polyProps, border, lty, col, colHoles, density, angle, bg, axes, tckLab, tck, tckMinor, isType, ...) .initPlotRegion(projection, xlim, ylim, plt)plotMap(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, ...) .plotMaps(polys, xlim, ylim, projection, plt, polyProps, border, lty, col, colHoles, density, angle, bg, axes, tckLab, tck, tckMinor, isType, ...) .initPlotRegion(projection, xlim, ylim, plt)
polys |
|
xlim |
|
ylim |
|
projection |
|
plt |
|
polyProps |
|
border |
|
lty |
|
col |
|
colHoles |
|
density |
|
angle |
|
bg |
|
axes |
|
tckLab |
|
tck |
|
tckMinor |
|
... |
|
isType |
|
This function plots a PolySet, where each unique (PID,
SID) describes a polygon. It connects each polygon's last
vertex to its first. The function supports both borders
(border, lty) and fills (col, density,
angle). When supplied with the appropriate arguments, it can
draw only borders or only fills . Unlike plotLines and
plotPolys, it uses the aspect ratio supplied in the
projection attribute of polys. If this attribute is
missing, it attempts to use its projection argument. In the
absence of both, it uses a default aspect ratio of 1:1. It clips
polys to xlim and ylim before plotting.
The function creates a blank plot when polys equals
NULL. In this case, the user must supply both xlim and
ylim arguments. Alternatively, it accepts the argument
type = "n" as part of ..., which is equivalent to specifying
polys = NULL, but requires a PolySet. In both cases,
the function's behaviour changes slightly. To resemble the
plot function, it plots the border, labels, and other
parts according to par parameters such as col.
For additional help on the arguments border, lty,
col, density, and angle, please see
polygon and par.
PolyData consisting of the PolyProps used to create the plot.
To satisfy the aspect ratio, this plotting routine resizes the plot
region. Consequently, par parameters such as
plt, mai, and mar will change. When the function
terminates, these changes persist to allow for additions to the plot.
Auxiliary dot function '.initPlotRegion' initialises the plot region,
accounting for the aspect ratio.
Nicholas M. Boers, Software Engineer, Jobber, Edmonton AB
Maintainer: Rowan Haigh, Program Head – Offshore Rockfish
Pacific Biological Station (PBS), Fisheries & Oceans Canada (DFO), Nanaimo BC
locus opus: Remote office, Vancouver BC
Last modified Rd: 2024-09-03
In package PBSmapping:addLabels,
addPolys,
addStipples,
clipPolys,
closePolys,
fixBound,
fixPOS,
locatePolys,
plotLines,
plotPoints,
thinPolys,
thickenPolys.
local(envir=.PBSmapEnv,expr={ oldpar = par(no.readonly=TRUE) #--- create a PolySet to plot polys <- data.frame(PID=rep(1,4),POS=1:4,X=c(0,1,1,0),Y=c(0,0,1,1)) #--- plot the PolySet plotMap(polys,xlim=c(-.5,1.5),ylim=c(-.5,1.5),density=0,projection=1) par(oldpar) })local(envir=.PBSmapEnv,expr={ oldpar = par(no.readonly=TRUE) #--- create a PolySet to plot polys <- data.frame(PID=rep(1,4),POS=1:4,X=c(0,1,1,0),Y=c(0,0,1,1)) #--- plot the PolySet plotMap(polys,xlim=c(-.5,1.5),ylim=c(-.5,1.5),density=0,projection=1) par(oldpar) })
Plot EventData/PolyData, where each unique EID or
(PID, SID) describes a point.
plotPoints (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, ...)plotPoints (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, ...)
data |
|
xlim |
range of X-coordinates. |
ylim |
range of Y-coordinates. |
projection |
desired projection when PolySet lacks a
|
plt |
four element numeric vector |
polyProps |
PolyData specifying which points to plot and their
properties. |
cex |
vector describing character expansion factors (cycled by
|
col |
vector describing colours (cycled by |
pch |
vector describing plotting characters (cycled by |
axes |
Boolean value; if |
tckLab |
Boolean vector (length 1 or 2); if |
tck |
numeric vector (length 1 or 2) describing the length
of tick marks as a fraction of the smallest dimension. If
|
tckMinor |
numeric vector (length 1 or 2) describing the length of tick marks as a fraction of the smallest dimension. These tick marks can not be automatically labelled. If given a two-element vector, the first element describes the tick marks on the x-axis and the second element describes those on the y-axis. |
... |
additional |
This function clips data to xlim and ylim before
plotting. It only adds PolyData containing X and
Y columns.
The function creates a blank plot when polys equals
NULL. In this case, the user must supply both xlim and
ylim arguments. Alternatively, it accepts the argument
type = "n" as part of ..., which is equivalent to specifying
polys = NULL, but requires a PolySet. In both cases,
the function's behaviour changes slightly. To resemble the
plot function, it plots the border, labels, and other
parts according to par parameters such as col.
For additional help on the arguments cex, col, and
pch, please see par.
PolyData consisting of the PolyProps used to create the plot.
To satisfy the aspect ratio, this plotting routine resizes the plot
region. Consequently, par parameters such as
plt, mai, and mar will change. When the function
terminates, these changes persist to allow for additions to the plot.
Nicholas M. Boers, Staff Software Engineer
Jobber, Edmonton AB
Last modified Rd: 2013-04-10
addPoints,
combineEvents,
convDP,
findPolys,
locateEvents.
local(envir=.PBSmapEnv,expr={ oldpar = par(no.readonly=TRUE) #--- load the data (if using R) if (!is.null(version$language) && (version$language=="R")) data(nepacLL,surveyData,envir=.PBSmapEnv) #--- plot a map plotMap(nepacLL, xlim=c(-136, -125), ylim=c(48, 57)) #--- add events addPoints(surveyData, col=1:7) par(oldpar) })local(envir=.PBSmapEnv,expr={ oldpar = par(no.readonly=TRUE) #--- load the data (if using R) if (!is.null(version$language) && (version$language=="R")) data(nepacLL,surveyData,envir=.PBSmapEnv) #--- plot a map plotMap(nepacLL, xlim=c(-136, -125), ylim=c(48, 57)) #--- add events addPoints(surveyData, col=1:7) par(oldpar) })
Plot a PolySet as polygons.
plotPolys (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, ...)plotPolys (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, ...)
polys |
PolySet to plot (required). |
xlim |
range of X-coordinates. |
ylim |
range of Y-coordinates. |
projection |
desired projection when PolySet lacks a
|
plt |
four element numeric vector |
polyProps |
PolyData specifying which polygons to
plot and their properties. |
border |
vector describing edge colours (cycled by |
lty |
vector describing line types (cycled by |
col |
vector describing fill colours (cycled by |
colHoles |
vector describing hole colours (cycled by |
density |
vector describing shading line densities (lines per
inch, cycled by |
angle |
vector describing shading line angles (degrees, cycled by
|
bg |
background colour of the plot. |
axes |
Boolean value; if |
tckLab |
Boolean vector (length 1 or 2); if |
tck |
numeric vector (length 1 or 2) describing the length
of tick marks as a fraction of the smallest dimension. If
|
tckMinor |
numeric vector (length 1 or 2) describing the length of tick marks as a fraction of the smallest dimension. These tick marks can not be automatically labelled. If given a two-element vector, the first element describes the tick marks on the x-axis and the second element describes those on the y-axis. |
... |
additional |
This function plots a PolySet, where each unique (PID,
SID) describes a polygon. It connects each polygon's last
vertex to its first. The function supports both borders
(border, lty) and fills (col, density,
angle). When supplied with the appropriate arguments, it can
draw only borders or only fills. Unlike plotMap, it
ignores the aspect ratio. It clips polys to xlim and
ylim before plotting.
This function creates a blank plot when polys equals
NULL. In this case, the user must supply both xlim and
ylim arguments. Alternatively, it accepts the argument
type = "n" as part of ..., which is equivalent to specifying
polys = NULL, but requires a PolySet. In both cases,
the function's behaviour changes slightly. To resemble the
plot function, it plots the border, labels, and other
parts according to par parameters such as col.
For additional help on the arguments border, lty,
col, density, and angle, please see
polygon and par.
PolyData consisting of the PolyProps used to create the plot.
To satisfy the aspect ratio, this plotting routine resizes the plot
region. Consequently, par parameters such as
plt, mai, and mar will change. When the function
terminates, these changes persist to allow for additions to the plot.
Nicholas M. Boers, Staff Software Engineer
Jobber, Edmonton AB
Last modified Rd: 2013-04-10
addLabels,
addPolys,
addStipples,
clipPolys,
closePolys,
fixBound,
fixPOS,
locatePolys,
plotLines,
plotMap,
plotPoints,
thinPolys,
thickenPolys.
local(envir=.PBSmapEnv,expr={ oldpar = par(no.readonly=TRUE) #--- create a PolySet to plot polys <- data.frame(PID=rep(1,4),POS=1:4,X=c(0,1,1,0),Y=c(0,0,1,1)) #--- plot the PolySet plotPolys(polys, xlim=c(-.5,1.5), ylim=c(-.5,1.5), density=0) par(oldpar) })local(envir=.PBSmapEnv,expr={ oldpar = par(no.readonly=TRUE) #--- create a PolySet to plot polys <- data.frame(PID=rep(1,4),POS=1:4,X=c(0,1,1,0),Y=c(0,0,1,1)) #--- plot the PolySet plotPolys(polys, xlim=c(-.5,1.5), ylim=c(-.5,1.5), density=0) par(oldpar) })
A PolyData object comprises a data frame that summarises information
for each polyline/polygon in a PolySet; each PolyData record is defined by
a unique PID or (PID, SID combination).
PBSmapping functions that expect PolyData will accept properly formatted data frames in their place (see 'Details').
as.PolyData attempts to coerce a data frame to an object with
class PolyData.
is.PolyData returns TRUE if its argument is of class
PolyData.
as.PolyData(x, projection = NULL, zone = NULL) is.PolyData(x, fullValidation = TRUE)as.PolyData(x, projection = NULL, zone = NULL) is.PolyData(x, fullValidation = TRUE)
x |
data frame to be coerced or tested. |
projection |
optional |
zone |
optional |
fullValidation |
Boolean value; if |
We define PolyData as a data frame with a first column named PID
and (optionally) a second column named SID. Unlike a
PolySet, where each contour has many records corresponding
to the vertices, a PolyData object must have only one record for each
PID or each (PID, SID) combination. Conceptually,
this object associates data with contours, where the data correspond to
additional fields in the data frame. The R language conveniently
allows data frames to contain fields of various atomic modes
("logical", "numeric", "complex",
"character", and "null"). For example, PolyData with the
fields (PID, PName) might assign character names to a set
of primary polygons. Additionally, if fields X and Y exist
(perhaps representing locations for placing labels), consider adding
attributes zone and projection. Inserting the string
"PolyData" as the class attribute's first element alters the
behaviour of some functions, including print (if
PBSprint is TRUE) and summary.
Our software particularly uses PolyData to set various plotting
characteristics. Consistent with graphical parameters used by the R/S
functions lines and polygon, column names
can specify graphical properties:
lty - line type in drawing the border and/or shading
lines;
col - line or fill colour;
border - border colour;
density - density of shading lines;
angle - angle of shading lines.
When drawing polylines (as opposed to closed polygons), only lty
and col have meaning.
The as.PolyData method returns an object with classes
"PolyData" and "data.frame", in that order.
Nicholas M. Boers, Staff Software Engineer
Jobber, Edmonton AB
Last modified Rd: 2015-04-23
PolySet, EventData, LocationSet
A PolySet object comprises a data frame that defines a collection of polygonal contours (i.e., line segments joined at vertices). These contours can be open-ended (polylines) or closed (polygons).
PBSmapping functions that expect PolySet's will accept properly formatted data frames in their place (see 'Details').
as.PolySet attempts to coerce a data frame to an object with
class PolySet.
is.PolySet returns TRUE if its argument is of class
PolySet.
as.PolySet(x, projection = NULL, zone = NULL) is.PolySet(x, fullValidation = TRUE)as.PolySet(x, projection = NULL, zone = NULL) is.PolySet(x, fullValidation = TRUE)
x |
data frame to be coerced or tested. |
projection |
optional |
zone |
optional |
fullValidation |
Boolean value; if |
In our software, a PolySet data frame defines a collection of polygonal contours (i.e., line segments joined at vertices), based on four or five numerical fields:
PID - the primary identification number for a contour;
SID - optional, the secondary identification number for a
contour;
POS - the position number associated with a vertex;
X - the horizontal coordinate at a vertex;
Y - the vertical coordinate at a vertex.
The simplest PolySet lacks an SID column, and each PID
corresponds to a different contour. By analogy with a child's
“follow the dots” game, the POS field enumerates the
vertices to be connected by straight lines. Coordinates (X,
Y) specify the location of each vertex. Thus, in familiar
mathematical notation, a contour consists of points
() with , where corresponds to the
POS index. A PolySet has two potential interpretations. The first
associates a line segment with each successive pair of points from 1 to
, giving a polyline (in GIS terminology) composed of the
sequential segments. The second includes a final line segment joining
points and 1, thus giving a polygon.
The secondary ID field allows us to define regions as composites of
polygons. From this point of view, each primary ID identifies a
collection of polygons distinguished by secondary IDs. For example, a
single management area (PID) might consist of two fishing areas,
each defined by a unique SID. A secondary polygon can also
correspond to an inner boundary, like the hole in a doughnut. We adopt
the convention that POS goes from 1 to along an outer
boundary, but from to 1 along an inner boundary, regardless of
rotational direction. This contrasts with other GIS software, such as
ArcView (ESRI 1996), in which outer and inner boundaries correspond to
clockwise and counter-clockwise directions, respectively.
The SID field in a PolySet with secondary IDs must have integer values
that appear in ascending order for a given PID. Furthermore,
inner boundaries must follow the outer boundary that encloses them. The
POS field for each contour (PID, SID) must
similarly appear as integers in strictly increasing or decreasing order,
for outer and inner boundaries respectively. If the POS field
erroneously contains floating-point numbers, fixPOS can
renumber them as sequential integers, thus simplifying the insertion of
a new point, such as point 3.5 between points 3 and 4.
A PolySet can have a projection attribute, which may be missing,
that specifies a map projection. In the current version of PBS Mapping,
projection can have character values "LL" or "UTM",
referring to “Longitude-Latitude” and “Universal
Transverse Mercator”. We explain these projections more completely
below. If projection is numeric, it specifies the aspect ratio ,
the number of units per unit. Thus, units of
on the graph occupy the same distance as one unit of
. Another optional attribute zone specifies the UTM zone
(if projection="UTM") or the preferred zone for conversion from
Longitude-Latitude (if projection="LL").
A data frame's class attribute by default contains the string
"data.frame". Inserting the string "PolySet" as the class
vector's first element alters the behaviour of some functions. For
example, the summary function will print details specific
to a PolySet. Also, when PBSprint is TRUE, the
print function will display a PolySet's summary rather than the contents
of the data frame.
The as.PolySet method returns an object with classes
"PolySet" and "data.frame", in that order.
Nicholas M. Boers, Staff Software Engineer
Jobber, Edmonton AB
Last modified Rd: 2022-09-06
Environmental Systems Research Institute (ESRI). (1996) ArcView GIS: The Geographic Information System for Everyone. ESRI Press, Redlands, California. 340 pp.
PolyData, EventData, LocationSet
This function displays information about a PBS Mapping object.
Functions 'summary.EventData', 'summary.LocationSet',
'summary.PolyData', and 'summary.PolySet'
produce an object with class 'summary.PBS'.
## S3 method for class 'EventData' print(x, ...) ## S3 method for class 'LocationSet' print(x, ...) ## S3 method for class 'PolyData' print(x, ...) ## S3 method for class 'PolySet' print(x, ...) ## S3 method for class 'summary.PBS' print(x, ...)## S3 method for class 'EventData' print(x, ...) ## S3 method for class 'LocationSet' print(x, ...) ## S3 method for class 'PolyData' print(x, ...) ## S3 method for class 'PolySet' print(x, ...) ## S3 method for class 'summary.PBS' print(x, ...)
x |
|
... |
|
Nicholas M. Boers, Staff Software Engineer
Jobber, Edmonton AB
Last modified Rd: 2019-03-14
In package PBSmapping:
Data structures:
EventData,
LocationSet,
PolyData,
PolySet
Functions:
PBSprint,
summary
local(envir=.PBSmapEnv,expr={ #--- load the data (if using R) if (!is.null(version$language) && (version$language=="R")) data(nepacLL,envir=.PBSmapEnv) #--- change to summary printing style PBSprint <- TRUE #--- print the PolySet print(nepacLL) })local(envir=.PBSmapEnv,expr={ #--- load the data (if using R) if (!is.null(version$language) && (version$language=="R")) data(nepacLL,envir=.PBSmapEnv) #--- change to summary printing style PBSprint <- TRUE #--- print the PolySet print(nepacLL) })
PolySet of shapes to prove Pythagoras' Theorem:
.
data(pythagoras)data(pythagoras)
4 column data frame: PID = primary polygon ID,
POS = position of each vertex within a given polyline, X
= X-coordinate, and Y = Y-coordinate. Attributes:
projection = 1.
In R, the data must be loaded using the data function.
Nicholas M. Boers, Staff Software Engineer
Jobber, Edmonton AB
Last modified Rd: 2022-09-06
An artificial construct to illustrate the proof of Pythagoras' Theorem using trigonometry.
addPolys,
plotPolys,
plotMap,
PolySet.
worldLL/worldLLhigh Data SetsRefocus the worldLL/worldLLhigh data sets, e.g., refocus
them so that Eastern Canada appears to the west of Western Europe.
refocusWorld (polys, xlim=NULL, ylim=NULL, clip.AN=TRUE)refocusWorld (polys, xlim=NULL, ylim=NULL, clip.AN=TRUE)
polys |
|
xlim |
|
ylim |
|
clip.AN |
|
This function accepts a PolySet containing one or more polygons with X-coordinates that collectively span approximately 360 degrees. The function effectively joins the PolySet into a cylinder and then splits it at an arbitrary longitude according to the user-specified limits. Modifications in the resulting PolySet are restricted to shifting X-coordinates by +/- multiples of 360 degrees, and instead of clipping polygons, the return value simply omits out-of-range polygons.
PolySet, likely a subset of the input PolySet, which
retains the same PID/SID values.
The Antarctic polygon is treated as a special case in that it is expanded longitudinally
by duplicating it to the West and East of the base polygon. The expanded Antarctica is
then clipped to the limits of the plot, or not if 'clip.NA=FALSE'.
Nicholas M. Boers, Staff Software Engineer
Jobber, Edmonton AB
Rowan Haigh, Program Head – Offshore Rockfish
Pacific Biological Station (PBS), Fisheries & Oceans Canada (DFO), Nanaimo BC
locus opus: Institute of Ocean Sciences (IOS), Sidney BC
Last modified Rd: 2018-10-26
local(envir=.PBSmapEnv,expr={ oldpar = par(no.readonly=TRUE) #--- load appropriate data data(worldLL,envir=.PBSmapEnv) #--- set limits xlim <- c(-100,25) ylim <- c(0,90) #--- refocus and plot the world polys <- refocusWorld(worldLL, xlim, ylim) plotMap(polys, xlim, ylim) par(oldpar) })local(envir=.PBSmapEnv,expr={ oldpar = par(no.readonly=TRUE) #--- load appropriate data data(worldLL,envir=.PBSmapEnv) #--- set limits xlim <- c(-100,25) ylim <- c(0,90) #--- refocus and plot the world polys <- refocusWorld(worldLL, xlim, ylim) plotMap(polys, xlim, ylim) par(oldpar) })
Convert RGB (red-green-blue) colours to RYB (red-yellow-blue) colours, and vice versa. Algorithm based on Sugita and Takahashi (2015, 2017)
RGB2RYB(RGBmat) RYB2RGB(RYBmat)RGB2RYB(RGBmat) RYB2RGB(RYBmat)
RGBmat |
|
RYBmat |
|
The RYB colour wheel is more commonly used by artists, and provides a more intuitive system when blending colours – red and yellow makes orange, yellow and blue makes green, blue and red makes purple. On the RYB colour wheel, red lies opposite green, but on the RGB colour wheel, red lies opposite cyan.
Matrix of RGB or RYB primary colour intensities, where rows are records and columns are primary colours.
Opposite colours calculated in RYB space (1-RYB) are not always what one expects.
For example the colour "purple", RGB {160, 32, 240}, might better be specified
as RGB {126, 0, 255} before converting to RYB and inverting.
Rowan Haigh, Program Head – Offshore Rockfish
Pacific Biological Station (PBS), Fisheries & Oceans Canada (DFO), Nanaimo BC
locus opus: Regional Headquarters, Vancouver BC
Last modified Rd: 2022-07-06
Sugita, J. and Takahashi, T. (2017) Computational RYB Color Model and its Applications. IIEEJ Transactions on Image Electronics and Visual Computing 5(2): 110-122.
Sugita, J. and Takahashi, T. (2015) RYB Color Compositing. International Workshop on Applications in Information Technology, October 8-10, 2015.
In package PBSmapping:addBubbles
In package grDevices:col2rgb,
rgb,
rgb2hsv
Rotate a PolySet (or an EventData set) clockwise by a specified angle around a fixed point.
rotatePolys(polys, angle=40, centroid=c(500,5700), proj.out, zone, xlim=c(-135,-121.5), ylim=c(47,56), plot=FALSE, keep.extra=FALSE, ...) rotateEvents(data, angle=40, centroid=c(500,5700), proj.out, zone, plot=FALSE, keep.extra=FALSE, ...)rotatePolys(polys, angle=40, centroid=c(500,5700), proj.out, zone, xlim=c(-135,-121.5), ylim=c(47,56), plot=FALSE, keep.extra=FALSE, ...) rotateEvents(data, angle=40, centroid=c(500,5700), proj.out, zone, plot=FALSE, keep.extra=FALSE, ...)
polys |
|
data |
|
angle |
|
centroid |
|
proj.out |
|
zone |
|
xlim |
|
ylim |
|
plot |
|
keep.extra |
|
... |
|
Map rotation returns coordinates that are no longer meaningful with respect to the original coordinate system.
When displaying rotated maps, the user might wish to turn off axis labels using xaxt="n" and yaxt="n".
Rotated PolySet or EventData set where 'X' and 'Y' are the rotated coordinates in the projection specified by 'proj.out'.
The returned object has an attribute list object named 'rotation' that contains:
angle – angle of clockwise rotation in degrees
radian – angle of rotation in radians: pi * (-angle)/180
centroid – fixed point in UTM coordinates (km) around which map is rotated in UTM projection
R – rotation matrix (2-dimensional)
xylim – list object to keep track of 'xlim', 'ylim' and a bounding box 'xybox'.
projection – projection of the rotated PolySet or EventData set
zone – zone of the rotated PolySet or EventData set
When keep.extra=TRUE, the returned object will contain additional fields calculated by the rotational algorithm:
(X0,Y0) – original coordinates of the input PolySet | EventData set
(uX0,uY0) – original coordinates converted to UTM (only if original projection is 'LL')
(aX,aY) – UTM coordinates adjusted by subtracting the UTM centroid
(tX,tY) – adjusted UTM coordinates transformed by multiplying the rotational matrix
(rX,rY) – rotated UTM coordinates re-centered by adding the UTM centroid
Note:
If proj.out="UTM", the coordinates c(rX, rY) are used as the final rotated coordinates.
If proj.out="LL", the coordinates c(rX, rY) are transformed back into projection 'LL' as the final rotated coordinates.
Additionally, 'xylim' in the 'rotation' list attribute contains intermediary bounding box objects.
For instance, if the input PolySet | EventData object has projection 'LL', the 'xylim' object contains:
LL – original (X,Y) limits ('xlim', 'ylim', 'xybox')
UTM – original (X,Y) limits transformed to UTM coordinates
rot – rotated UTM (X,Y) limits
out – final projection (X,Y) limits
The map rotation algorithm is not heavily tested at this time. Report any issues to the package maintainer.
Rowan Haigh, Program Head – Offshore Rockfish
Pacific Biological Station (PBS), Fisheries & Oceans Canada (DFO), Nanaimo BC
locus opus: Institute of Ocean Sciences (IOS), Sidney BC
Last modified Rd: 2019-03-14
Academo – 2D Rotation about a point
In package PBSmapping:as.PolySet in PolySet,
clipPolys,
nepacLL,
plotMap,
plotPoints,
refocusWorld,
surveyData
summary method for PBS Mapping classes.
## S3 method for class 'EventData' summary(object, ...) ## S3 method for class 'LocationSet' summary(object, ...) ## S3 method for class 'PolyData' summary(object, ...) ## S3 method for class 'PolySet' summary(object, ...)## S3 method for class 'EventData' summary(object, ...) ## S3 method for class 'LocationSet' summary(object, ...) ## S3 method for class 'PolyData' summary(object, ...) ## S3 method for class 'PolySet' summary(object, ...)
object |
|
... |
|
After creating a list of summary statistics, this function assigns the
class 'summary.PBS' to the output in order to accomplish
formatted printing via print.summary.PBS.
A list of summary statistics.
Nicholas M. Boers, Staff Software Engineer
Jobber, Edmonton AB
Last modified Rd: 2019-03-14
In package PBSmapping:
EventData,
LocationSet,
PolyData,
PolySet,
PBSprint
local(envir=.PBSmapEnv,expr={ #--- load the data (if using R) if (!is.null(version$language) && (version$language=="R")) data(surveyData,envir=.PBSmapEnv) print(summary(surveyData)) })local(envir=.PBSmapEnv,expr={ #--- load the data (if using R) if (!is.null(version$language) && (version$language=="R")) data(surveyData,envir=.PBSmapEnv) print(summary(surveyData)) })
EventData of Pacific ocean perch (POP) tow information (1966-89).
data(surveyData)data(surveyData)
Data frame consisting of 9 columns: PID = primary polygon ID,
POS = position of each vertex within a given polygon, X
= longitude coordinate, Y = latitude coordinate, trip
= trip ID, tow = tow number in trip, catch = catch of
POP (kg), effort = tow effort (minutes), depth = fishing
depth (m), and year = year of survey trip. Attributes:
projection = "LL", zone = 9.
In R, the data must be loaded using the data function.
Nicholas M. Boers, Staff Software Engineer
Jobber, Edmonton AB
Last modified Rd: 2008-09-03
The GFBio database, maintained at the Pacific Biological Station (Fisheries and Oceans Canada, Nanaimo, BC V9T 6N7), archives catches and related biological data from commercial groundfish fishing trips and research/assessment cruises off the west coast of British Columbia (BC).
The POP (Sebastes alutus) survey data were extracted from GFBio. The data extraction covers bottom trawl surveys that focus primarily on POP biomass estimation: 1966-89 for the central BC coast and 1970-85 for the west coast of Vancouver Island. Additionally, a 1989 cruise along the entire BC coast concentrated on the collection of biological samples. Schnute et al. (2001) provide a more comprehensive history of POP surveys including the subset of data presented here.
Schnute, J.T., Haigh, R., Krishka, B.A. and Starr, P. (2001) Pacific ocean perch assessment for the west coast of Canada in 2001. Canadian Science Advisory Secretariat, Research Document 2001/138, 90 pp.
addPoints,
combineEvents,
EventData,
findPolys,
makeGrid,
plotPoints.
Thicken a PolySet, where each unique (PID, SID)
describes a polygon.
thickenPolys (polys, tol = 1, filter = 3, keepOrig = TRUE, close = TRUE)thickenPolys (polys, tol = 1, filter = 3, keepOrig = TRUE, close = TRUE)
polys |
PolySet to thicken. |
tol |
tolerance (in kilometres when |
filter |
minimum number of vertices per result polygon. |
keepOrig |
Boolean value; if |
close |
Boolean value; if |
This function thickens each polygon within polys according to
the input arguments.
If keepOrig = TRUE, all of the original vertices appear in the
result. It calculates the distance between two sequential original
vertices, and if that distance exceeds tol, it adds a
sufficient number of vertices spaced evenly between the two original
vertices so that the distance between vertices no longer exceeds
tol. If close = TRUE, it adds intermediate vertices
between the last and first vertices when necessary.
If keepOrig = FALSE, only the first vertex of each polygon is
guaranteed to appear in the results. From this first vertex, the
algorithm walks the polygon summing the distance between vertices.
When this cumulative distance exceeds tol, it adds a vertex on
the line segment under inspection. After doing so, it resets the
distance sum, and walks the polygon from this new vertex. If
close = TRUE, it will walk the line segment from the last
vertex to the first.
PolySet containing the thickened data. The function
recalculates the POS values for each polygon.
Nicholas M. Boers, Staff Software Engineer
Jobber, Edmonton AB
Last modified Rd: 2013-04-10
local(envir=.PBSmapEnv,expr={ oldpar = par(no.readonly=TRUE) #--- load the data (if using R) if (!is.null(version$language) && (version$language=="R")) data(nepacLL,envir=.PBSmapEnv) #--- plot Vancouver Island plotMap(nepacLL[nepacLL$PID == 33, ]) #--- calculate a thickened version using a 30 kilometres tolerance, #--- without keeping the original points p <- thickenPolys(nepacLL[nepacLL$PID == 33, ], tol = 30, keepOrig = FALSE) #--- convert the PolySet to EventData by dropping the PID column and #--- renaming POS to EID p <- p[-1]; names(p)[1] <- "EID" #--- convert the now invalid PolySet into a data frame, and then into #--- EventData p <- as.EventData(as.data.frame(p), projection="LL") #--- plot the results addPoints(p, col=2, pch=19) par(oldpar) })local(envir=.PBSmapEnv,expr={ oldpar = par(no.readonly=TRUE) #--- load the data (if using R) if (!is.null(version$language) && (version$language=="R")) data(nepacLL,envir=.PBSmapEnv) #--- plot Vancouver Island plotMap(nepacLL[nepacLL$PID == 33, ]) #--- calculate a thickened version using a 30 kilometres tolerance, #--- without keeping the original points p <- thickenPolys(nepacLL[nepacLL$PID == 33, ], tol = 30, keepOrig = FALSE) #--- convert the PolySet to EventData by dropping the PID column and #--- renaming POS to EID p <- p[-1]; names(p)[1] <- "EID" #--- convert the now invalid PolySet into a data frame, and then into #--- EventData p <- as.EventData(as.data.frame(p), projection="LL") #--- plot the results addPoints(p, col=2, pch=19) par(oldpar) })
Thin a PolySet, where each unique (PID, SID)
describes a polygon.
thinPolys (polys, tol = 1, filter = 3)thinPolys (polys, tol = 1, filter = 3)
polys |
PolySet to thin. |
tol |
tolerance (in kilometres when |
filter |
minimum number of vertices per result polygon. |
This function executes the Douglas-Peuker line simplification
algorithm on each polygon within polys.
PolySet containing the thinned data. The function recalculates
the POS values for each polygon.
Nicholas M. Boers, Staff Software Engineer
Jobber, Edmonton AB
Last modified Rd: 2013-04-10
local(envir=.PBSmapEnv,expr={ oldpar = par(no.readonly=TRUE) #--- load the data (if using R) if (!is.null(version$language) && (version$language=="R")) data(nepacLL,envir=.PBSmapEnv) #--- plot a thinned version of Vancouver Island (3 km tolerance) plotMap(thinPolys(nepacLL[nepacLL$PID == 33, ], tol = 3)) #--- add the original Vancouver Island in a different line type to #--- emphasize the difference addPolys(nepacLL[nepacLL$PID == 33, ], border=2, lty=8, density=0) par(oldpar) })local(envir=.PBSmapEnv,expr={ oldpar = par(no.readonly=TRUE) #--- load the data (if using R) if (!is.null(version$language) && (version$language=="R")) data(nepacLL,envir=.PBSmapEnv) #--- plot a thinned version of Vancouver Island (3 km tolerance) plotMap(thinPolys(nepacLL[nepacLL$PID == 33, ], tol = 3)) #--- add the original Vancouver Island in a different line type to #--- emphasize the difference addPolys(nepacLL[nepacLL$PID == 33, ], border=2, lty=8, density=0) par(oldpar) })
PolyData of tow information for a longspine thornyhead survey (2001).
data(towData)data(towData)
Data frame consisting of 8 columns: PID = primary polygon ID,
POS = position of each vertex within a given polygon, X
= longitude coordinate, Y = latitude coordinate, depth
= fishing depth (m), effort = tow effort (minutes), distance
= tow track distance (km), catch = catch of longspine
thornyhead (kg), and year = year of survey. Attributes:
projection = "LL", zone = 9.
In R, the data must be loaded using the data function.
Nicholas M. Boers, Staff Software Engineer
Jobber, Edmonton AB
Last modified Rd: 2008-09-03
The GFBio database, maintained at the Pacific Biological Station (Fisheries and Oceans Canada, Nanaimo, BC V9T 6N7), archives catches and related biological data from commercial groundfish fishing trips and research/assessment cruises off the west coast of British Columbia (BC). The longspine thornyhead (Sebastolobus altivelis) survey data were extracted from GFBio. Information on the first 45 tows from the 2001 survey (Starr et al. 2002) are included here. Effort is time (minutes) from winch lock-up to winch release.
Starr, P.J., Krishka, B.A. and Choromanski, E.M. (2002) Trawl survey for thornyhead biomass estimation off the west coast of Vancouver Island, September 15 - October 2, 2001. Canadian Technical Report of Fisheries and Aquatic Sciences 2421, 60 pp.
makeProps,
PolyData,
towTracks.
PolySet of geo-referenced polyline tow track data from a longspine thornyhead survey (2001).
data(towTracks)data(towTracks)
Data frame consisting of 4 columns: PID = primary polygon ID,
POS = position of each vertex within a given polyline, X
= longitude coordinate, and Y = latitude
coordinate. Attributes: projection = "LL", zone = 9.
In R, the data must be loaded using the data function.
Nicholas M. Boers, Staff Software Engineer
Jobber, Edmonton AB
Last modified Rd: 2008-09-03
The longspine thornyhead (Sebastolobus altivelis) tow track spatial coordinates are available at the Pacific Biological Station (Fisheries and Oceans Canada, Nanaimo, BC V9T 6N7). The geo-referenced coordinates of the first 45 tows from the 2001 survey (Starr et al. 2002) are included here. Coordinates are recorded once per minute between winch lock-up and winch release.
Starr, P.J., Krishka, B.A. and Choromanski, E.M. (2002) Trawl survey for thornyhead biomass estimation off the west coast of Vancouver Island, September 15 - October 2, 2001. Canadian Technical Report of Fisheries and Aquatic Sciences 2421, 60 pp.
addLines,
calcLength,
clipLines,
plotLines,
PolySet,
towData.