Package 'rgeos'

Converts a bounding box into a SpatialPolygons object.

Description

Converts a bounding box into a SpatialPolygons object.

Usage

bbox2SP(n,s,w,e,bbox=NA,proj4string=CRS("+init=epsg:4326"))

Arguments

n	the top north latitude
s	the bottom south latitude
w	the most western longitude
e	the most eastern longitude
bbox	a bounding box 2 x 2 matrix as produced by bbox
proj4string	a coordinate reference system as defined in CRS

Details

This function converts a set of coordinates limiting a bounding box into a SpatialPolygons. It can be used for instance to clip a subset of a larger spatial object (e.g. using gIntersection)

Value

An object of SpatialPolygons class.

Examples

library(sp)
run <- FALSE
if (require(rgdal, quietly=TRUE)) run <- TRUE
if (run) {
cities <- readOGR(dsn=system.file("vectors", package = "rgdal")[1], layer="cities")
n<-75
s<-30
w<--40
e<-32
myPoly<-bbox2SP(n,s,e,w)
}
if (run) {
plot(cities)
plot(myPoly,border="red",add=TRUE)
}
if (run) {
bb<-bbox(cities)
myPoly<-bbox2SP(bbox=bb,proj4string=CRS(proj4string(cities)))
plot(myPoly,add=TRUE,border="blue")
}

---

Area of Geometry

Description

Calculates the area of the given geometry.

Usage

gArea(spgeom, byid=FALSE)

Arguments

spgeom	sp object as defined in package sp
byid	Logical determining if the function should be applied across subgeometries (TRUE) or the entire object (FALSE)

Value

Returns the area of the geometry in the units of the current projection. By definition non-[MULTI]POLYGON geometries have an area of 0. The area of a POLYGON is the area of its shell less the area of any holes. Note that this value may be different from the area slot of the Polygons class as this value does not subtract the area of any holes in the geometry.

Author(s)

Roger Bivand & Colin Rundel

See Also

gLength

Examples

gArea(readWKT("POINT(1 1)"))
	gArea(readWKT("LINESTRING(0 0,1 1,2 2)"))
	gArea(readWKT("LINEARRING(0 0,3 0,3 3,0 3,0 0)"))
	
	
	p1 = readWKT("POLYGON((0 0,3 0,3 3,0 3,0 0))")
	p2 = readWKT("POLYGON((0 0,3 0,3 3,0 3,0 0),(1 1,2 1,2 2,1 2,1 1))")
	
	gArea(p1)
	p1@polygons[[1]]@area
	
	gArea(p2)
	p2@polygons[[1]]@area

---

Boundary of Geometry

Description

Function for determinging the Boundary of the given geometry as defined by SFS Section 2.1.13.1

Usage

gBoundary(spgeom, byid=FALSE, id = NULL)

Arguments

spgeom	sp object as defined in package sp
byid	Logical determining if the function should be applied across subgeometries (TRUE) or the entire object (FALSE)
id	Character vector defining id labels for the resulting geometries, if unspecified returned geometries will be labeled based on their parent geometries' labels.

Value

Depending of the class of the spgeom the returned results will differ.

Based on the documentation of JTS (on which GEOS is based) the following outputs are expected:

Point	empty GeometryCollection
MultiPoint	empty GeometryCollection
LineString	if closed: empty MultiPoint if not closed: MultiPoint containing the two endpoints.
MultiLineString	MultiPoint obtained by applying the Mod-2 rule to the boundaries of the element LineStrings
LinearRing	empty MultiPoint
Polygon	MultiLineString containing the LinearRings of the shell and holes, in that order (SFS 2.1.10)
MultiPolygon	MultiLineString containing the LinearRings for the boundaries of the element polygons, in the same order as they occur in the MultiPolygon (SFS 2.1.12/JTS)
GeometryCollection	The boundary of an arbitrary collection of geometries whose interiors are disjoint consist of geometries drawn from the boundaries of the element geometries by application of the Mod-2 rule (SFS Section 2.1.13.1)

The mod-2 rule states that for a multiline a point is on the boundary if and only if it on the boundary of an odd number of subgeometries of the multiline (See example below).

Author(s)

Roger Bivand & Colin Rundel

See Also

gCentroid gConvexHull gEnvelope gPointOnSurface

Examples

x = readWKT("POLYGON((0 0,10 0,10 10,0 10,0 0))")
b = gBoundary(x)

plot(x,col='black')
plot(b,col='red',lwd=3,add=TRUE)

# mod-2 rule example
x1 = readWKT("MULTILINESTRING((2 2,2 0),(2 2,0 2))")
x2 = readWKT("MULTILINESTRING((2 2,2 0),(2 2,0 2),(2 2,4 2))")
x3 = readWKT("MULTILINESTRING((2 2,2 0),(2 2,0 2),(2 2,4 2),(2 2,2 4))")
x4 = readWKT("MULTILINESTRING((2 2,2 0),(2 2,0 2),(2 2,4 2),(2 2,2 4),(2 2,4 4))")

b1 = gBoundary(x1)
b2 = gBoundary(x2)
b3 = gBoundary(x3)
b4 = gBoundary(x4)

par(mfrow=c(2,2))
plot(x1); plot(b1,pch=16,col='red',add=TRUE)
plot(x2); plot(b2,pch=16,col='red',add=TRUE)
plot(x3); plot(b3,pch=16,col='red',add=TRUE)
plot(x4); plot(b4,pch=16,col='red',add=TRUE)

---

Centroid of Geometry

Description

Function calculates the centroid of the given geometry.

Usage

gCentroid(spgeom, byid=FALSE, id = NULL)

Arguments

spgeom	sp object as defined in package sp
byid	Logical determining if the function should be applied across subgeometries (TRUE) or the entire object (FALSE)
id	Character vector defining id labels for the resulting geometries, if unspecified returned geometries will be labeled based on their parent geometries' labels.

Details

Returns a SpatialPoints object of the centroid(s) for spgeom.

Author(s)

Roger Bivand & Colin Rundel

See Also

gBoundary gConvexHull gEnvelope gPointOnSurface

Examples

x = readWKT(paste("GEOMETRYCOLLECTION(POLYGON((0 0,10 0,10 10,0 10,0 0)),",
 "POLYGON((15 0,25 15,35 0,15 0)))"))

# Centroids of both the square and circle independently
c1 = gCentroid(x,byid=TRUE) 
# Centroid of square and circle together
c2 = gCentroid(x)

plot(x)
plot(c1,col='red',add=TRUE)
plot(c2,col='blue',add=TRUE)

---

Geometry Relationships - Contains and Within

Description

Functions for testing whether one geometry contains or is contained within another geometry

Usage

gContains(spgeom1, spgeom2 = NULL, byid = FALSE, prepared=TRUE,
 returnDense=TRUE, STRsubset=FALSE, checkValidity=FALSE)
gContainsProperly(spgeom1, spgeom2 = NULL, byid = FALSE, returnDense=TRUE,
 checkValidity=FALSE)
gCovers(spgeom1, spgeom2 = NULL, byid = FALSE, returnDense=TRUE, checkValidity=FALSE)
gCoveredBy(spgeom1, spgeom2 = NULL, byid = FALSE, returnDense=TRUE, checkValidity=FALSE)
gWithin(spgeom1, spgeom2 = NULL, byid = FALSE, returnDense=TRUE, checkValidity=FALSE)

Arguments

spgeom1, spgeom2	sp objects as defined in package sp. If spgeom2 is NULL then spgeom1 is compared to itself.
byid	Logical vector determining if the function should be applied across ids (TRUE) or the entire object (FALSE) for spgeom1 and spgeom2
prepared	Logical determining if prepared geometry (spatially indexed) version of the GEOS function should be used. In general prepared geometries should be faster than the alternative.
returnDense	default TRUE, if false returns a list of the length of spgeom1 of integer vectors listing the 1:length(spgeom2) indices which would be TRUE in the dense logical matrix representation; useful when the sizes of the byid=TRUE returned matrix is very large and it is sparse; essential when the returned matrix would be too large
checkValidity	default FALSE; error meesages from GEOS do not say clearly which object fails if a topology exception is encountered. If this argument is TRUE, gIsValid is run on each in turn in an environment in which object names are available. If objects are invalid, this is reported and those affected are named
STRsubset	logical argument for future use

Value

gContains returns TRUE if none of the point of spgeom2 is outside of spgeom1 and at least one point of spgeom2 falls within spgeom1.

gContainsProperly returns TRUE under the same conditions as gContains with the additional requirement that spgeom2 does not intersect with the boundary of spgeom1. As such any given geometry will Contain itself but will not ContainProperly itself.

gCovers returns TRUE if no point in spgeom2 is outside of spgeom1. This is slightly different from gContains as it does not require a point within spgeom1 which can be an issue as boundaries are not considered to be "within" a geometry, see gBoundary for specifics of geometry boundaries.

gCoveredBy is the converse of gCovers and is equivalent to swapping spgeom1 and spgeom2.

gWithin is the converse of gContains and is equivalent to swapping spgeom1 and spgeom2.

Note

Error messages from GEOS, in particular topology exceptions, report 0-based object order, so geom 0 is spgeom1, and geom 1 is spgeom2.

Author(s)

Roger Bivand & Colin Rundel

References

Helpful information on the subtle differences between these functions: http://lin-ear-th-inking.blogspot.com/2007/06/subtleties-of-ogc-covers-spatial.html

See Also

gCrosses gDisjoint gEquals gEqualsExact gIntersects gOverlaps gRelate gTouches

Examples

l1 = readWKT("LINESTRING(0 3,1 1,2 2,3 0)")
	l2 = readWKT("LINESTRING(1 3.5,3 3,2 1)")
	l3 = readWKT("LINESTRING(1 3.5,3 3,4 1)")
	
	pt1 = readWKT("MULTIPOINT(1 1,3 0)")
	pt2 = readWKT("MULTIPOINT(0 3,3 0)")
	pt3 = readWKT("MULTIPOINT(1 1,2 2,3 1)")
	
	p1 = readWKT("POLYGON((0 0,0 2,1 3.5,3 3,4 1,3 0,0 0))")
	p2 = readWKT("POLYGON((1 1,1 2,2 2,2 1,1 1))")

	
	par(mfrow=c(2,3))
	plot(l1,col='blue');plot(pt1,add=TRUE,pch=16)
	title(paste("Contains:",gContains(l1,pt1),
				"\nContainsProperly:",gContainsProperly(l1,pt1),
				"\nCovers:",gCovers(l1,pt1)))
	
	plot(l1,col='blue');plot(pt2,add=TRUE,pch=16)
	title(paste("Contains:",gContains(l1,pt2),
				"\nContainsProperly:",gContainsProperly(l1,pt2),
				"\nCovers:",gCovers(l1,pt2)))
	
	plot(p1,col='blue',border='blue');plot(pt3,add=TRUE,pch=16)
	title(paste("Contains:",gContains(p1,pt3),
				"\nContainsProperly:",gContainsProperly(p1,pt3),
				"\nCovers:",gCovers(p1,pt3)))
	
	plot(p1,col='blue',border='blue');plot(l2,lwd=2,add=TRUE,pch=16)
	title(paste("Contains:",gContains(p1,l2),
				"\nContainsProperly:",gContainsProperly(p1,l2),
				"\nCovers:",gCovers(p1,l2)))

	plot(p1,col='blue',border='blue');plot(l3,lwd=2,add=TRUE,pch=16)
	title(paste("Contains:",gContains(p1,l3),
				"\nContainsProperly:",gContainsProperly(p1,l3),
				"\nCovers:",gCovers(p1,l3)))

	plot(p1,col='blue',border='blue');plot(p2,col='black',add=TRUE,pch=16)
	title(paste("Contains:",gContains(p1,p2),
				"\nContainsProperly:",gContainsProperly(p1,p2),
				"\nCovers:",gCovers(p1,p2)))

---

Convex Hull of Geometry

Description

Function produces the Convex Hull of the given geometry, the smallest convex polygon that contains all subgeometries

Usage

gConvexHull(spgeom, byid=FALSE, id = NULL)

Arguments

spgeom	sp object as defined in package sp
byid	Logical determining if the function should be applied across subgeometries (TRUE) or the entire object (FALSE)
id	Character vector defining id labels for the resulting geometries, if unspecified returned geometries will be labeled based on their parent geometries' labels.

Details

Returns the convex hull as a SpatialPolygons object.

Author(s)

Roger Bivand & Colin Rundel

See Also

gBoundary gCentroid gEnvelope gPointOnSurface

Examples

x = readWKT(paste("POLYGON((0 40,10 50,0 60,40 60,40 100,50 90,60 100,60",
 "60,100 60,90 50,100 40,60 40,60 0,50 10,40 0,40 40,0 40))"))

ch = gConvexHull(x)

plot(x,col='blue',border='blue')
plot(ch,add=TRUE)

---

Coverage Union

Description

GEOSCoverageUnion is an optimized union algorithm for polygonal inputs that are correctly noded and do not overlap. It will not generate an error (return NULL) for inputs that do not satisfy this constraint.

Usage

gCoverageUnion(spgeom, byid=FALSE, id = NULL)

Arguments

spgeom	sp object as defined in package sp
byid	Logical determining if the function should be applied across subgeometries (TRUE) or the entire object (FALSE)
id	Character vector defining id labels for the resulting geometries, if unspecified returned geometries will be labeled based on their parent geometries' labels.

Author(s)

Roger Bivand

See Also

gUnaryUnion

Examples

run <- FALSE
if (require(maptools)) run <- TRUE
if (run) {
    nc1 <- readShapePoly(system.file("shapes/sids.shp", package="maptools")[1],
        proj4string=CRS("+proj=longlat +datum=NAD27"))
}
if (run) {
    print(system.time(oU <- gUnionCascaded(nc1)))
    if (version_GEOS0() >= "3.8.0") {
        print(system.time(oCU <- gCoverageUnion(nc1)))
    }
}

---

Geometry Relationships - Crosses and Overlaps

Description

Functions for testing whether geometries share some but not all interior points

Usage

gCrosses(spgeom1, spgeom2 = NULL, byid = FALSE, returnDense=TRUE,
 checkValidity=FALSE)
	gOverlaps(spgeom1, spgeom2 = NULL, byid = FALSE, returnDense=TRUE,
 checkValidity=FALSE)

Arguments

spgeom1, spgeom2	sp objects as defined in package sp. If spgeom2 is NULL then spgeom1 is compared to itself.
byid	Logical vector determining if the function should be applied across ids (TRUE) or the entire object (FALSE) for spgeom1 and spgeom2
returnDense	default TRUE, if false returns a list of the length of spgeom1 of integer vectors listing the 1:length(spgeom2) indices which would be TRUE in the dense logical matrix representation; useful when the sizes of the byid=TRUE returned matrix is very large and it is sparse; essential when the returned matrix would be too large
checkValidity	default FALSE; error meesages from GEOS do not say clearly which object fails if a topology exception is encountered. If this argument is TRUE, gIsValid is run on each in turn in an environment in which object names are available. If objects are invalid, this is reported and those affected are named

Value

gCrosses returns TRUE when the geometries share some but not all interior points, and the dimension of the intersection is less than that of at least one of the geometries.

gOverlaps returns TRUE when the geometries share some but not all interior points, and the intersection has the same dimension as the geometries themselves.

Note

Error messages from GEOS, in particular topology exceptions, report 0-based object order, so geom 0 is spgeom1, and geom 1 is spgeom2.

Author(s)

Roger Bivand & Colin Rundel

See Also

gContains gContainsProperly gCovers gCoveredBy gDisjoint gEquals gEqualsExact gIntersects gRelate gTouches gWithin

Examples

l1 = readWKT("LINESTRING(0 3,1 1,2 2,3 0)")
	l2 = readWKT("LINESTRING(0 0.5,1 1,2 2,3 2.5)")
	l3 = readWKT("LINESTRING(1 3,1.5 1,2.5 2)")
	
	pt1 = readWKT("MULTIPOINT(1 1,3 0)")
	pt2 = readWKT("MULTIPOINT(1 1,3 0,1 2)")

	p1 = readWKT("POLYGON((0 0,0 2,1 3.5,3 3,4 1,3 0,0 0))")
	p2 = readWKT("POLYGON((2 2,3 4,4 1,4 0,2 2))")
	
	par(mfrow=c(2,3))
	plot(l1,col='blue');plot(pt1,add=TRUE,pch=16)
	title(paste("Crosses:",gCrosses(l1,pt1),
				"\nOverlaps:",gOverlaps(l1,pt1)))
	
	plot(l1,col='blue');plot(pt2,add=TRUE,pch=16)
	title(paste("Crosses:",gCrosses(l1,pt2),
				"\nOverlaps:",gOverlaps(l1,pt2)))
				
	plot(l1,col='blue');plot(l2,add=TRUE)
	title(paste("Crosses:",gCrosses(l1,l2),
				"\nOverlaps:",gOverlaps(l1,l2)))
				
	plot(l1,col='blue');plot(l3,add=TRUE)
	title(paste("Crosses:",gCrosses(l1,l3),
				"\nOverlaps:",gOverlaps(l1,l3)))
				
	plot(p1,border='blue',col='blue');plot(l1,add=TRUE)
	title(paste("Crosses:",gCrosses(p1,l1),
				"\nOverlaps:",gOverlaps(p1,l1)))

	plot(p1,border='blue',col='blue');plot(p2,add=TRUE)
	title(paste("Crosses:",gCrosses(p1,p2),
				"\nOverlaps:",gOverlaps(p1,p2)))

---

Compute Delaunay triangulation between points

Description

Function to compute the Delaunay triangulation between points; only available for GEOS >= 3.4.0.

Usage

gDelaunayTriangulation(spgeom, tolerance=0.0, onlyEdges=FALSE)

Arguments

spgeom	sp points object as defined in package sp
tolerance	Numerical tolerance value to be used in triangulation
onlyEdges	Logical, default returns triangles as polygons, if TRUE, returns a SpatialLines object with a single MULTILINESTRING

Details

When onlyEdges is TRUE, the SpatialLines object may be de-merged to identify the input points that are touched by each edge, making it possible to identify spatial neighbours.

Value

Either a SpatialPolygons object or a SpatialLines object containing a single Lines object of the undirected edges in the triangulation.

Author(s)

Roger Bivand

References

https://en.wikipedia.org/wiki/Delaunay_triangulation

Examples

if (version_GEOS0() > "3.4.0") {
library(sp)
data(meuse)
coordinates(meuse) <- c("x", "y")
plot(gDelaunayTriangulation(meuse))
points(meuse)
out <- gDelaunayTriangulation(meuse, onlyEdges=TRUE)
lns <- slot(slot(out, "lines")[[1]], "Lines")
out1 <- SpatialLines(lapply(seq(along=lns), function(i) Lines(list(lns[[i]]),
 ID=as.character(i))))
out2 <- lapply(1:length(out1), function(i) which(gTouches(meuse, out1[i],
 byid=TRUE)))
out3 <- do.call("rbind", out2)
o <- order(out3[,1], out3[,2])
out4 <- out3[o,]
out5 <- data.frame(from=out4[,1], to=out4[,2], weight=1)
head(out5)
## Not run: 
if (require(spdep)) {
class(out5) <- c("spatial.neighbour", class(out5))
attr(out5, "n") <- length(meuse)
attr(out5, "region.id") <- as.character(1:length(meuse))
nb1 <- sn2listw(out5)$neighbours
nb2 <- make.sym.nb(nb1)
}

## End(Not run)
}

---

Distance between geometries

Description

Calculates the distance between the given geometries

Usage

gDistance(spgeom1, spgeom2=NULL, byid=FALSE, hausdorff=FALSE, densifyFrac = NULL)
gWithinDistance(spgeom1, spgeom2=NULL, dist, byid=FALSE,
 hausdorff=FALSE, densifyFrac=NULL)

Arguments

spgeom1, spgeom2	sp objects as defined in package sp. If spgeom2 is NULL then spgeom1 is compared to itself.
byid	Logical vector determining if the function should be applied across ids (TRUE) or the entire object (FALSE) for spgeom1 and spgeom2
hausdorff	Logical determining if the discrete Hausdorff distance should be calculated
densifyFrac	Numerical value between 0 and 1 that determines the fraction by which to densify each segment of the geometry.
dist	Numerical value that determines cutoff distance

Details

Discrete Hausdorff distance is essentially a measure of the similarity or dissimilarity of the two geometries, see references below for more detailed explanations / descriptions.

If hausdorff is TRUE and densifyFrac is specified then the geometries' segments are densified by splitting each segment into equal length subsegments whose fraction of the total length is equal to densifyFrac.

Value

gDistance by default returns the cartesian minimum distance between the two geometries in the units of the current projection. If hausdorff is TRUE then the Hausdorff distance is returned for the two geometries.

gWithinDistance returns TRUE if returned distance is less than or equal to the specified dist.

Author(s)

Roger Bivand & Colin Rundel

References

Hausdorff Differences: https://en.wikipedia.org/wiki/Hausdorff_distance http://lin-ear-th-inking.blogspot.com/2009/01/computing-geometric-similarity.html

See Also

gWithinDistance

Examples

pt1 = readWKT("POINT(0.5 0.5)")
pt2 = readWKT("POINT(2 2)")

p1 = readWKT("POLYGON((0 0,1 0,1 1,0 1,0 0))")
p2 = readWKT("POLYGON((2 0,3 1,4 0,2 0))")

gDistance(pt1,pt2)
gDistance(p1,pt1)
gDistance(p1,pt2)
gDistance(p1,p2)


p3 = readWKT("POLYGON((0 0,2 0,2 2,0 2,0 0))")
p4 = readWKT("POLYGON((0 0,2 0,2 1.9,1.9 2,0 2,0 0))")
p5 = readWKT("POLYGON((0 0,2 0,2 1.5,1.5 2,0 2,0 0))")
p6 = readWKT("POLYGON((0 0,2 0,2 1,1 2,0 2,0 0))")
p7 = readWKT("POLYGON((0 0,2 0,0 2,0 0))")

gDistance(p3,hausdorff=TRUE)
gDistance(p3,p4,hausdorff=TRUE)
gDistance(p3,p5,hausdorff=TRUE)
gDistance(p3,p6,hausdorff=TRUE)
gDistance(p3,p7,hausdorff=TRUE)

---

Envelope of Geometry

Description

Function calculates the rectangular bounding box for the given geometry

Usage

gEnvelope(spgeom, byid=FALSE, id = NULL)

Arguments

spgeom	sp object as defined in package sp
byid	Logical determining if the function should be applied across subgeometries (TRUE) or the entire object (FALSE)
id	Character vector defining id labels for the resulting geometries, if unspecified returned geometries will be labeled based on their parent geometries' labels.

Details

Returns the rectangular bounding box as a SpatialPolygons object. If spgeom is a degenerate case (horizontal/vertical line, single point) then the function may return an object with lower dimension (SpatialLines or SpatialPoints) or an invalid polygon.

Author(s)

Roger Bivand & Colin Rundel

See Also

gBoundary gCentroid gConvexHull gPointOnSurface

Examples

x = readWKT(paste("POLYGON((0 40,10 50,0 60,40 60,40 100,50 90,60 100,60",
 "60,100 60,90 50,100 40,60 40,60 0,50 10,40 0,40 40,0 40))"))
env = gEnvelope(x)

plot(x,col='blue',border='blue')
plot(env,add=TRUE)

#Degenerate Cases
gEnvelope(readWKT("POINT(1 1)")) #returns SpatialPoints
gEnvelope(readWKT("LINESTRING(1 1,1 2)")) #invalid polygon
gEnvelope(readWKT("LINESTRING(1 1,2 1)")) #invalid polygon

---

Geometry Relationships - Equality

Description

Function for testing equivalence of the given geometries

Usage

gEquals(spgeom1, spgeom2 = NULL, byid = FALSE, returnDense=TRUE,
 checkValidity=FALSE)
	gEqualsExact(spgeom1, spgeom2 = NULL, tol=0.0, byid = FALSE,
 returnDense=TRUE, checkValidity=FALSE)

Arguments

spgeom1, spgeom2	sp objects as defined in package sp. If spgeom2 is NULL then spgeom1 is compared to itself.
byid	Logical vector determining if the function should be applied across ids (TRUE) or the entire object (FALSE) for spgeom1 and spgeom2
tol	Numerical value of tolerance to use when assessing equivalence
returnDense	default TRUE, if false returns a list of the length of spgeom1 of integer vectors listing the 1:length(spgeom2) indices which would be TRUE in the dense logical matrix representation; useful when the sizes of the byid=TRUE returned matrix is very large and it is sparse; essential when the returned matrix would be too large
checkValidity	default FALSE; error meesages from GEOS do not say clearly which object fails if a topology exception is encountered. If this argument is TRUE, gIsValid is run on each in turn in an environment in which object names are available. If objects are invalid, this is reported and those affected are named

Value

gEquals returns TRUE if geometries are "spatially equivalent" which requires that spgeom1 is within spgeom2 and spgeom2 is within spgeom1, this ignores ordering of points within the geometries. Note that it is possible for geometries with different coordinates to be "spatially equivalent".

gEqualsExact returns TRUE if geometries are "exactly equivalent" which requires that spgeom1 and spgeom1 are "spatially equivalent" and that their constituent points are in the same order.

Note

Error messages from GEOS, in particular topology exceptions, report 0-based object order, so geom 0 is spgeom1, and geom 1 is spgeom2.

Author(s)

Roger Bivand & Colin Rundel

See Also

gContains gContainsProperly gCovers gCoveredBy gCrosses gDisjoint gEqualsExact gIntersects gOverlaps gRelate gTouches gWithin

Examples

# p1 and p2 are spatially identical but not exactly identical due to point ordering
p1=readWKT("POLYGON((0 0,1 0,1 1,0 1,0 0))")
p2=readWKT("POLYGON((1 1,0 1,0 0,1 0,1 1))")
p3=readWKT("POLYGON((0.01 0.01,1.01 0.01,1.01 1.01,0.01 1.01,0.01 0.01))")

gEquals(p1,p2)
gEquals(p1,p3)
gEqualsExact(p1,p2)
gEqualsExact(p1,p3,tol=0)
gEqualsExact(p1,p3,tol=0.1)

# pt1 and p2t are spatially identical but not exactly identical due to point ordering
pt1 = readWKT("MULTIPOINT(1 1,2 2,3 3)")
pt2 = readWKT("MULTIPOINT(3 3,2 2,1 1)")
pt3 = readWKT("MULTIPOINT(1.01 1.01,2.01 2.01,3.01 3.01)")

gEquals(pt1,pt2)
gEquals(pt1,pt3)
gEqualsExact(pt1,pt2)
gEqualsExact(pt1,pt3,tol=0)
gEqualsExact(pt1,pt3,tol=0.1)


# l2 contains a point that l1 does not
l1 = readWKT("LINESTRING (10 10, 20 20)")
l2 = readWKT("LINESTRING (10 10, 15 15,20 20)")
gEquals(l1,l2)
gEqualsExact(l1,l2)

---

Interpolate Points along Line Geometry

Description

Return points at specified distances along a line.

Usage

gInterpolate(spgeom, d, normalized = FALSE)

Arguments

spgeom	SpatialLines or SpatialLinesDataFrame object
d	Numeric vector specifying the distance along the line geometry
normalized	Logical determining if normalized distances should be used

Details

If normalized=TRUE, the distances will be interpreted as fractions of the line length.

Value

SpatialPoints object

Author(s)

Rainer Stuetz

See Also

gInterpolate

Examples

gInterpolate(readWKT("LINESTRING(25 50, 100 125, 150 190)"),
             d=seq(0, 1, by = 0.2), normalized = TRUE)

---

Geometry Relationships - Intersects and Disjoint

Description

Function for testing if the geometries have at least one point in common or no points in common

Usage

gIntersects(spgeom1, spgeom2 = NULL, byid = FALSE, prepared=TRUE,
 returnDense=TRUE, checkValidity=FALSE)
gDisjoint(spgeom1, spgeom2 = NULL, byid = FALSE, returnDense=TRUE,
 checkValidity=FALSE)

Arguments

spgeom1, spgeom2	sp objects as defined in package sp. If spgeom2 is NULL then spgeom1 is compared to itself.
byid	Logical vector determining if the function should be applied across ids (TRUE) or the entire object (FALSE) for spgeom1 and spgeom2
prepared	Logical determining if prepared geometry (spatially indexed) version of the GEOS function should be used. In general prepared geometries should be faster than the alternative.
returnDense	default TRUE, if false returns a list of the length of spgeom1 of integer vectors listing the 1:length(spgeom2) indices which would be TRUE in the dense logical matrix representation; useful when the sizes of the byid=TRUE returned matrix is very large and it is sparse; essential when the returned matrix would be too large
checkValidity	default FALSE; error meesages from GEOS do not say clearly which object fails if a topology exception is encountered. If this argument is TRUE, gIsValid is run on each in turn in an environment in which object names are available. If objects are invalid, this is reported and those affected are named

Value

gIntersects returns TRUE if spgeom1 and spgeom2 have at least one point in common.

gDisjoint returns TRUE if spgeom1 and spgeom2 have no points in common.

Both return a conforming logical matrix if byid = TRUE.

Note

Error messages from GEOS, in particular topology exceptions, report 0-based object order, so geom 0 is spgeom1, and geom 1 is spgeom2.

Author(s)

Roger Bivand & Colin Rundel

See Also

gContains gContainsProperly gCovers gCoveredBy gCrosses gEquals gEqualsExact gOverlaps gRelate gTouches gWithin

Examples

p1 = readWKT("POLYGON((0 0,1 0,1 1,0 1,0 0))")
	p2 = readWKT("POLYGON((0.5 1,0 2,1 2,0.5 1))")
	p3 = readWKT("POLYGON((0.5 0.5,0 1.5,1 1.5,0.5 0.5))")

	l1 = readWKT("LINESTRING(0 3,1 1,2 2,3 0)")
	l2 = readWKT("LINESTRING(1 3.5,3 3,2 1)")
	l3 = readWKT("LINESTRING(-0.1 0,-0.1 1.1,1 1.1)")
	
	pt1 = readWKT("MULTIPOINT(1 1,3 0,2 1)")
	pt2 = readWKT("MULTIPOINT(0 3,3 0,2 1)")
	pt3 = readWKT("MULTIPOINT(-0.2 0,1 -0.2,1.2 1,0 1.2)")

	
	par(mfrow=c(3,2))
	plot(p1,col='blue',border='blue',ylim=c(0,2.5));plot(p2,col='black',add=TRUE,pch=16)
	title(paste("Intersects:",gIntersects(p1,p2),
				"\nDisjoint:",gDisjoint(p1,p2)))
	
	plot(p1,col='blue',border='blue',ylim=c(0,2.5));plot(p3,col='black',add=TRUE,pch=16)
	title(paste("Intersects:",gIntersects(p1,p3),
				"\nDisjoint:",gDisjoint(p1,p3)))
				
	plot(l1,col='blue');plot(pt1,add=TRUE,pch=16)
	title(paste("Intersects:",gIntersects(l1,pt1),
				"\nDisjoint:",gDisjoint(l1,pt1)))
	
	plot(l1,col='blue');plot(pt2,add=TRUE,pch=16)
	title(paste("Intersects:",gIntersects(l1,pt2),
				"\nDisjoint:",gDisjoint(l1,pt2)))
	
	plot(p1,col='blue',border='blue',xlim=c(-0.5,2),ylim=c(0,2.5))
	plot(l3,lwd=2,col='black',add=TRUE)
	title(paste("Intersects:",gIntersects(p1,l3),
				"\nDisjoint:",gDisjoint(p1,l3)))
	
	plot(p1,col='blue',border='blue',xlim=c(-0.5,2),ylim=c(-0.5,2))
	plot(pt3,pch=16,col='black',add=TRUE)
	title(paste("Intersects:",gIntersects(p1,pt3),
				"\nDisjoint:",gDisjoint(p1,pt3)))

# Michael Chirico bug report and fix 2019-08-16
	SP1 = SpatialPoints(
	  cbind(c(.25, .75), c(.75, .25))	)
	SP2 = SpatialPolygons(list(
	  Polygons(list(Polygon(cbind(c(0, 0, 1, 1), c(0, 1, 1, 0)))), ID = 'a'),
	  Polygons(list(Polygon(cbind(c(1, 1, 2, 2), c(1, 2, 2, 1)))), ID = 'b')
	))

	gIntersects(SP1, SP2, byid = c(TRUE, FALSE))
	gIntersects(SP1, SP2, byid = c(TRUE, TRUE))
        gIntersects(SP1, SP2, byid = c(FALSE, TRUE))
        gIntersects(SP1, SP2, byid = c(FALSE, FALSE))

---

Is Geometry Empty?

Description

Tests if the given geometry is empty

Usage

gIsEmpty(spgeom, byid = FALSE)

Arguments

spgeom	sp object as defined in package sp
byid	Logical determining if the function should be applied across subgeometries (TRUE) or the entire object (FALSE)

Details

Because no sp Spatial object may be empty, the function exists but cannot work, as readWKT is not permitted to create an empty object.

Value

Returns TRUE if the given geometry is empty, FALSE otherwise.

Author(s)

Roger Bivand & Colin Rundel

See Also

gIsRing gIsSimple gIsValid

Examples

try(gIsEmpty(readWKT("POINT EMPTY")))
gIsEmpty(readWKT("POINT(1 1)"))
try(gIsEmpty(readWKT("LINESTRING EMPTY")))
gIsEmpty(readWKT("LINESTRING(0 0,1 1)"))
try(gIsEmpty(readWKT("POLYGON EMPTY")))
gIsEmpty(readWKT("POLYGON((0 0,1 0,1 1,0 1,0 0))"))

---

Is Geometry a Ring?

Description

Tests if the given geometry is a ring

Usage

gIsRing(spgeom, byid = FALSE)

Arguments

spgeom	sp object as defined in package sp
byid	Logical determining if the function should be applied across subgeometries (TRUE) or the entire object (FALSE)

Value

Returns TRUE if the geometry is a LINEARRING.

Returns TRUE if the geometry is a LINESTRING that is both Simple (gIsSimple) and Closed (end points intersect), FALSE otherwise.

Returns FALSE if the geometry is a [MULTI]POINT, MULTILINESTRING, or [MULTI]POLYGON.

Author(s)

Roger Bivand & Colin Rundel

See Also

gIsEmpty gIsSimple gIsValid

Examples

l1 = readWKT("LINESTRING(0 0, 1 1, 1 0, 0 1, 0 0)")
l2 = readWKT("LINESTRING(0 0,1 0,1 1,0 1,0 0)")
r1 = readWKT("LINEARRING(0 0, 1 1, 1 0, 0 1, 0 0)")
r2 = readWKT("LINEARRING(0 0,1 0,1 1,0 1,0 0)")
p1 = readWKT("POLYGON((0 0, 1 1, 1 0, 0 1, 0 0))")
p2 = readWKT("POLYGON((0 0,1 0,1 1,0 1,0 0))")

par(mfrow=c(3,2))
plot(l1);title(paste("LINESTRING\nRing:",gIsRing(l1)))
plot(l2);title(paste("LINESTRING\nRing:",gIsRing(l2)))
plot(r1);title(paste("LINEARRING\nRing:",gIsRing(r1)))
plot(r2);title(paste("LINEARRING\nRing:",gIsRing(r2)))
plot(p1);title(paste("POLYGON\nRing:",gIsRing(p1)))
plot(p2);title(paste("POLYGON\nRing:",gIsRing(p2)))

---

Is Geometry Simple?

Description

Function tests if the given geometry is simple

Usage

gIsSimple(spgeom, byid = FALSE)

Arguments

spgeom	sp object as defined in package sp
byid	Logical determining if the function should be applied across subgeometries (TRUE) or the entire object (FALSE)

Details

Simplicity is used in reference to 0 and 1-dimensional geometries ([MULTI]POINT and [MULTI]LINESTRING) whereas Validity (gIsValid) is used in reference to 2-dimensional geometries ([MULTI]POLYGON).

A POINT is always simple.

A MULTIPOINT is simple if no two points are identical.

A LINESTRING is simple if it does not pass through the same point twice (self intersection) except at the end points, in which case it is a ring (gIsRing).

A MULTILINESTRING is simple if all of its subgeometries are simple and none of the subgeometries intersect except at end points.

A [MULTI]POLYGON is simple by definition.

Many of the functions in rgeos expect simple/valid geometries and may exhibit unpredictable behavior if given an invalid geometry. Checking of validity/simplicity can be computationally expensive for complex geometries and so is not done by default, any new geometries should be checked.

Value

Returns TRUE if the given geometry does not contain anomalous points, such as self intersection or self tangency.

Author(s)

Roger Bivand & Colin Rundel

References

Validity Details: http://postgis.net/docs/manual-2.0/using_postgis_dbmanagement.html#OGC_Validity

See Also

gIsEmpty gIsRing gIsValid

Examples

# MULTIPOINT examples
gIsSimple(readWKT("MULTIPOINT(1 1,2 2,3 3)"))
gIsSimple(readWKT("MULTIPOINT(1 1,2 2,1 1)"))

# LINESTRING examples
l1 = readWKT("LINESTRING(0 5,3 4,2 3,5 2)")
l2 = readWKT("LINESTRING(0 5,4 2,5 4,0 1)")
l3 = readWKT("LINESTRING(3 5,0 4,0 2,2 0,5 1,4 4,4 5,3 5)")
l4 = readWKT("LINESTRING(3 5,0 4,4 3,5 2,3 0,1 2,4 5,3 5)")

par(mfrow=c(2,2))
plot(l1);title(paste("Simple:",gIsSimple(l1)))
plot(l2);title(paste("Simple:",gIsSimple(l2)))
plot(l3);title(paste("Simple:",gIsSimple(l3)))
plot(l4);title(paste("Simple:",gIsSimple(l4)))

# MULTILINESTRING examples
ml1 = readWKT("MULTILINESTRING((0 5,1 2,5 0),(3 5,5 4,4 1))")
ml2 = readWKT("MULTILINESTRING((0 5,1 2,5 0),(0 5,5 4,4 1))")
ml3 = readWKT("MULTILINESTRING((0 5,1 2,5 0),(3 5,5 4,2 0))")

par(mfrow=c(1,3))
plot(ml1);title(paste("Simple:",gIsSimple(ml1)))
plot(ml2);title(paste("Simple:",gIsSimple(ml2)))
plot(ml3);title(paste("Simple:",gIsSimple(ml3)))

---

Is Geometry Valid?

Description

Function tests if the given geometry is valid

Usage

gIsValid(spgeom, byid = FALSE, reason=FALSE)

Arguments

spgeom	sp object as defined in package sp
byid	Logical determining if the function should be applied across subgeometries (TRUE) or the entire object (FALSE)
reason	Logical determining if the function should return a character string describing why the geometry is invalid

Details

Validity is used in reference to 2-dimensional geometries (LINEARRING and [MULTI]POLYGON) whereas Simplicity (gIsSimple) is used in reference to 0 and 1-dimensional geometries ([MULTI]POINT and [MULTI]LINESTRING).

A LINEARRING is valid if it does not intersect itself.

A POLYGON is valid if no two rings in the boundary (made up of an exterior ring and interior rings) cross. The boundary of a POLYGON may intersect at a POINT but only as a tangent (i.e. not on a line). A POLYGON may not have cut lines or spikes and the interior rings must be contained entirely within the exterior ring.

A MULTIPOLYGON is valid if and only if all of its elements are valid and the interiors of no two elements intersect. The boundaries of any two elements may touch, but only at a finite number of POINTs.

Many of the functions in rgeos expect simple/valid geometries and may exhibit unpredictable behavior if given an invalid geometry. Checking of validity/simplicity can be computationally expensive for complex geometries and so is not done by default, any new geometries should be checked.

Value

By default will return TRUE if the given geometry is well formed, FALSE otherwise. If reason is set to TRUE then a character string is returned describing the geometry, "Valid Geometry" if it is valid or details of the specific issue. Any given geometry may have multiple issues that make it invalid, gIsValid will only return the first, once it has been corrected additionally checking is necessary to confirm that additional issues do not remain.

Author(s)

Roger Bivand & Colin Rundel

References

Validity Details: http://postgis.net/docs/manual-2.0/using_postgis_dbmanagement.html#OGC_Validity

See Also

gIsEmpty gIsRing gIsSimple

Examples

library(sp)
#LINEARRING Example
l = readWKT("LINEARRING(0 0, 100 100, 100 0, 0 100, 0 0)")
plot(l);title(paste("Valid:",gIsValid(l),"\n",gIsValid(l,reason=TRUE)))

#POLYGON and MULTIPOLYGON Examples
p1 = readWKT("POLYGON ((0 60, 0 0, 60 0, 60 60, 0 60), (20 40, 20 20, 40 20, 40 40, 20 40))")
p2 = readWKT("POLYGON ((0 60, 0 0, 60 0, 60 60, 0 60), (20 40, 20 20, 60 20, 20 40))")
p3 = readWKT(paste("POLYGON ((0 120, 0 0, 140 0, 140 120, 0 120),",
 "(100 100, 100 20, 120 20, 120 100, 100 100),",
 "(20 100, 20 40, 100 40, 20 100))"))

p4 = readWKT("POLYGON ((0 40, 0 0, 40 40, 40 0, 0 40))")
p5 = readWKT("POLYGON ((-10 50, 50 50, 50 -10, -10 -10, -10 50), (0 40, 0 0, 40 40, 40 0, 0 40))")
p6 = readWKT("POLYGON ((0 60, 0 0, 60 0, 60 20, 100 20, 60 20, 60 60, 0 60))")
p7 = readWKT(paste("POLYGON ((40 300, 40 20, 280 20, 280 300, 40 300),",
 "(120 240, 80 180, 160 220, 120 240),",
 "(220 240, 160 220, 220 160, 220 240),",
 "(160 100, 80 180, 100 80, 160 100),",
 "(160 100, 220 160, 240 100, 160 100))"))
p8 = readWKT(paste("POLYGON ((40 320, 340 320, 340 20, 40 20, 40 320),",
 "(100 120, 40 20, 180 100, 100 120),",
 "(200 200, 180 100, 240 160, 200 200),",
 "(260 260, 240 160, 300 200, 260 260),",
 "(300 300, 300 200, 340 320, 300 300))"))
p9 = readWKT(paste("MULTIPOLYGON (((20 380, 420 380, 420 20, 20 20, 20 380),",
 "(220 340, 180 240, 60 200, 200 180, 340 60, 240 220, 220 340)),",
 "((60 200, 340 60, 220 340, 60 200)))"))

par(mfrow=c(3,3))
plot(p1,col='black',pbg='white');title(paste("Valid:",gIsValid(p1),"\n",gIsValid(p1,reason=TRUE)))
plot(p2,col='black',pbg='white');title(paste("Valid:",gIsValid(p2),"\n",gIsValid(p2,reason=TRUE)))
plot(p3,col='black',pbg='white');title(paste("Valid:",gIsValid(p3),"\n",gIsValid(p3,reason=TRUE)))
plot(p4,col='black',pbg='white');title(paste("Valid:",gIsValid(p4),"\n",gIsValid(p4,reason=TRUE)))
plot(p5,col='black',pbg='white');title(paste("Valid:",gIsValid(p5),"\n",gIsValid(p5,reason=TRUE)))
plot(p6,col='black',pbg='white');title(paste("Valid:",gIsValid(p6),"\n",gIsValid(p6,reason=TRUE)))
plot(p7,col='black',pbg='white');title(paste("Valid:",gIsValid(p7),"\n",gIsValid(p7,reason=TRUE)))
plot(p8,col='black',pbg='white');title(paste("Valid:",gIsValid(p8),"\n",gIsValid(p8,reason=TRUE)))
plot(p9,col='black',pbg='white')
title(paste("Valid:",gIsValid(p9),"\n",gIsValid(p9,reason=TRUE)))

---

Length of Geometry

Description

Calculates the length of the given geometry.

Usage

gLength(spgeom, byid=FALSE)

Arguments

spgeom	sp object as defined in package sp
byid	Logical determining if the function should be applied across subgeometries (TRUE) or the entire object (FALSE)

Value

Returns the length of the geometry in the units of the current projection. By definition [MULTI]POINTs have a length of 0. The length of POLYGONs is the sum of the length of their shell and their hole(s).

Author(s)

Roger Bivand & Colin Rundel

See Also

gArea

Examples

gLength(readWKT("POINT(1 1)"))

gLength(readWKT("LINESTRING(0 0,1 1,2 2)"))
gLength(readWKT("LINESTRING(0 0,1 1,2 0,3 1)"))

gLength(readWKT("POLYGON((0 0,3 0,3 3,0 3,0 0))"))
gLength(readWKT("POLYGON((0 0,3 0,3 3,0 3,0 0),(1 1,2 1,2 2,1 2,1 1))"))

---

Make invalid geometries valid

Description

Function returns a valid geometry, not available before GEOS 3.8.0; from 3.10.0 two correction strategies

Usage

gMakeValid(spgeom, byid=FALSE, id = NULL, original=NULL, keepCollapsed=NULL)

Arguments

spgeom	sp object as defined in package sp
byid	Logical determining if the function should be applied across subgeometries (TRUE) or the entire object (FALSE)
id	Character vector defining id labels for the resulting geometries, if unspecified returned geometries will be labeled based on their parent geometries' labels.
original	default NULL; if GEOS < 3.10.0, TRUE, if GEOS >= 3.10.0, either set directly or taken from environment variable GEOS_MAKE_VALID with values LINEWORK for the original algorithm or STRUCTURE for the new algorithm. If GEOS_MAKE_VALID is STRUCTURE, original is FALSE
keepCollapsed	default NULL; If GEOS >= 3.10.0 and original is FALSE, the argument is used to control whether objects with dimension reduction are discarded (keepCollapsed is FALSE) or retained (TRUE). It may be set directly, or if NULL controlled by environment variable GEOS_MAKE_VALID_KEEPCOLLAPSED with values FALSE or TRUE. If GEOS_MAKE_VALID_KEEPCOLLAPSED is TRUE, keepCollapsed is TRUE

Details

Returns a valid geometry or collection of geometries of different types. For details on the buffered geometry fixer, see links from https://github.com/r-spatial/sf/issues/1655.

Author(s)

Roger Bivand

See Also

SpatialCollections-class

Examples

# Based on test geometries from sf
run <- FALSE
if (version_GEOS0() >= "3.8.0") run <- TRUE
if (run) {
X <- readWKT("POLYGON ((0 0, 0.5 0, 0.5 0.5, 0.5 0, 1 0, 1 1, 0 1, 0 0))")
gIsValid(X)
}
if (run) {
class(X)
}
if (run) {
row.names(X)
}
if (run) {
Y <- gMakeValid(X)
}
if (run) {
gIsValid(Y)
}
if (run) {
class(Y)
}
if (run) {
plot(slot(Y, "polyobj"))
plot(slot(Y, "lineobj"), add=TRUE, col="red")
}
if (run) {
row.names(slot(Y, "polyobj"))
}
if (run) {
row.names(slot(Y, "lineobj"))
}
run <- FALSE
if (version_GEOS0() >= "3.10.0") run <- TRUE
if (run) {
JTSplot <- function(x, fill="grey90", pts="black", main="", xlim, ylim) {
  if (inherits(x, "SpatialCollections")) {
    xl <- xp <- xpl <- NULL
    if (!is.null(slot(x, "lineobj"))) xl <- slot(x, "lineobj")
    if (!is.null(slot(x, "pointobj"))) xp <- slot(x, "pointobj")
    if (!is.null(slot(x, "polyobj"))) xpl <- slot(x, "polyobj")
    if (is.null(xl)) xl <- as(xpl, "SpatialLines")
    else xl <- rbind(xl, as(xpl, "SpatialLines"))
    if (is.null(xp)) xp <- as(xl, "SpatialPoints")
    else xp <- rbind(xp, as(xl, "SpatialPoints"))
  } else {
    xl <- as(x, "SpatialLines")
    xp <- as(xl, "SpatialPoints")
  }
  plot(coordinates(xp), type="n", main=main, xlab="", ylab="", axes=FALSE, xlim=xlim, ylim=ylim)
  plot(x, col=fill, border="transparent", add=TRUE, xlim=xlim, ylim=ylim)
  plot(xl, col=pts, add=TRUE, xlim=xlim, ylim=ylim)
  plot(xp, col=pts, add=TRUE, pch=14, xlim=xlim, ylim=ylim)
}
}
if (run) {
X <- readWKT(paste0("POLYGON ((10 70, 90 70, 90 50, 30 50, 30 30, ",
 "50 30, 50 90, 70 90, 70 10, 10 10, 10 70))", sep=""))
bb <- bbox(X)
opar <- par(mfrow=c(2,2))
JTSplot(X, fill="#dfdfff", pts="#00007c", main="input", xlim=bb[1,], ylim=bb[2,])
JTSplot(gMakeValid(X, original=TRUE), fill="#ffffa4", pts="#78b400",
 main="linework", xlim=bb[1,], ylim=bb[2,])
JTSplot(gMakeValid(X, original=FALSE), fill="#ffffa4", pts="#78b400",
 main="structure", xlim=bb[1,], ylim=bb[2,])
JTSplot(gMakeValid(X, original=FALSE, keepCollapsed=TRUE), fill="#ffffa4",
 pts="#78b400", main="structure + keep collapsed", xlim=bb[1,], ylim=bb[2,])
par(opar)
}
if (run) {
X <- readWKT(paste0("POLYGON ((10 70, 90 70, 90 10, 10 10, 10 70), ",
 "(60 80, 50 30, 0 20, 10 30, 60 80), (30 40, 60 30, 40 0, 30 40))", sep=""))
bb <- rbind(c(0, 90), c(0, 80))
opar <- par(mfrow=c(2,2))
JTSplot(X, fill="#dfdfff", pts="#00007c", main="input", xlim=bb[1,], ylim=bb[2,])
JTSplot(gMakeValid(X, original=TRUE), fill="#ffffa4", pts="#78b400",
 main="linework", xlim=bb[1,], ylim=bb[2,])
JTSplot(gMakeValid(X, original=FALSE), fill="#ffffa4", pts="#78b400",
 main="structure", xlim=bb[1,], ylim=bb[2,])
JTSplot(gMakeValid(X, original=FALSE, keepCollapsed=TRUE), fill="#ffffa4",
 pts="#78b400", main="structure + keep collapsed", xlim=bb[1,], ylim=bb[2,])
par(opar)
}
if (run) {
X <- readWKT(paste0("MULTIPOLYGON (((10 90, 60 90, 60 40, 10 40, 10 90)), ",
 "((40 60, 70 60, 70 30, 40 30, 40 60)), ((20 50, 50 50, 50 20, 20 20, 20 50)))", sep=""))
bb <- bbox(X)
opar <- par(mfrow=c(2,2))
JTSplot(X, fill="#dfdfff", pts="#00007c", main="input", xlim=bb[1,], ylim=bb[2,])
JTSplot(gMakeValid(X, original=TRUE), fill="#ffffa4", pts="#78b400",
 main="linework", xlim=bb[1,], ylim=bb[2,])
JTSplot(gMakeValid(X, original=FALSE), fill="#ffffa4", pts="#78b400",
 main="structure", xlim=bb[1,], ylim=bb[2,])
JTSplot(gMakeValid(X, original=FALSE, keepCollapsed=TRUE), fill="#ffffa4",
 pts="#78b400", main="structure + keep collapsed", xlim=bb[1,], ylim=bb[2,])
par(opar)
}
if (run) {
X <- readWKT(paste0("MULTIPOLYGON (((10 40, 40 40, 40 10, 10 10, 10 40), ", 
"(16 17, 31 32, 24 25, 16 17)), ((50 40, 50 40, 50 40, 50 40, 50 40)))", sep=""))
bb <- bbox(X)
opar <- par(mfrow=c(2,2))
JTSplot(X, fill="#dfdfff", pts="#00007c", main="input", xlim=bb[1,], ylim=bb[2,])
JTSplot(gMakeValid(X, original=TRUE), fill="#ffffa4", pts="#78b400",
 main="linework", xlim=bb[1,], ylim=bb[2,])
JTSplot(gMakeValid(X, original=FALSE), fill="#ffffa4", pts="#78b400",
 main="structure", xlim=bb[1,], ylim=bb[2,])
JTSplot(gMakeValid(X, original=FALSE, keepCollapsed=TRUE), fill="#ffffa4",
 pts="#78b400", main="structure + keep collapsed", xlim=bb[1,], ylim=bb[2,])
par(opar)
}
if (run) {
X <- readWKT(paste0("MULTIPOLYGON (((10 40, 40 40, 40 10, 10 10, 10 40), ", 
"(16 17, 31 32, 24 25, 16 17)), ((50 40, 50 40, 50 40, 50 40, 50 40)))", sep=""))
bb <- bbox(X)
opar <- par(mfrow=c(2,2))
JTSplot(X, fill="#dfdfff", pts="#00007c", main="input", xlim=bb[1,], ylim=bb[2,])
Sys.setenv("GEOS_MAKE_VALID"="LINEWORK")
JTSplot(gMakeValid(X), fill="#ffffa4", pts="#78b400", main="linework (env-var)",
 xlim=bb[1,], ylim=bb[2,])
Sys.setenv("GEOS_MAKE_VALID"="STRUCTURE")
Sys.setenv("GEOS_MAKE_VALID_KEEPCOLLAPSED"="FALSE")
JTSplot(gMakeValid(X), fill="#ffffa4", pts="#78b400", main="structure (env-vars)",
 xlim=bb[1,], ylim=bb[2,])
Sys.setenv("GEOS_MAKE_VALID_KEEPCOLLAPSED"="TRUE")
JTSplot(gMakeValid(X), fill="#ffffa4", pts="#78b400", main="structure + keep collapsed (env-vars)",
 xlim=bb[1,], ylim=bb[2,])
par(opar)
Sys.unsetenv("GEOS_MAKE_VALID")
Sys.unsetenv("GEOS_MAKE_VALID_KEEPCOLLAPSED")

}

---

Centroid of Geometry

Description

.

Usage

gMaximumInscribedCircle(spgeom, byid=FALSE, id = NULL, tol=.Machine$double.eps^(1/2),
 checkValidity=NULL)

Arguments

spgeom	sp object as defined in package sp
byid	Logical determining if the function should be applied across subgeometries (TRUE) or the entire object (FALSE)
id	Character vector defining id labels for the resulting geometries, if unspecified returned geometries will be labeled based on their parent geometries' labels.
tol	Numerical tolerance value
checkValidity	default FALSE; error meesages from GEOS do not say clearly which object fails if a topology exception is encountered. If this argument is TRUE, gIsValid is run on each in turn in an environment in which object names are available. If objects are invalid, this is reported and those affected are named

Details

.

Author(s)

Roger Bivand & Colin Rundel

See Also

gBoundary gConvexHull gEnvelope gPointOnSurface

Examples

if (version_GEOS0() >= "3.9.0") {
x = readWKT(paste("GEOMETRYCOLLECTION(POLYGON((0 0,10 0,10 10,0 10,0 0)),",
 "POLYGON((15 0,25 15,35 0,15 0)))"))

# Maximum inscribed circles of both the square and circle independently
c1 = gMaximumInscribedCircle(x, byid=TRUE) 
c1_sp <- as(c1, "SpatialPoints") # coercion of straight line segments to points
c1_sp1 <- NULL
if ((length(c1_sp) %% 2) == 0) c1_sp1 <- c1_sp[seq(1, length(c1_sp), 2)]
if (!is.null(c1_sp1)) c1_circ <- gBuffer(c1_sp1, byid=TRUE,
 width=gLength(c1, byid=TRUE))
# Maximum inscribed circle of square and circle together, needs gUnaryUnion(), 
# inscribes circle in the component permitting the largest circle
c2 = gMaximumInscribedCircle(gUnaryUnion(x))

opar <- par(mfrow=c(2,1))
plot(x)
plot(c1, col='red', add=TRUE, lwd=2)
if (!is.null(c1_sp1)) plot(c1_circ, border="red", add=TRUE)
plot(x)
plot(c2, col='blue', add=TRUE)
par(opar)
}

---

Centroid of Geometry

Description

Returns the minimum rotated rectangular POLYGON which encloses the input geometry.

Usage

gMinumumRotatedRectangle(spgeom, byid=FALSE, id = NULL)

Arguments

spgeom	sp object as defined in package sp
byid	Logical determining if the function should be applied across subgeometries (TRUE) or the entire object (FALSE)
id	Character vector defining id labels for the resulting geometries, if unspecified returned geometries will be labeled based on their parent geometries' labels.

Details

Returns the minimum rotated rectangular POLYGON which encloses the input geometry. The rectangle has width equal to the minimum diameter, and a longer length. If the convex hill of the input is degenerate (a line or point) a LINESTRING or POINT is returned. The minimum rotated rectangle can be used as an extremely generalized representation for the given geometry.

Author(s)

Roger Bivand & Colin Rundel

See Also

gBoundary gConvexHull gEnvelope gPointOnSurface

Examples

if (version_GEOS0() >= "3.9.0") {
x = readWKT(paste("GEOMETRYCOLLECTION(POLYGON((0 0,10 1,10 11,0 10,0 0)),",
 "POLYGON((15 0,25 12,35 3,15 0)))"))

# Minimum rotated rectangles of both the square and circle independently
c1 = gMinumumRotatedRectangle(x,byid=TRUE) 
# Minimum rotated rectangle of square and circle together
c2 = gMinumumRotatedRectangle(x)

opar <- par(mfrow=c(2,1))
plot(c1, border='red', lwd=2, lty=2)
plot(x,add=TRUE)
plot(c2, border='blue', lwd=2, lty=3)
plot(x,add=TRUE)
par(opar)
}

---

Closest Points of two Geometries

Description

Return closest points of two geometries.

Usage

gNearestPoints(spgeom1, spgeom2)

Arguments

spgeom1, spgeom2

sp objects as defined in package sp.

Value

The closest points of the two geometries or NULL on exception. The first point comes from spgeom1 geometry and the second point comes from spgeom2.

Author(s)

Rainer Stuetz

See Also

gDistance

Examples

g1 <- readWKT("MULTILINESTRING((34 54, 60 34), (0 10, 50 10, 100 50))")
g2 <- readWKT("MULTIPOINT(30 0, 100 30)")
plot(g1, pch=4, axes=TRUE)
plot(g2, add=TRUE)
plot(gNearestPoints(g1, g2), add=TRUE, col="red", pch=7)
gDistance(g1, g2)

---

Linestring Noder

Description

Function attempts to node a linestring object, inserting coordinates at intersection points; only available for GEOS >= 3.4.0.

Usage

gNode(spgeom);

Arguments

spgeom

an sp object inheriting from SpatialLines

Details

Because gPolygonize expects linestrings to be fully noded, as such they must not cross and must touch only at endpoints. gNodee takes an object inheriting from SpatialLines and attempts to add omitted nodes. Issue reported by Nicola Farina 21 March 2014.

Value

Returns a noded linestring object.

Author(s)

Roger Bivand

See Also

gPolygonize

Examples

library(sp)
pol1 <- readWKT(paste("POLYGON((39.936 43.446, 39.94 43.446, 39.94 43.45,",
 "39.936 43.45, 39.936 43.446))"))
pol2 <- readWKT(paste("POLYGON((39.9417 43.45, 39.9395 43.4505,",
 "39.9385 43.4462, 39.9343 43.4452, 39.9331 43.4469, 39.9417 43.45))"))
plot(pol2, axes=TRUE)
plot(pol1, add=TRUE, border="blue")
gIsValid(pol1)
gIsValid(pol2)
try(res <- gUnion(pol1, pol2))
if (version_GEOS0() > "3.4.0") {
pol2a <- gPolygonize(gNode(as(pol2, "SpatialLines")))
try(res <- gUnion(pol1, pol2a))
plot(res, add=TRUE, border="red", lty=2, lwd=2)
set.seed(1)
# rw from Jim Holtman's R-help posting 2010-12-2
n <- 1000
rw <- matrix(0, ncol = 2, nrow = n)
indx <- cbind(seq(n), sample(c(1, 2), n, TRUE))
rw[indx] <- sample(c(-1, 1), n, TRUE)
rw[,1] <- cumsum(rw[, 1])
rw[, 2] <- cumsum(rw[, 2])
slrw <- SpatialLines(list(Lines(list(Line(rw)), "1")))
res0 <- gNode(slrw)
print(length(slrw))
print(length(res0))
res <- gPolygonize(res0)
print(summary(res))
print(length(res))
plot(res0, axes=TRUE)
plot(res, add=TRUE, col=sample(rainbow(length(res))))
# library(spatstat)
# set.seed(0)
# X <- psp(runif(100), runif(100), runif(100), runif(100), window=owin())
# library(maptools)
# sppsp <- as(X, "SpatialLines")
# writeLines(writeWKT(sppsp, byid=FALSE), con="sppsp.wkt")
sppsp <- readWKT(readLines(system.file("wkts/sppsp.wkt", package="rgeos")))
plot(sppsp, axes=TRUE)
res0 <- gNode(sppsp)
res <- gPolygonize(res0)
plot(res, add=TRUE, col=sample(rainbow(length(res))))
}

---

Class "gpc.poly"

Description

A class for representing polygons composed of multiple contours, some of which may be holes.

Objects from the Class

Objects can be created by calls of the form new("gpc.poly", ...) or by reading in from a file using read.polyfile.

Slots

pts

Object of class “list”. Actually, pts is a list of lists with length equal to the number of contours in the "gpc.poly" object. Each element of pts is a list of length 3 with names x, y, and hole. x and y are vectors containing the x and y coordinates, respectively, while hole is a logical indicating whether or not the contour is a hole.

Methods

[

signature(x = "gpc.poly"): ...

append.poly

signature(x = "gpc.poly", y = "gpc.poly"): ...

area.poly

signature(object = "gpc.poly"): ...

coerce

signature(from = "matrix", to = "gpc.poly"): ...

coerce

signature(from = "data.frame", to = "gpc.poly"): ...

coerce

signature(from = "numeric", to = "gpc.poly"): ...

coerce

signature(from = "list", to = "gpc.poly"): ...

coerce

signature(from = "SpatialPolygons", to = "gpc.poly"): ...

coerce

signature(from = "gpc.poly", to = "matrix"): ...

coerce

signature(from = "gpc.poly", to = "numeric"): ...

coerce

signature(from = "gpc.poly", to = "SpatialPolygons"): ...

get.bbox

signature(x = "gpc.poly"): ...

get.pts

signature(object = "gpc.poly"): ...

intersect

signature(x = "gpc.poly", y = "gpc.poly"): ...

plot

signature(x = "gpc.poly"): The argument poly.args can be used to pass a list of additional arguments to be passed to the underlying polygon call.

scale.poly

signature(x = "gpc.poly"): ...

setdiff

signature(x = "gpc.poly", y = "gpc.poly"): ...

show

signature(object = "gpc.poly"): Scale x and y coordinates by amount xscale and yscale. By default xscale equals yscale.

symdiff

signature(x = "gpc.poly", y = "gpc.poly"): ...

union

signature(x = "gpc.poly", y = "gpc.poly"): ...

tristrip

signature(x = "gpc.poly"): ...

triangulate

signature(x = "gpc.poly"): ...

Note

The class "gpc.poly.nohole" is identical to "gpc.poly" except the hole flag for each contour of a "gpc.poly.nohole" object is always FALSE.

Author(s)

Roger D. Peng

Examples

## Make some random polygons
set.seed(100)
a <- cbind(rnorm(100), rnorm(100))
a <- a[chull(a), ]

## Convert `a' from matrix to "gpc.poly"
a <- as(a, "gpc.poly")

b <- cbind(rnorm(100), rnorm(100))
b <- as(b[chull(b), ], "gpc.poly")

## More complex polygons with an intersection
p1 <- read.polyfile(system.file("poly-ex-gpc/ex-poly1.txt", package = "rgeos"))
p2 <- read.polyfile(system.file("poly-ex-gpc/ex-poly2.txt", package = "rgeos"))

## Plot both polygons and highlight their intersection in red
plot(append.poly(p1, p2))
plot(intersect(p1, p2), poly.args = list(col = 2), add = TRUE)

## Highlight the difference p1 \ p2 in green
plot(setdiff(p1, p2), poly.args = list(col = 3), add = TRUE)

## Highlight the difference p2 \ p1 in blue
plot(setdiff(p2, p1), poly.args = list(col = 4), add = TRUE)

## Plot the union of the two polygons
plot(union(p1, p2))

## Take the non-intersect portions and create a new polygon
## combining the two contours
p.comb <- append.poly(setdiff(p1, p2), setdiff(p2, p1))
plot(p.comb, poly.args = list(col = 2, border = 0))

## Coerce from a matrix
x <- 
structure(c(0.0934073560027759, 0.192713393476752, 0.410062456627342, 
0.470020818875781, 0.41380985426787, 0.271408743927828, 0.100902151283831, 
0.0465648854961832, 0.63981588032221, 0.772382048331416,
0.753739930955121, 0.637744533947066, 0.455466052934407,
0.335327963176065, 0.399539700805524, 
0.600460299194476), .Dim = c(8, 2))
y <- 
structure(c(0.404441360166551, 0.338861901457321, 0.301387925052047, 
0.404441360166551, 0.531852879944483, 0.60117973629424, 0.625537820957668, 
0.179976985040276, 0.341542002301496, 0.445109321058688,
0.610817031070196, 0.596317606444189, 0.459608745684695,
0.215189873417722), .Dim = c(7, 2))

x1 <- as(x, "gpc.poly")
y1 <- as(y, "gpc.poly")

plot(append.poly(x1, y1))
plot(intersect(x1, y1), poly.args = list(col = 2), add = TRUE)

## Show the triangulation
#plot(append.poly(x1, y1))
#triangles <- triangulate(append.poly(x1,y1))
#for (i in 0:(nrow(triangles)/3 - 1)) 
#    polygon(triangles[3*i + 1:3,], col="lightblue")

---

Class "gpc.poly.nohole"

Description

A class for representing polygons with multiple contours but without holes.

Objects from the Class

Objects can be created by calls of the form ‘new("gpc.poly.nohole", ...) or by calling read.polyfile’.

Slots

pts

Object of class “list”. See the help for “gpc.poly” for details.

Extends

Class “gpc.poly”, directly.

Methods

coerce

signature(from = "numeric", to = "gpc.poly.nohole"): ...

Note

This class is identical to “"gpc.poly"” and is needed because the file formats for polygons without holes is slightly different from the file format for polygons with holes. For a “gpc.poly.nohole” object, the ⁠hole⁠ flag for each contour is always FALSE.

Also, write.polyfile will write the correct file format, depending on whether the object is of class “gpc.poly” or “gpc.poly.nohole”.

Author(s)

Roger D. Peng

See Also

gpc.poly-class

Examples

## None

---

Point on Surface of Geometry

Description

Function returns a point on the surface of the given geometry

Usage

gPointOnSurface(spgeom, byid=FALSE, id = NULL)

Arguments

spgeom	sp object as defined in package sp
byid	Logical determining if the function should be applied across subgeometries (TRUE) or the entire object (FALSE)
id	Character vector defining id labels for the resulting geometries, if unspecified returned geometries will be labeled based on their parent geometries' labels.

Details

Returns a SpatialPoints object with a point that intersects with the geometry

Author(s)

Roger Bivand & Colin Rundel

See Also

gBoundary gCentroid gConvexHull gEnvelope

Examples

# Based on test geometries from JTS
g1 = readWKT(paste("MULTIPOINT (60 300, 200 200, 240 240, 200 300, 40 140,",
 "80 240, 140 240, 100 160, 140 200, 60 200)"))
g2 = readWKT("LINESTRING (0 0, 7 14)")
g3 = readWKT("LINESTRING (0 0, 3 15, 6 2, 11 14, 16 5, 16 18, 2 22)")
g4 = readWKT(paste("MULTILINESTRING ((60 240, 140 300, 180 200, 40 140, 100 100, 120 220),",
 "(240 80, 260 160, 200 240, 180 340, 280 340, 240 180, 180 140, 40 200, 140 260))"))
g5 = readWKT("POLYGON ((0 0, 0 10, 10 10, 10 0, 0 0))")
g6 = readWKT(paste("MULTIPOLYGON (((50 260, 240 340, 260 100, 20 60, 90 140, 50 260),",
 "(200 280, 140 240, 180 160, 240 140, 200 280)),",
 "((380 280, 300 260, 340 100, 440 80, 380 280),",
 "(380 220, 340 200, 400 100, 380 220)))"))

par(mfrow=c(2,3))
plot(g1); plot(gPointOnSurface(g1),col='red',add=TRUE)
plot(g2); plot(gPointOnSurface(g2),col='red',add=TRUE)
plot(g3); plot(gPointOnSurface(g3),col='red',add=TRUE)
plot(g4); plot(gPointOnSurface(g4),col='red',add=TRUE)
plot(g5); plot(gPointOnSurface(g5),col='red',add=TRUE)
plot(g6); plot(gPointOnSurface(g6),col='red',add=TRUE)

---

Linestring Polygonizer

Description

Function attempts to polygonize the given list of linestrings. If the linestrings are not noded (coordinates inserted at intersection points), the function may fail. gNode may be tried to insert such missing points

Usage

gPolygonize( splist, getCutEdges=FALSE);

Arguments

splist	a list of sp SpatialLines objects
getCutEdges	Logical vector indicating if cut edges should be returned

Details

Polygonization is the process of forming polygons from linestrings which enclose an area. Linestrings are expected to be fully noded, as such they must not cross and must touch only at endpoints. gPolygonize takes a list of fully noded linestrings and forms all the polygons which are enclosed by the lines. Polygonization errors such as dangling lines or cut lines can be identified and reported.

Value

By default returns polygons generated by polygonizing the given linestrings. If getCutEdges is TRUE then any cut edges are returned.

Author(s)

Roger Bivand & Colin Rundel

See Also

gNode

Examples

library(sp)
linelist = list(readWKT("LINESTRING (0 0 , 10 10)"),
				readWKT("LINESTRING (185 221, 100 100)"),
				readWKT("LINESTRING (185 221, 88 275, 180 316)"),
				readWKT("LINESTRING (185 221, 292 281, 180 316)"), 
				readWKT("LINESTRING (189 98, 83 187, 185 221)"), 
				readWKT("LINESTRING (189 98, 325 168, 185 221)"))

par(mfrow=c(1,2))
plot(linelist[[1]],xlim=c(0,350),ylim=c(0,350))
title("Linestrings with nodes")
plot(as(linelist[[1]],"SpatialPoints"),pch=16,add=TRUE)

for(i in 2:length(linelist)) {
	plot(linelist[[i]],add=TRUE)
	plot(as(linelist[[i]],"SpatialPoints"),pch=16,add=TRUE)
}

plot(gPolygonize(linelist),xlim=c(0,350),ylim=c(0,350))
title("Polygonized Results")

gPolygonize(linelist,getCutEdges=TRUE) # no cut edges



linelist2 = list(readWKT("LINESTRING(1 3, 3 3, 3 1, 1 1, 1 3)"),
				 readWKT("LINESTRING(1 3, 3 3, 3 1, 1 1, 1 3)"))

gPolygonize(linelist2,getCutEdges=FALSE) # NULL
gPolygonize(linelist2,getCutEdges=TRUE) # Contains LineStrings
# bug fix 130206
LS = list(
readWKT("LINESTRING (425963 576719, 425980 576703)"),
readWKT("LINESTRING (425963 576719, 425882 577073)"),
readWKT("LINESTRING (425980 576703, 426082 577072)"),
readWKT("LINESTRING (425882 577073, 426082 577072)"),
readWKT("LINESTRING (426138 577068, 426082 577072)"),
readWKT("LINESTRING (426138 577068, 426420 577039)"),
readWKT("LINESTRING (426420 577039, 426554 576990)"),
readWKT("LINESTRING (426751 576924, 426776 576823)"),
readWKT("LINESTRING (426751 576924, 426783 576919)"),
readWKT("LINESTRING (426751 576924, 426714 576953)"),
readWKT("LINESTRING (426776 576823, 426783 576919)"),
readWKT("LINESTRING (426658 576966, 426554 576990)"),
readWKT("LINESTRING (426658 576966, 426667 577031)"),
readWKT("LINESTRING (426658 576966, 426714 576953)"),
readWKT("LINESTRING (426667 577031, 426714 576953)")
)
plot(gPolygonize(LS))

---

Project Points to Line Geometry

Description

Return distances along geometry to points nearest the specified points.

Usage

gProject(spgeom, sppoint, normalized = FALSE)

Arguments

spgeom	SpatialLines or SpatialLinesDataFrame object
sppoint	SpatialPoints or SpatialPointsDataFrame object
normalized	Logical determining if normalized distances should be used

Details

If normalized=TRUE, distances normalized to the length of the geometry are returned, i.e., values between 0 and 1.

Value

a numeric vector containing the distances along the line to points nearest to the specified points

Author(s)

Rainer Stuetz

See Also

gInterpolate

Examples

l <- readWKT("LINESTRING(0 1, 3 4, 5 6)")
p1 <- readWKT("MULTIPOINT(3 2, 3 5)")
frac <- gProject(l, p1)
p2 <- gInterpolate(l, frac)
plot(l, axes=TRUE)
plot(p1, col = "blue", add = TRUE)
plot(p2, col = "red", add = TRUE)

---

Geometry Relationships - Intersection Matrix Pattern (DE-9IM)

Description

Determines the relationships between two geometries by comparing the intersection of Interior, Boundary and Exterior of both geometries to each other. The results are summarized by the Dimensionally Extended 9-Intersection Matrix or DE-9IM.

Usage

gRelate(spgeom1, spgeom2 = NULL, pattern = NULL, byid = FALSE,
 checkValidity=FALSE)

Arguments

spgeom1, spgeom2	sp objects as defined in package sp. If spgeom2 is NULL then spgeom1 is compared to itself.
byid	Logical vector determining if the function should be applied across ids (TRUE) or the entire object (FALSE) for spgeom1 and spgeom2
pattern	Character string containing intersection matrix pattern to match against DE-9IM for given geometries. Wild card * or * characters allowed.
checkValidity	default FALSE; error meesages from GEOS do not say clearly which object fails if a topology exception is encountered. If this argument is TRUE, gIsValid is run on each in turn in an environment in which object names are available. If objects are invalid, this is reported and those affected are named

Details

Each geometry is decomposed into an interior, a boundary, and an exterior region, all the resulting geometries are then tested by intersection against one another resulting in 9 total comparisons. These comparisons are summarized by the dimensions of the intersection region, as such intersection at point(s) is labeled 0, at linestring(s) is labeled 1, at polygons(s) is labeled 2, and if they do not intersect labeled F.

If a pattern is specified then limited matching with wildcards is possible, * matches any character whereas T matches any non-F character.

See references for additional details.

Value

By default returns a 9 character string that represents the DE-9IM.

If pattern returns TRUE if the pattern matches the DE-9IM.

Note

Error messages from GEOS, in particular topology exceptions, report 0-based object order, so geom 0 is spgeom1, and geom 1 is spgeom2.

Author(s)

Roger Bivand & Colin Rundel

References

Documentation of Intersection Matrix Patterns: https://docs.geotools.org/stable/userguide/library/jts/dim9.html

See Also

gContains gContainsProperly gCovers gCoveredBy gCrosses gDisjoint gEquals gEqualsExact gIntersects gOverlaps gTouches gWithin

Examples

x = readWKT("POLYGON((1 0,0 1,1 2,2 1,1 0))")
	x.inter = x
	x.bound = gBoundary(x)
	
	y = readWKT("POLYGON((2 0,1 1,2 2,3 1,2 0))")
	y.inter = y
	y.bound = gBoundary(y)
	
	
	xy.inter = gIntersection(x,y)
	xy.inter.bound = gBoundary(xy.inter)
	
	xy.union = gUnion(x,y)
	bbox = gBuffer(gEnvelope(xy.union),width=0.5,joinStyle='mitre',mitreLimit=3)

	x.exter = gDifference(bbox,x)
	y.exter = gDifference(bbox,y)

	# geometry decomposition
	par(mfrow=c(2,3))
	plot(bbox,border='grey');plot(x,col="black",add=TRUE);title("Interior",ylab = "Polygon X")
	plot(bbox,border='grey');plot(x.bound,col="black",add=TRUE);title("Boundary")
	plot(bbox,border='grey');plot(x.exter,col="black",pbg='white',add=TRUE);title("Exterior")
	plot(bbox,border='grey');plot(y,col="black",add=TRUE);title(ylab = "Polygon Y")
	plot(bbox,border='grey');plot(y.bound,col="black",add=TRUE)
	plot(bbox,border='grey');plot(y.exter,col="black",pbg='white',add=TRUE)
	
	
	
	defaultplot = function() {
		plot(bbox,border='grey')
		plot(x,add=TRUE,col='red1',border="red3")
		plot(y,add=TRUE,col='blue1',border="blue3")		
		plot(xy.inter,add=TRUE,col='orange1',border="orange3")
	}
	
	# Dimensionally Extended 9-Intersection Matrix
	pat = gRelate(x,y)
	patchars = strsplit(pat,"")[[1]]
	
	par(mfrow=c(3,3))
	defaultplot(); plot(gIntersection(x.inter,y.inter),add=TRUE,col='black')
	title(paste("dim:",patchars[1]))
	defaultplot(); plot(gIntersection(x.bound,y.inter),add=TRUE,col='black',lwd=2)
	title(paste("dim:",patchars[2]))
	defaultplot(); plot(gIntersection(x.exter,y.inter),add=TRUE,col='black')
	title(paste("dim:",patchars[3]))
	
	defaultplot(); plot(gIntersection(x.inter,y.bound),add=TRUE,col='black',lwd=2)
	title(paste("dim:",patchars[4]))
	defaultplot(); plot(gIntersection(x.bound,y.bound),add=TRUE,col='black',pch=16)
	title(paste("dim:",patchars[5]))
	defaultplot(); plot(gIntersection(x.exter,y.bound),add=TRUE,col='black',lwd=2)
	title(paste("dim:",patchars[6]))
	
	defaultplot(); plot(gIntersection(x.inter,y.exter),add=TRUE,col='black')
	title(paste("dim:",patchars[7]))
	defaultplot(); plot(gIntersection(x.bound,y.exter),add=TRUE,col='black',lwd=2)
	title(paste("dim:",patchars[8]))
	defaultplot(); plot(gIntersection(x.exter,y.exter),add=TRUE,col='black')
	title(paste("dim:",patchars[9]))

---

Simplify Geometry

Description

Function simplifies the given geometry using the Douglas-Peuker algorithm

Usage

gSimplify(spgeom, tol, topologyPreserve=FALSE)

Arguments

spgeom	sp object as defined in package sp
tol	Numerical tolerance value to be used by the Douglas-Peuker algorithm
topologyPreserve	Logical determining if the algorithm should attempt to preserve the topology of the original geometry

Details

When applied to lines it is possible for the resulting geometry to lose simplicity (gIsSimple). If topologyPreserve is not specified it is also possible that the resulting geometries may no longer be valid (gIsValid). Remember to check the resulting geometry as many other functions rely on simplicity and or validity when performing their calculations.

Value

Returns a simplified version of the given geometry when applied to [MULTI]LINEs or [MULTI]POLYGONs.

Author(s)

Roger Bivand & Colin Rundel

References

Douglas-Peuker Algorithm: https://en.wikipedia.org/wiki/Ramer-Douglas-Peucker_algorithm

Examples

p = readWKT(paste("POLYGON((0 40,10 50,0 60,40 60,40 100,50 90,60 100,60",
 "60,100 60,90 50,100 40,60 40,60 0,50 10,40 0,40 40,0 40))"))
l = readWKT("LINESTRING(0 7,1 6,2 1,3 4,4 1,5 7,6 6,7 4,8 6,9 4)")

par(mfrow=c(2,4))
plot(p);title("Original")
plot(gSimplify(p,tol=10));title("tol: 10")
plot(gSimplify(p,tol=20));title("tol: 20")
plot(gSimplify(p,tol=25));title("tol: 25")

plot(l);title("Original")
plot(gSimplify(l,tol=3));title("tol: 3")
plot(gSimplify(l,tol=5));title("tol: 5")
plot(gSimplify(l,tol=7));title("tol: 7")
par(mfrow=c(1,1))

---

Geometry Relationships - Touches

Description

Functions for testing if the geometries have at least one boundary point in common, but no interior points

Usage

gTouches(spgeom1, spgeom2 = NULL, byid = FALSE, returnDense=TRUE,
 checkValidity=FALSE)

Arguments

spgeom1, spgeom2	sp objects as defined in package sp. If spgeom2 is NULL then spgeom1 is compared to itself.
byid	Logical vector determining if the function should be applied across ids (TRUE) or the entire object (FALSE) for spgeom1 and spgeom2
returnDense	default TRUE, if false returns a list of the length of spgeom1 of integer vectors listing the 1:length(spgeom2) indices which would be TRUE in the dense logical matrix representation; useful when the sizes of the byid=TRUE returned matrix is very large and it is sparse; essential when the returned matrix would be too large
checkValidity	default FALSE; error meesages from GEOS do not say clearly which object fails if a topology exception is encountered. If this argument is TRUE, gIsValid is run on each in turn in an environment in which object names are available. If objects are invalid, this is reported and those affected are named

Value

Returns TRUE if the intersection of the boundaries of the two geometries is not empty.

Note

Error messages from GEOS, in particular topology exceptions, report 0-based object order, so geom 0 is spgeom1, and geom 1 is spgeom2.

Author(s)

Roger Bivand & Colin Rundel

See Also

gContains gContainsProperly gCovers gCoveredBy gCrosses gDisjoint gEquals gEqualsExact gIntersects gOverlaps gRelate gWithin

Examples

p1 = readWKT("POLYGON((0 0,1 0,1 1,0 1,0 0))")
	p2 = readWKT("POLYGON((0 1,0.5 2,1 1,0 1))")
	p3 = readWKT("POLYGON((0.5 1,0 2,1 2,0.5 1))")
	p4 = readWKT("POLYGON((0.5 0.5,0 1.5,1 1.5,0.5 0.5))")
	
	l0 = readWKT("LINESTRING(0 1,0.5 2,1 1)")
	
	l1 = readWKT("LINESTRING(0 0,2 2)")
	l2 = readWKT("LINESTRING(1 1,2 0)")
	l3 = readWKT("LINESTRING(0 2,2 0)")
	
	
	par(mfrow=c(2,3))
	plot(p1,col='blue',border='blue',ylim=c(0,2.5));plot(p2,col='black',add=TRUE,pch=16)
	title(paste("Touches:",gTouches(p1,p2)))
	
	plot(p1,col='blue',border='blue',ylim=c(0,2.5));plot(p3,col='black',add=TRUE,pch=16)
	title(paste("Touches:",gTouches(p1,p3)))
	
	plot(p1,col='blue',border='blue',ylim=c(0,2.5));plot(p4,col='black',add=TRUE,pch=16)
	title(paste("Touches:",gTouches(p1,p4)))
	
	plot(p1,col='blue',border='blue',ylim=c(0,2.5));plot(l0,lwd=2,col='black',add=TRUE,pch=16)
	title(paste("Touches:",gTouches(p1,l0)))
	
	plot(l1,lwd=2,col='blue');plot(l2,lwd=2,col='black',add=TRUE,pch=16)
	title(paste("Touches:",gTouches(l1,l2)))
	
	plot(l1,lwd=2,col='blue');plot(l3,lwd=2,col='black',add=TRUE,pch=16)
	title(paste("Touches:",gTouches(l1,l3)))

---

Geometry Union

Description

Functions for joining intersecting geometries.

Usage

gUnionCascaded(spgeom, id = NULL)
gUnaryUnion(spgeom, id = NULL, checkValidity=NULL)
gLineMerge(spgeom, byid=FALSE, id = NULL)

Arguments

byid	Logical vector determining if the function should be applied across ids (TRUE) or the entire object (FALSE) for spgeom1 and spgeom2
id	Character vector defining id labels for the resulting geometries, if unspecified returned geometries will be labeled based on their parent geometries' labels; it may contain NA values for input objects not included in the union; it should define the memberships of the output Polygons objects
spgeom	sp Polygon(s) or Line(s) depending on the function used
checkValidity	default NULL, integer 0L (no action), 1L (check), 2L (check and try to buffer by zero distance to repair). If NULL, a value set to 0L for GEOS < 3.7.2, 1L for GEOS >= 3.7.2 is read from values assigned on load. Error meesages from GEOS do not say clearly which object fails if a topology exception is encountered. If this argument is > 0L, gIsValid is run on each in turn

Details

gUnionCascaded expects a single sp object of class SpatialPolygons with subgeometries which it unions together. gUnionCascaded can only dissolve MultiPolygon objects, so GeometryCollection objects to be dissolved, here a SpatialPolygons object, must be flattened a Polygons object; if GEOS version 3.3.0 is available, use gUnaryUnion.

gUnaryUnion expects a single sp object of class SpatialPolygons with subgeometries which it unions together; introduced in GEOS version 3.3.0, and handles GeometryCollection objects. If the id argument is used, it should be a character vector defining the memberships of the output Polygons objects, equal in length to the length of the polygons slot of spgeom.

gLineMerge is similar to gUnionCascaded but is written to work with lines, specifically it joins line segments with intersecting end points.

Note

Error messages from GEOS, in particular topology exceptions, report 0-based object order, so geom 0 is spgeom1, and geom 1 is spgeom2.

Author(s)

Roger Bivand & Colin Rundel

See Also

gDifference gIntersection gSymdifference

Examples

if (require(maptools)) {
nc1 <- readShapePoly(system.file("shapes/sids.shp", package="maptools")[1],
 proj4string=CRS("+proj=longlat +datum=NAD27"))
lps <- coordinates(nc1)
ID <- cut(lps[,1], quantile(lps[,1]), include.lowest=TRUE)
if (version_GEOS0() < "3.3.0") {
   reg4 <- gUnionCascaded(nc1, ID)
} else {
   reg4 <- gUnaryUnion(nc1, ID)
}
row.names(reg4)
par(mfrow=c(2,1))
plot(nc1)
plot(reg4)
par(mfrow=c(1,1))
}
gt <- GridTopology(c(0.05,0.05), c(0.1,0.1), c(2,2))
set.seed(1)
xv <- rnorm(length(coordinates(gt)[,1]))
xvs <- ifelse(xv > 0.2,1,0)
grd <- SpatialGridDataFrame(gt, data.frame(xvs))
spix <- as(grd, "SpatialPixelsDataFrame")
spol <- as(spix, "SpatialPolygonsDataFrame")
image(grd, axes=TRUE)
if (version_GEOS0() < "3.3.0") {
  spol1 <- gUnionCascaded(spol, as.character(spol$xvs))
} else {
  spol1 <- gUnaryUnion(spol, as.character(spol$xvs))
}
plot(spol1, add=TRUE)

---

Generics/Methods for polygon objects

Description

Some generic functions and methods for polygon objects

Usage

append.poly(x, y)
area.poly(object, ...)
get.pts(object)
get.bbox(x)
scale.poly(x, ...)
tristrip(x)
triangulate(x)

Arguments

x, object	A polygon object
y	A polygon object
...	Other arguments passed to methods

Details

The result of tristrip(x) is a list of two-column matrices. Each matrix is a tristrip, i.e. the rows are vertices of triangles, with each overlapping triple of rows corresponding to a separate triangle.

The result of triangulate(x) is a single two-column matrix. The rows are vertices of triangles, taken in non-overlapping triples.

Methods

append.poly

signature(x = "gpc.poly", y = "gpc.poly"): Combine all contours of two "gpc.poly" objects and return the combined polygon as a "gpc.poly" object.

area.poly

signature(object = "gpc.poly"): Compute and return the sum of the areas of all contours in a "gpc.poly" object.

scale.poly

signature(x = "gpc.poly"): Scale (divide) the x and y coordinates of a "gpc.poly" object by the amount xscale and yscale, respectively. Return a scaled "gpc.poly" object.

get.pts

signature(object = "gpc.poly"): Return the list of x and y coordinates of the vertices of a "gpc.poly" object.

get.bbox

signature(x = "gpc.poly"): Return the bounding box for a "gpc.poly" object.

tristrip

signature(x = "gpc.poly"): Return a tristrip list for a "gpc.poly" object.

triangulate

signature(x = "gpc.poly"): Return a matrix of vertices of a triangulation of a "gpc.poly" object.

Author(s)

Roger D. Peng; GPC Library by Alan Murta; tristrip additions by Duncan Murdoch

See Also

"gpc.poly" class documentation.

Examples

holepoly <- read.polyfile(system.file("poly-ex-gpc/hole-poly.txt",
 package ="rgeos"), nohole = FALSE)
area.poly(holepoly)
stopifnot(area.poly(holepoly) == 8)

---

Find spatial join or intersections

Description

Find spatial join or intersections

Usage

overGeomGeom(x, y, returnList = FALSE, fn = NULL, ..., minDimension = -1)
overGeomGeomDF(x, y, returnList = FALSE, fn = NULL, ..., minDimension = -1)

Arguments

x	see over
y	see over
returnList	see over
fn	see over
...	see over
minDimension	integer; if -1, gIntersects is used to find geometry intersections; if 0, 1 or 2, gRelate is used. if gRelate is used, intersecting geometries are ordered by the dimension of the intersection (2: area overlap; 1: line in common; 0: point in common), and minDimension determines the mininum dimension of intersection required.

Value

see over

Note

gRelate (minDimension > -1) is likely to be substantially slower than gIntersects.

Author(s)

Edzer Pebesma

Examples

p1 = readWKT("POLYGON((0 1,0.95 0.31,0.59 -0.81,-0.59 -0.81,-0.95 0.31,0 1))")
p2 = readWKT("POLYGON((2 2,-2 2,-2 -2,2 -2,2 2),(1 1,-1 1,-1 -1,1 -1,1 1))")

overGeomGeom(p1,p2)

---