Title: | Simple Features for R |
---|---|
Description: | Support for simple features, a standardized way to encode spatial vector data. Binds to 'GDAL' for reading and writing data, to 'GEOS' for geometrical operations, and to 'PROJ' for projection conversions and datum transformations. Uses by default the 's2' package for spherical geometry operations on ellipsoidal (long/lat) coordinates. |
Authors: | Edzer Pebesma [aut, cre] , Roger Bivand [ctb] , Etienne Racine [ctb], Michael Sumner [ctb], Ian Cook [ctb], Tim Keitt [ctb], Robin Lovelace [ctb], Hadley Wickham [ctb], Jeroen Ooms [ctb] , Kirill Müller [ctb], Thomas Lin Pedersen [ctb], Dan Baston [ctb], Dewey Dunnington [ctb] |
Maintainer: | Edzer Pebesma <[email protected]> |
License: | GPL-2 | MIT + file LICENSE |
Version: | 1.0-20 |
Built: | 2024-12-21 11:07:01 UTC |
Source: | https://github.com/r-spatial/sf |
sf
objectaggregate an sf
object, possibly union-ing geometries
## S3 method for class 'sf' aggregate( x, by, FUN, ..., do_union = TRUE, simplify = TRUE, join = st_intersects )
## S3 method for class 'sf' aggregate( x, by, FUN, ..., do_union = TRUE, simplify = TRUE, join = st_intersects )
x |
object of class sf |
by |
either a list of grouping vectors with length equal to |
FUN |
function passed on to aggregate, in case |
... |
arguments passed on to |
do_union |
logical; should grouped geometries be unioned using st_union? See details. |
simplify |
logical; see aggregate |
join |
logical spatial predicate function to use if |
In case do_union
is FALSE
, aggregate
will simply combine geometries using c.sfg. When polygons sharing a boundary are combined, this leads to geometries that are invalid; see https://github.com/r-spatial/sf/issues/681.
an sf
object with aggregated attributes and geometries; additional grouping variables having the names of names(ids)
or are named Group.i
for ids[[i]]
; see aggregate.
Does not work using the formula notation involving ~
defined in aggregate.
m1 = cbind(c(0, 0, 1, 0), c(0, 1, 1, 0)) m2 = cbind(c(0, 1, 1, 0), c(0, 0, 1, 0)) pol = st_sfc(st_polygon(list(m1)), st_polygon(list(m2))) set.seed(1985) d = data.frame(matrix(runif(15), ncol = 3)) p = st_as_sf(x = d, coords = 1:2) plot(pol) plot(p, add = TRUE) (p_ag1 = aggregate(p, pol, mean)) plot(p_ag1) # geometry same as pol # works when x overlaps multiple objects in 'by': p_buff = st_buffer(p, 0.2) plot(p_buff, add = TRUE) (p_ag2 = aggregate(p_buff, pol, mean)) # increased mean of second # with non-matching features m3 = cbind(c(0, 0, -0.1, 0), c(0, 0.1, 0.1, 0)) pol = st_sfc(st_polygon(list(m3)), st_polygon(list(m1)), st_polygon(list(m2))) (p_ag3 = aggregate(p, pol, mean)) plot(p_ag3) # In case we need to pass an argument to the join function: (p_ag4 = aggregate(p, pol, mean, join = function(x, y) st_is_within_distance(x, y, dist = 0.3)))
m1 = cbind(c(0, 0, 1, 0), c(0, 1, 1, 0)) m2 = cbind(c(0, 1, 1, 0), c(0, 0, 1, 0)) pol = st_sfc(st_polygon(list(m1)), st_polygon(list(m2))) set.seed(1985) d = data.frame(matrix(runif(15), ncol = 3)) p = st_as_sf(x = d, coords = 1:2) plot(pol) plot(p, add = TRUE) (p_ag1 = aggregate(p, pol, mean)) plot(p_ag1) # geometry same as pol # works when x overlaps multiple objects in 'by': p_buff = st_buffer(p, 0.2) plot(p_buff, add = TRUE) (p_ag2 = aggregate(p_buff, pol, mean)) # increased mean of second # with non-matching features m3 = cbind(c(0, 0, -0.1, 0), c(0, 0.1, 0.1, 0)) pol = st_sfc(st_polygon(list(m3)), st_polygon(list(m1)), st_polygon(list(m2))) (p_ag3 = aggregate(p, pol, mean)) plot(p_ag3) # In case we need to pass an argument to the join function: (p_ag4 = aggregate(p, pol, mean, join = function(x, y) st_is_within_distance(x, y, dist = 0.3)))
Spatial*
and Spatial*DataFrame
objectsas_Spatial()
allows to convert sf
and sfc
to Spatial*DataFrame
and
Spatial*
for sp
compatibility. You can also use as(x, "Spatial")
To transform
sp
objects to sf
and sfc
with as(x, "sf")
.
as_Spatial(from, cast = TRUE, IDs = paste0("ID", seq_along(from)))
as_Spatial(from, cast = TRUE, IDs = paste0("ID", seq_along(from)))
from |
object of class |
cast |
logical; if |
IDs |
character vector with IDs for the |
Package sp
supports three dimensions for POINT
and MULTIPOINT
(SpatialPoint*
).
Other geometries must be two-dimensional (XY
). Dimensions can be dropped using
st_zm()
with what = "M"
or what = "ZM"
.
For converting simple features (i.e., sf
objects) to their Spatial
counterpart, use as(obj, "Spatial")
geometry-only object deriving from Spatial
, of the appropriate class
nc <- st_read(system.file("shape/nc.shp", package="sf")) if (require(sp, quietly = TRUE)) { # convert to SpatialPolygonsDataFrame spdf <- as_Spatial(nc) # identical to spdf <- as(nc, "Spatial") # convert to SpatialPolygons as(st_geometry(nc), "Spatial") # back to sf as(spdf, "sf") }
nc <- st_read(system.file("shape/nc.shp", package="sf")) if (require(sp, quietly = TRUE)) { # convert to SpatialPolygonsDataFrame spdf <- as_Spatial(nc) # identical to spdf <- as(nc, "Spatial") # convert to SpatialPolygons as(st_geometry(nc), "Spatial") # back to sf as(spdf, "sf") }
Bind rows (features) of sf objects
Bind columns (variables) of sf objects
## S3 method for class 'sf' rbind(..., deparse.level = 1) ## S3 method for class 'sf' cbind(..., deparse.level = 1, sf_column_name = NULL) st_bind_cols(...)
## S3 method for class 'sf' rbind(..., deparse.level = 1) ## S3 method for class 'sf' cbind(..., deparse.level = 1, sf_column_name = NULL) st_bind_cols(...)
... |
objects to bind; note that for the rbind and cbind methods, all objects have to be of class |
deparse.level |
integer; see rbind |
sf_column_name |
character; specifies active geometry; passed on to st_sf |
both rbind
and cbind
have non-standard method dispatch (see cbind): the rbind
or cbind
method for sf
objects is only called when all arguments to be binded are of class sf
.
If you need to cbind
e.g. a data.frame
to an sf
, use data.frame directly and use st_sf on its result, or use bind_cols; see examples.
st_bind_cols
is deprecated; use cbind
instead.
cbind
called with multiple sf
objects warns about multiple geometry columns present when the geometry column to use is not specified by using argument sf_column_name
; see also st_sf.
crs = st_crs(3857) a = st_sf(a=1, geom = st_sfc(st_point(0:1)), crs = crs) b = st_sf(a=1, geom = st_sfc(st_linestring(matrix(1:4,2))), crs = crs) c = st_sf(a=4, geom = st_sfc(st_multilinestring(list(matrix(1:4,2)))), crs = crs) rbind(a,b,c) rbind(a,b) rbind(a,b) rbind(b,c) cbind(a,b,c) # warns if (require(dplyr, quietly = TRUE)) dplyr::bind_cols(a,b) c = st_sf(a=4, geomc = st_sfc(st_multilinestring(list(matrix(1:4,2)))), crs = crs) cbind(a,b,c, sf_column_name = "geomc") df = data.frame(x=3) st_sf(data.frame(c, df)) if (require(dplyr, quietly = TRUE)) dplyr::bind_cols(c, df)
crs = st_crs(3857) a = st_sf(a=1, geom = st_sfc(st_point(0:1)), crs = crs) b = st_sf(a=1, geom = st_sfc(st_linestring(matrix(1:4,2))), crs = crs) c = st_sf(a=4, geom = st_sfc(st_multilinestring(list(matrix(1:4,2)))), crs = crs) rbind(a,b,c) rbind(a,b) rbind(a,b) rbind(b,c) cbind(a,b,c) # warns if (require(dplyr, quietly = TRUE)) dplyr::bind_cols(a,b) c = st_sf(a=4, geomc = st_sfc(st_multilinestring(list(matrix(1:4,2)))), crs = crs) cbind(a,b,c, sf_column_name = "geomc") df = data.frame(x=3) st_sf(data.frame(c, df)) if (require(dplyr, quietly = TRUE)) dplyr::bind_cols(c, df)
TRUE
by defaultDrivers for which update should be TRUE
by default
db_drivers
db_drivers
An object of class character
of length 12.
Determine database type for R vector
Determine database type for R vector
## S4 method for signature 'PostgreSQLConnection,sf' dbDataType(dbObj, obj) ## S4 method for signature 'DBIObject,sf' dbDataType(dbObj, obj)
## S4 method for signature 'PostgreSQLConnection,sf' dbDataType(dbObj, obj) ## S4 method for signature 'DBIObject,sf' dbDataType(dbObj, obj)
dbObj |
DBIObject driver or connection. |
obj |
Object to convert |
sf
object to DatabaseWrite sf
object to Database
Write sf
object to Database
## S4 method for signature 'PostgreSQLConnection,character,sf' dbWriteTable( conn, name, value, ..., row.names = FALSE, overwrite = FALSE, append = FALSE, field.types = NULL, binary = TRUE ) ## S4 method for signature 'DBIObject,character,sf' dbWriteTable( conn, name, value, ..., row.names = FALSE, overwrite = FALSE, append = FALSE, field.types = NULL, binary = TRUE )
## S4 method for signature 'PostgreSQLConnection,character,sf' dbWriteTable( conn, name, value, ..., row.names = FALSE, overwrite = FALSE, append = FALSE, field.types = NULL, binary = TRUE ) ## S4 method for signature 'DBIObject,character,sf' dbWriteTable( conn, name, value, ..., row.names = FALSE, overwrite = FALSE, append = FALSE, field.types = NULL, binary = TRUE )
conn |
DBIObject |
name |
character vector of names (table names, fields, keywords). |
value |
a data.frame. |
... |
placeholder for future use. |
row.names |
Add a |
overwrite |
Will try to |
append |
Append rows to existing table; default |
field.types |
default |
binary |
Send geometries serialized as Well-Known Binary (WKB);
if |
Map extension to driver
extension_map
extension_map
An object of class list
of length 26.
add or remove overviews to/from a raster image
gdal_addo( file, overviews = c(2, 4, 8, 16), method = "NEAREST", layers = integer(0), options = character(0), config_options = character(0), clean = FALSE, read_only = FALSE )
gdal_addo( file, overviews = c(2, 4, 8, 16), method = "NEAREST", layers = integer(0), options = character(0), config_options = character(0), clean = FALSE, read_only = FALSE )
file |
character; file name |
overviews |
integer; overview levels |
method |
character; method to create overview; one of: nearest, average, rms, gauss, cubic, cubicspline, lanczos, average_mp, average_magphase, mode |
layers |
integer; layers to create overviews for (default: all) |
options |
character; dataset opening options |
config_options |
named character vector with GDAL config options, like |
clean |
logical; if |
read_only |
logical; if |
TRUE
, invisibly, on success
gdal_utils for access to other gdal utilities that have a C API
Native interface to gdal utils
gdal_utils( util = "info", source, destination, options = character(0), quiet = !(util %in% c("info", "gdalinfo", "ogrinfo", "vectorinfo", "mdiminfo")) || ("-multi" %in% options), processing = character(0), colorfilename = character(0), config_options = character(0), read_only = FALSE )
gdal_utils( util = "info", source, destination, options = character(0), quiet = !(util %in% c("info", "gdalinfo", "ogrinfo", "vectorinfo", "mdiminfo")) || ("-multi" %in% options), processing = character(0), colorfilename = character(0), config_options = character(0), read_only = FALSE )
util |
character; one of |
source |
character; name of input layer(s); for |
destination |
character; name of output layer |
options |
character; options for the utility |
quiet |
logical; if |
processing |
character; processing options for |
colorfilename |
character; name of color file for |
config_options |
named character vector with GDAL config options, like |
read_only |
logical; only for |
info
returns a character vector with the raster metadata; all other utils return (invisibly) a logical indicating success (i.e., TRUE
); in case of failure, an error is raised.
gdal_addo for adding overlays to a raster file; st_layers to query geometry type(s) and crs from layers in a (vector) data source
if (sf_extSoftVersion()["GDAL"] > "2.1.0") { # info utils can be used to list information about a raster # dataset. More info: https://gdal.org/programs/gdalinfo.html in_file <- system.file("tif/geomatrix.tif", package = "sf") gdal_utils("info", in_file, options = c("-mm", "-proj4")) # vectortranslate utils can be used to convert simple features data between # file formats. More info: https://gdal.org/programs/ogr2ogr.html in_file <- system.file("shape/storms_xyz.shp", package="sf") out_file <- paste0(tempfile(), ".gpkg") gdal_utils( util = "vectortranslate", source = in_file, destination = out_file, # output format must be specified for GDAL < 2.3 options = c("-f", "GPKG") ) # The parameters can be specified as c("name") or c("name", "value"). The # vectortranslate utils can perform also various operations during the # conversion process. For example, we can reproject the features during the # translation. gdal_utils( util = "vectortranslate", source = in_file, destination = out_file, options = c( "-f", "GPKG", # output file format for GDAL < 2.3 "-s_srs", "EPSG:4326", # input file SRS "-t_srs", "EPSG:2264", # output file SRS "-overwrite" ) ) st_read(out_file) # The parameter s_srs had to be specified because, in this case, the in_file # has no associated SRS. st_read(in_file) }
if (sf_extSoftVersion()["GDAL"] > "2.1.0") { # info utils can be used to list information about a raster # dataset. More info: https://gdal.org/programs/gdalinfo.html in_file <- system.file("tif/geomatrix.tif", package = "sf") gdal_utils("info", in_file, options = c("-mm", "-proj4")) # vectortranslate utils can be used to convert simple features data between # file formats. More info: https://gdal.org/programs/ogr2ogr.html in_file <- system.file("shape/storms_xyz.shp", package="sf") out_file <- paste0(tempfile(), ".gpkg") gdal_utils( util = "vectortranslate", source = in_file, destination = out_file, # output format must be specified for GDAL < 2.3 options = c("-f", "GPKG") ) # The parameters can be specified as c("name") or c("name", "value"). The # vectortranslate utils can perform also various operations during the # conversion process. For example, we can reproject the features during the # translation. gdal_utils( util = "vectortranslate", source = in_file, destination = out_file, options = c( "-f", "GPKG", # output file format for GDAL < 2.3 "-s_srs", "EPSG:4326", # input file SRS "-t_srs", "EPSG:2264", # output file SRS "-overwrite" ) ) st_read(out_file) # The parameter s_srs had to be specified because, in this case, the in_file # has no associated SRS. st_read(in_file) }
Perform geometric set operations with simple feature geometry collections
st_intersection(x, y, ...) ## S3 method for class 'sfc' st_intersection(x, y, ...) ## S3 method for class 'sf' st_intersection(x, y, ...) st_difference(x, y, ...) ## S3 method for class 'sfc' st_difference(x, y, ...) st_sym_difference(x, y, ...) st_snap(x, y, tolerance)
st_intersection(x, y, ...) ## S3 method for class 'sfc' st_intersection(x, y, ...) ## S3 method for class 'sf' st_intersection(x, y, ...) st_difference(x, y, ...) ## S3 method for class 'sfc' st_difference(x, y, ...) st_sym_difference(x, y, ...) st_snap(x, y, tolerance)
x |
object of class |
y |
object of class |
... |
arguments passed on to s2_options |
tolerance |
tolerance values used for |
When using GEOS and not using s2, a spatial index is built on argument x
; see https://r-spatial.org/r/2017/06/22/spatial-index.html. The reference for the STR tree algorithm is: Leutenegger, Scott T., Mario A. Lopez, and Jeffrey Edgington. "STR: A simple and efficient algorithm for R-tree packing." Data Engineering, 1997. Proceedings. 13th international conference on. IEEE, 1997. For the pdf, search Google Scholar.
When called with missing y
, the sfc
method for st_intersection
returns all non-empty intersections of the geometries of x
; an attribute idx
contains a list-column with the indexes of contributing geometries.
when called with a missing y
, the sf
method for st_intersection
returns an sf
object with attributes taken from the contributing feature with lowest index; two fields are added: n.overlaps
with the number of overlapping features in x
, and a list-column origins
with indexes of all overlapping features.
When st_difference
is called with a single argument,
overlapping areas are erased from geometries that are indexed at greater
numbers in the argument to x
; geometries that are empty
or contained fully inside geometries with higher priority are removed entirely.
The st_difference.sfc
method with a single argument returns an object with
an "idx"
attribute with the original index for returned geometries.
st_snap
snaps the vertices and segments of a geometry to another geometry's vertices. If y
contains more than one geometry, its geometries are merged into a collection before snapping to that collection.
(from the GEOS docs:) "A snap distance tolerance is used to control where snapping is performed. Snapping one geometry to another can improve robustness for overlay operations by eliminating nearly-coincident edges (which cause problems during noding and intersection calculation). Too much snapping can result in invalid topology being created, so the number and location of snapped vertices is decided using heuristics to determine when it is safe to snap. This can result in some potential snaps being omitted, however."
The intersection, difference or symmetric difference between two sets of geometries.
The returned object has the same class as that of the first argument (x
) with the non-empty geometries resulting from applying the operation to all geometry pairs in x
and y
. In case x
is of class sf
, the matching attributes of the original object(s) are added. The sfc
geometry list-column returned carries an attribute idx
, which is an n
-by-2 matrix with every row the index of the corresponding entries of x
and y
, respectively.
To find whether pairs of simple feature geometries intersect, use
the function st_intersects
instead of st_intersection
.
When using GEOS and not using s2 polygons contain their boundary. When using s2 this is determined by the model
defaults of s2_options, which can be overridden via the ... argument, e.g. model = "closed"
to force DE-9IM compliant behaviour of polygons (and reproduce GEOS results).
st_union for the union of simple features collections; intersect and setdiff for the base R set operations.
set.seed(131) library(sf) m = rbind(c(0,0), c(1,0), c(1,1), c(0,1), c(0,0)) p = st_polygon(list(m)) n = 100 l = vector("list", n) for (i in 1:n) l[[i]] = p + 10 * runif(2) s = st_sfc(l) plot(s, col = sf.colors(categorical = TRUE, alpha = .5)) title("overlapping squares") d = st_difference(s) # sequential differences: s1, s2-s1, s3-s2-s1, ... plot(d, col = sf.colors(categorical = TRUE, alpha = .5)) title("non-overlapping differences") i = st_intersection(s) # all intersections plot(i, col = sf.colors(categorical = TRUE, alpha = .5)) title("non-overlapping intersections") summary(lengths(st_overlaps(s, s))) # includes self-counts! summary(lengths(st_overlaps(d, d))) summary(lengths(st_overlaps(i, i))) sf = st_sf(s) i = st_intersection(sf) # all intersections plot(i["n.overlaps"]) summary(i$n.overlaps - lengths(i$origins)) # A helper function that erases all of y from x: st_erase = function(x, y) st_difference(x, st_union(st_combine(y))) poly = st_polygon(list(cbind(c(0, 0, 1, 1, 0), c(0, 1, 1, 0, 0)))) lines = st_multilinestring(list( cbind(c(0, 1), c(1, 1.05)), cbind(c(0, 1), c(0, -.05)), cbind(c(1, .95, 1), c(1.05, .5, -.05)) )) snapped = st_snap(poly, lines, tolerance=.1) plot(snapped, col='red') plot(poly, border='green', add=TRUE) plot(lines, lwd=2, col='blue', add=TRUE)
set.seed(131) library(sf) m = rbind(c(0,0), c(1,0), c(1,1), c(0,1), c(0,0)) p = st_polygon(list(m)) n = 100 l = vector("list", n) for (i in 1:n) l[[i]] = p + 10 * runif(2) s = st_sfc(l) plot(s, col = sf.colors(categorical = TRUE, alpha = .5)) title("overlapping squares") d = st_difference(s) # sequential differences: s1, s2-s1, s3-s2-s1, ... plot(d, col = sf.colors(categorical = TRUE, alpha = .5)) title("non-overlapping differences") i = st_intersection(s) # all intersections plot(i, col = sf.colors(categorical = TRUE, alpha = .5)) title("non-overlapping intersections") summary(lengths(st_overlaps(s, s))) # includes self-counts! summary(lengths(st_overlaps(d, d))) summary(lengths(st_overlaps(i, i))) sf = st_sf(s) i = st_intersection(sf) # all intersections plot(i["n.overlaps"]) summary(i$n.overlaps - lengths(i$origins)) # A helper function that erases all of y from x: st_erase = function(x, y) st_difference(x, st_union(st_combine(y))) poly = st_polygon(list(cbind(c(0, 0, 1, 1, 0), c(0, 1, 1, 0, 0)))) lines = st_multilinestring(list( cbind(c(0, 1), c(1, 1.05)), cbind(c(0, 1), c(0, -.05)), cbind(c(1, .95, 1), c(1.05, .5, -.05)) )) snapped = st_snap(poly, lines, tolerance=.1) plot(snapped, col='red') plot(poly, border='green', add=TRUE) plot(lines, lwd=2, col='blue', add=TRUE)
Geometric binary predicates on pairs of simple feature geometry sets
st_intersects(x, y, sparse = TRUE, ...) st_disjoint(x, y = x, sparse = TRUE, prepared = TRUE, ...) st_touches(x, y, sparse = TRUE, prepared = TRUE, ...) st_crosses(x, y, sparse = TRUE, prepared = TRUE, ...) st_within(x, y, sparse = TRUE, prepared = TRUE, ...) st_contains(x, y, sparse = TRUE, prepared = TRUE, ..., model = "open") st_contains_properly(x, y, sparse = TRUE, prepared = TRUE, ...) st_overlaps(x, y, sparse = TRUE, prepared = TRUE, ...) st_equals( x, y, sparse = TRUE, prepared = FALSE, ..., retain_unique = FALSE, remove_self = FALSE ) st_covers(x, y, sparse = TRUE, prepared = TRUE, ..., model = "closed") st_covered_by(x, y = x, sparse = TRUE, prepared = TRUE, ..., model = "closed") st_equals_exact(x, y, par, sparse = TRUE, prepared = FALSE, ...) st_is_within_distance(x, y = x, dist, sparse = TRUE, ..., remove_self = FALSE)
st_intersects(x, y, sparse = TRUE, ...) st_disjoint(x, y = x, sparse = TRUE, prepared = TRUE, ...) st_touches(x, y, sparse = TRUE, prepared = TRUE, ...) st_crosses(x, y, sparse = TRUE, prepared = TRUE, ...) st_within(x, y, sparse = TRUE, prepared = TRUE, ...) st_contains(x, y, sparse = TRUE, prepared = TRUE, ..., model = "open") st_contains_properly(x, y, sparse = TRUE, prepared = TRUE, ...) st_overlaps(x, y, sparse = TRUE, prepared = TRUE, ...) st_equals( x, y, sparse = TRUE, prepared = FALSE, ..., retain_unique = FALSE, remove_self = FALSE ) st_covers(x, y, sparse = TRUE, prepared = TRUE, ..., model = "closed") st_covered_by(x, y = x, sparse = TRUE, prepared = TRUE, ..., model = "closed") st_equals_exact(x, y, par, sparse = TRUE, prepared = FALSE, ...) st_is_within_distance(x, y = x, dist, sparse = TRUE, ..., remove_self = FALSE)
x |
object of class |
y |
object of class |
sparse |
logical; should a sparse index list be returned ( |
... |
Arguments passed on to
|
prepared |
logical; prepare geometry for |
model |
character; polygon/polyline model; one of "open", "semi-open" or "closed"; see Details. |
retain_unique |
logical; if |
remove_self |
logical; if |
par |
numeric; parameter used for "equals_exact" (margin); |
dist |
distance threshold; geometry indexes with distances smaller or equal to this value are returned; numeric value or units value having distance units. |
If prepared
is TRUE
, and x
contains POINT geometries and y
contains polygons, then the polygon geometries are prepared, rather than the points.
For most predicates, a spatial index is built on argument x
; see https://r-spatial.org/r/2017/06/22/spatial-index.html.
Specifically, st_intersects
, st_disjoint
, st_touches
st_crosses
, st_within
, st_contains
, st_contains_properly
, st_overlaps
, st_equals
, st_covers
and st_covered_by
all build spatial indexes for more efficient geometry calculations. st_relate
, st_equals_exact
, and do not; st_is_within_distance
uses a spatial index for geographic coordinates when sf_use_s2()
is true.
If y
is missing, st_predicate(x, x)
is effectively called, and a square matrix is returned with diagonal elements st_predicate(x[i], x[i])
.
Sparse geometry binary predicate (sgbp
) lists have the following attributes: region.id
with the row.names
of x
(if any, else 1:n
), ncol
with the number of features in y
, and predicate
with the name of the predicate used.
for model
, see https://github.com/r-spatial/s2/issues/32
st_contains_properly(A,B)
is true if A intersects B's interior, but not its edges or exterior; A contains A, but A does not properly contain A.
See also st_relate and https://en.wikipedia.org/wiki/DE-9IM for a more detailed description of the underlying algorithms.
st_equals_exact
returns true for two geometries of the same type and their vertices corresponding by index are equal up to a specified tolerance.
If sparse=FALSE
, st_predicate
(with predicate
e.g. "intersects") returns a dense logical matrix with element i,j
equal to TRUE
when predicate(x[i], y[j])
(e.g., when geometry of feature i and j intersect); if sparse=TRUE
, an object of class sgbp
is returned, which is a sparse list representation of the same matrix, with list element i
an integer vector with all indices j
for which predicate(x[i],y[j])
is TRUE
(and hence a zero-length integer vector if none of them is TRUE
). From the dense matrix, one can find out if one or more elements intersect by apply(mat, 1, any)
, and from the sparse list by lengths(lst) > 0
, see examples below.
For intersection on pairs of simple feature geometries, use
the function st_intersection
instead of st_intersects
.
pts = st_sfc(st_point(c(.5,.5)), st_point(c(1.5, 1.5)), st_point(c(2.5, 2.5))) pol = st_polygon(list(rbind(c(0,0), c(2,0), c(2,2), c(0,2), c(0,0)))) (lst = st_intersects(pts, pol)) (mat = st_intersects(pts, pol, sparse = FALSE)) # which points fall inside a polygon? apply(mat, 1, any) lengths(lst) > 0 # which points fall inside the first polygon? st_intersects(pol, pts)[[1]] # remove duplicate geometries: p1 = st_point(0:1) p2 = st_point(2:1) p = st_sf(a = letters[1:8], geom = st_sfc(p1, p1, p2, p1, p1, p2, p2, p1)) st_equals(p) st_equals(p, remove_self = TRUE) (u = st_equals(p, retain_unique = TRUE)) # retain the records with unique geometries: p[-unlist(u),]
pts = st_sfc(st_point(c(.5,.5)), st_point(c(1.5, 1.5)), st_point(c(2.5, 2.5))) pol = st_polygon(list(rbind(c(0,0), c(2,0), c(2,2), c(0,2), c(0,0)))) (lst = st_intersects(pts, pol)) (mat = st_intersects(pts, pol, sparse = FALSE)) # which points fall inside a polygon? apply(mat, 1, any) lengths(lst) > 0 # which points fall inside the first polygon? st_intersects(pol, pts)[[1]] # remove duplicate geometries: p1 = st_point(0:1) p2 = st_point(2:1) p = st_sf(a = letters[1:8], geom = st_sfc(p1, p1, p2, p1, p1, p2, p2, p1)) st_equals(p) st_equals(p, remove_self = TRUE) (u = st_equals(p, retain_unique = TRUE)) # retain the records with unique geometries: p[-unlist(u),]
Combine several feature geometries into one, without unioning or resolving internal boundaries
st_combine(x) st_union(x, y, ..., by_feature = FALSE, is_coverage = FALSE)
st_combine(x) st_union(x, y, ..., by_feature = FALSE, is_coverage = FALSE)
x |
object of class |
y |
object of class |
... |
ignored |
by_feature |
logical; if |
is_coverage |
logical; if |
st_combine
combines geometries without resolving borders, using c.sfg (analogous to c for ordinary vectors).
If st_union
is called with a single argument, x
, (with y
missing) and by_feature
is FALSE
all geometries are unioned together and an sfg
or single-geometry sfc
object is returned.
If by_feature
is TRUE
each feature geometry is unioned individually.
This can for instance be used to resolve internal boundaries after polygons were combined using st_combine
.
If y
is provided, all elements of x
and y
are unioned, pairwise if by_feature
is TRUE, or else as the Cartesian product of both sets.
Unioning a set of overlapping polygons has the effect of merging the areas (i.e. the same effect as iteratively unioning all individual polygons together). Unioning a set of LineStrings has the effect of fully noding and dissolving the input linework. In this context "fully noded" means that there will be a node or endpoint in the output for every endpoint or line segment crossing in the input. "Dissolved" means that any duplicate (e.g. coincident) line segments or portions of line segments will be reduced to a single line segment in the output. Unioning a set of Points has the effect of merging all identical points (producing a set with no duplicates).
st_combine
returns a single, combined geometry, with no resolved boundaries; returned geometries may well be invalid.
If y
is missing, st_union(x)
returns a single geometry with resolved boundaries, else the geometries for all unioned pairs of x[i]
and y[j]
.
st_intersection, st_difference, st_sym_difference
nc = st_read(system.file("shape/nc.shp", package="sf")) st_combine(nc) plot(st_union(nc))
nc = st_read(system.file("shape/nc.shp", package="sf")) st_combine(nc) plot(st_union(nc))
Compute Euclidean or great circle distance between pairs of geometries; compute, the area or the length of a set of geometries.
st_area(x, ...) ## S3 method for class 'sfc' st_area(x, ...) st_length(x, ...) st_perimeter(x, ...) st_distance( x, y, ..., dist_fun, by_element = FALSE, which = ifelse(isTRUE(st_is_longlat(x)), "Great Circle", "Euclidean"), par = 0, tolerance = 0 )
st_area(x, ...) ## S3 method for class 'sfc' st_area(x, ...) st_length(x, ...) st_perimeter(x, ...) st_distance( x, y, ..., dist_fun, by_element = FALSE, which = ifelse(isTRUE(st_is_longlat(x)), "Great Circle", "Euclidean"), par = 0, tolerance = 0 )
x |
object of class |
... |
passed on to s2_distance, s2_distance_matrix, or s2_perimeter |
y |
object of class |
dist_fun |
deprecated |
by_element |
logical; if |
which |
character; for Cartesian coordinates only: one of |
par |
for |
tolerance |
ignored if |
great circle distance calculations use by default spherical distances (s2_distance or s2_distance_matrix); if sf_use_s2()
is FALSE
, ellipsoidal distances are computed using st_geod_distance which uses function geod_inverse
from GeographicLib (part of PROJ); see Karney, Charles FF, 2013, Algorithms for geodesics, Journal of Geodesy 87(1), 43–55
If the coordinate reference system of x
was set, these functions return values with unit of measurement; see set_units.
st_area returns the area of each feature geometry, computed in the coordinate reference system used. In case x
has geodetic coordinates (unprojected), then if sf_use_s2()
is FALSE
st_geod_area is used for area calculation, if it is TRUE
then s2_area is used: the former assumes an ellipsoidal shape, the latter a spherical shape of the Earth. In case of projected data, areas are computed in flat space. The argument ...
can be used to specify radius
to s2_area, to modify the Earth radius.
st_length returns the length of a LINESTRING
or MULTILINESTRING
geometry, using the coordinate reference system. POINT
, MULTIPOINT
, POLYGON
or MULTIPOLYGON
geometries return zero.
If by_element
is FALSE
st_distance
returns a dense numeric matrix of dimension length(x) by length(y); otherwise it returns a numeric vector the same length as x
and y
with an error raised if the lengths of x
and y
are unequal. Distances involving empty geometries are NA
.
st_dimension, st_cast to convert geometry types
b0 = st_polygon(list(rbind(c(-1,-1), c(1,-1), c(1,1), c(-1,1), c(-1,-1)))) b1 = b0 + 2 b2 = b0 + c(-0.2, 2) x = st_sfc(b0, b1, b2) st_area(x) line = st_sfc(st_linestring(rbind(c(30,30), c(40,40))), crs = 4326) st_length(line) outer = matrix(c(0,0,10,0,10,10,0,10,0,0),ncol=2, byrow=TRUE) hole1 = matrix(c(1,1,1,2,2,2,2,1,1,1),ncol=2, byrow=TRUE) hole2 = matrix(c(5,5,5,6,6,6,6,5,5,5),ncol=2, byrow=TRUE) poly = st_polygon(list(outer, hole1, hole2)) mpoly = st_multipolygon(list( list(outer, hole1, hole2), list(outer + 12, hole1 + 12) )) st_length(st_sfc(poly, mpoly)) st_perimeter(poly) st_perimeter(mpoly) p = st_sfc(st_point(c(0,0)), st_point(c(0,1)), st_point(c(0,2))) st_distance(p, p) st_distance(p, p, by_element = TRUE)
b0 = st_polygon(list(rbind(c(-1,-1), c(1,-1), c(1,1), c(-1,1), c(-1,-1)))) b1 = b0 + 2 b2 = b0 + c(-0.2, 2) x = st_sfc(b0, b1, b2) st_area(x) line = st_sfc(st_linestring(rbind(c(30,30), c(40,40))), crs = 4326) st_length(line) outer = matrix(c(0,0,10,0,10,10,0,10,0,0),ncol=2, byrow=TRUE) hole1 = matrix(c(1,1,1,2,2,2,2,1,1,1),ncol=2, byrow=TRUE) hole2 = matrix(c(5,5,5,6,6,6,6,5,5,5),ncol=2, byrow=TRUE) poly = st_polygon(list(outer, hole1, hole2)) mpoly = st_multipolygon(list( list(outer, hole1, hole2), list(outer + 12, hole1 + 12) )) st_length(st_sfc(poly, mpoly)) st_perimeter(poly) st_perimeter(mpoly) p = st_sfc(st_point(c(0,0)), st_point(c(0,1)), st_point(c(0,2))) st_distance(p, p) st_distance(p, p, by_element = TRUE)
Dimension, simplicity, validity or is_empty queries on simple feature geometries
st_dimension(x, NA_if_empty = TRUE) st_is_simple(x) st_is_empty(x)
st_dimension(x, NA_if_empty = TRUE) st_is_simple(x) st_is_empty(x)
x |
object of class |
NA_if_empty |
logical; if TRUE, return NA for empty geometries |
st_dimension returns a numeric vector with 0 for points, 1 for lines, 2 for surfaces, and, if NA_if_empty
is TRUE
, NA
for empty geometries.
st_is_simple returns a logical vector, indicating for each geometry whether it is simple (e.g., not self-intersecting)
st_is_empty returns for each geometry whether it is empty
x = st_sfc( st_point(0:1), st_linestring(rbind(c(0,0),c(1,1))), st_polygon(list(rbind(c(0,0),c(1,0),c(0,1),c(0,0)))), st_multipoint(), st_linestring(), st_geometrycollection()) st_dimension(x) st_dimension(x, FALSE) ls = st_linestring(rbind(c(0,0), c(1,1), c(1,0), c(0,1))) st_is_simple(st_sfc(ls, st_point(c(0,0)))) ls = st_linestring(rbind(c(0,0), c(1,1), c(1,0), c(0,1))) st_is_empty(st_sfc(ls, st_point(), st_linestring()))
x = st_sfc( st_point(0:1), st_linestring(rbind(c(0,0),c(1,1))), st_polygon(list(rbind(c(0,0),c(1,0),c(0,1),c(0,0)))), st_multipoint(), st_linestring(), st_geometrycollection()) st_dimension(x) st_dimension(x, FALSE) ls = st_linestring(rbind(c(0,0), c(1,1), c(1,0), c(0,1))) st_is_simple(st_sfc(ls, st_point(c(0,0)))) ls = st_linestring(rbind(c(0,0), c(1,1), c(1,0), c(0,1))) st_is_empty(st_sfc(ls, st_point(), st_linestring()))
Geometric unary operations on simple feature geometries. These are all generics, with methods for sfg
, sfc
and sf
objects, returning an object of the same class. All operations work on a per-feature basis, ignoring all other features.
st_buffer( x, dist, nQuadSegs = 30, endCapStyle = "ROUND", joinStyle = "ROUND", mitreLimit = 1, singleSide = FALSE, ... ) st_boundary(x) st_convex_hull(x) st_concave_hull(x, ratio, ..., allow_holes) st_simplify(x, preserveTopology, dTolerance = 0) st_triangulate(x, dTolerance = 0, bOnlyEdges = FALSE) st_triangulate_constrained(x) st_inscribed_circle(x, dTolerance, ...) st_minimum_rotated_rectangle(x, ...) st_minimum_bounding_circle(x, ...) st_voronoi( x, envelope, dTolerance = 0, bOnlyEdges = FALSE, point_order = FALSE ) st_polygonize(x) st_line_merge(x, ..., directed = FALSE) st_centroid(x, ..., of_largest_polygon = FALSE) st_point_on_surface(x) st_reverse(x) st_node(x) st_segmentize(x, dfMaxLength, ...) st_exterior_ring(x, ...)
st_buffer( x, dist, nQuadSegs = 30, endCapStyle = "ROUND", joinStyle = "ROUND", mitreLimit = 1, singleSide = FALSE, ... ) st_boundary(x) st_convex_hull(x) st_concave_hull(x, ratio, ..., allow_holes) st_simplify(x, preserveTopology, dTolerance = 0) st_triangulate(x, dTolerance = 0, bOnlyEdges = FALSE) st_triangulate_constrained(x) st_inscribed_circle(x, dTolerance, ...) st_minimum_rotated_rectangle(x, ...) st_minimum_bounding_circle(x, ...) st_voronoi( x, envelope, dTolerance = 0, bOnlyEdges = FALSE, point_order = FALSE ) st_polygonize(x) st_line_merge(x, ..., directed = FALSE) st_centroid(x, ..., of_largest_polygon = FALSE) st_point_on_surface(x) st_reverse(x) st_node(x) st_segmentize(x, dfMaxLength, ...) st_exterior_ring(x, ...)
x |
object of class |
dist |
numeric or object of class |
nQuadSegs |
integer; number of segments per quadrant (fourth of a circle), for all or per-feature; see details |
endCapStyle |
character; style of line ends, one of 'ROUND', 'FLAT', 'SQUARE'; see details |
joinStyle |
character; style of line joins, one of 'ROUND', 'MITRE', 'BEVEL'; see details |
mitreLimit |
numeric; limit of extension for a join if |
singleSide |
logical; if |
... |
in |
ratio |
numeric; fraction convex: 1 returns the convex hulls, 0 maximally concave hulls |
allow_holes |
logical; if |
preserveTopology |
logical; carry out topology preserving
simplification? May be specified for each, or for all feature geometries.
Note that topology is preserved only for single feature geometries, not for
sets of them. If not specified (i.e. the default), then it is internally
set equal to |
dTolerance |
numeric; tolerance parameter, specified for all or for each
feature geometry. If you run |
bOnlyEdges |
logical; if |
envelope |
object of class |
point_order |
logical; preserve point order if TRUE and GEOS version >= 3.12; overrides bOnlyEdges |
directed |
logical; if |
of_largest_polygon |
logical; for |
dfMaxLength |
maximum length of a line segment. If |
st_buffer
computes a buffer around this geometry/each geometry. Depending on the spatial
coordinate system, a different engine (GEOS or S2) can be used, which have different function
arguments. The nQuadSegs
, endCapsStyle
, joinStyle
, mitreLimit
and
singleSide
parameters only work if the GEOS engine is used (i.e. projected coordinates or
when sf_use_s2()
is set to FALSE
). See postgis.net/docs/ST_Buffer.html
for details. The max_cells
and min_level
parameters (s2::s2_buffer_cells()
) work with the S2
engine (i.e. geographic coordinates) and can be used to change the buffer shape (e.g. smoothing). If
a negative buffer returns empty polygons instead of shrinking, set sf_use_s2()
to FALSE
.
st_boundary
returns the boundary of a geometry
st_convex_hull
creates the convex hull of a set of points
st_concave_hull
creates the concave hull of a geometry
st_simplify
simplifies lines by removing vertices.
st_triangulate
triangulates set of points (not constrained). st_triangulate
requires GEOS version 3.4 or above
st_triangulate_constrained
returns the constrained delaunay triangulation of polygons; requires GEOS version 3.10 or above
st_inscribed_circle
returns the maximum inscribed circle for polygon geometries.
For st_inscribed_circle
, if nQuadSegs
is 0 a 2-point LINESTRING is returned with the
center point and a boundary point of every circle, otherwise a circle (buffer) is returned where
nQuadSegs
controls the number of points per quadrant to approximate the circle.
st_inscribed_circle
requires GEOS version 3.9 or above
st_minimum_rotated_rectangle
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.
st_minimum_bounding_circle
returns a geometry which represents the "minimum bounding circle",
the smallest circle that contains the input.
st_voronoi
creates voronoi tessellation. st_voronoi
requires GEOS version 3.5 or above
st_polygonize
creates a polygon from lines that form a closed ring. In case of st_polygonize
, x
must be an object of class LINESTRING
or MULTILINESTRING
, or an sfc
geometry list-column object containing these
st_line_merge
merges lines. In case of st_line_merge
, x
must be an object of class MULTILINESTRING
, or an sfc
geometry list-column object containing these
st_centroid
gives the centroid of a geometry
st_point_on_surface
returns a point guaranteed to be on the (multi)surface.
st_reverse
reverses the nodes in a line
st_node
adds nodes to linear geometries at intersections without a node, and only works on individual linear geometries
st_segmentize
adds points to straight lines
an object of the same class of x
, with manipulated geometry.
chull for a more efficient algorithm for calculating the convex hull
## st_buffer, style options (taken from rgeos gBuffer) l1 = st_as_sfc("LINESTRING(0 0,1 5,4 5,5 2,8 2,9 4,4 6.5)") op = par(mfrow=c(2,3)) plot(st_buffer(l1, dist = 1, endCapStyle="ROUND"), reset = FALSE, main = "endCapStyle: ROUND") plot(l1,col='blue',add=TRUE) plot(st_buffer(l1, dist = 1, endCapStyle="FLAT"), reset = FALSE, main = "endCapStyle: FLAT") plot(l1,col='blue',add=TRUE) plot(st_buffer(l1, dist = 1, endCapStyle="SQUARE"), reset = FALSE, main = "endCapStyle: SQUARE") plot(l1,col='blue',add=TRUE) plot(st_buffer(l1, dist = 1, nQuadSegs=1), reset = FALSE, main = "nQuadSegs: 1") plot(l1,col='blue',add=TRUE) plot(st_buffer(l1, dist = 1, nQuadSegs=2), reset = FALSE, main = "nQuadSegs: 2") plot(l1,col='blue',add=TRUE) plot(st_buffer(l1, dist = 1, nQuadSegs= 5), reset = FALSE, main = "nQuadSegs: 5") plot(l1,col='blue',add=TRUE) par(op) l2 = st_as_sfc("LINESTRING(0 0,1 5,3 2)") op = par(mfrow = c(2, 3)) plot(st_buffer(l2, dist = 1, joinStyle="ROUND"), reset = FALSE, main = "joinStyle: ROUND") plot(l2, col = 'blue', add = TRUE) plot(st_buffer(l2, dist = 1, joinStyle="MITRE"), reset = FALSE, main = "joinStyle: MITRE") plot(l2, col= 'blue', add = TRUE) plot(st_buffer(l2, dist = 1, joinStyle="BEVEL"), reset = FALSE, main = "joinStyle: BEVEL") plot(l2, col= 'blue', add=TRUE) plot(st_buffer(l2, dist = 1, joinStyle="MITRE" , mitreLimit=0.5), reset = FALSE, main = "mitreLimit: 0.5") plot(l2, col = 'blue', add = TRUE) plot(st_buffer(l2, dist = 1, joinStyle="MITRE",mitreLimit=1), reset = FALSE, main = "mitreLimit: 1") plot(l2, col = 'blue', add = TRUE) plot(st_buffer(l2, dist = 1, joinStyle="MITRE",mitreLimit=3), reset = FALSE, main = "mitreLimit: 3") plot(l2, col = 'blue', add = TRUE) par(op) nc = st_read(system.file("shape/nc.shp", package="sf")) nc_g = st_geometry(nc) plot(st_convex_hull(nc_g)) plot(nc_g, border = grey(.5), add = TRUE) pt = st_combine(st_sfc(st_point(c(0,80)), st_point(c(120,80)), st_point(c(240,80)))) st_convex_hull(pt) # R2 st_convex_hull(st_set_crs(pt, 'OGC:CRS84')) # S2 set.seed(131) if (compareVersion(sf_extSoftVersion()[["GEOS"]], "3.11.0") > -1) { pts = cbind(runif(100), runif(100)) m = st_multipoint(pts) co = sf:::st_concave_hull(m, 0.3) coh = sf:::st_concave_hull(m, 0.3, allow_holes = TRUE) plot(co, col = 'grey') plot(coh, add = TRUE, border = 'red') plot(m, add = TRUE) } # st_simplify examples: op = par(mfrow = c(2, 3), mar = rep(0, 4)) plot(nc_g[1]) plot(st_simplify(nc_g[1], dTolerance = 1e3)) # 1000m plot(st_simplify(nc_g[1], dTolerance = 5e3)) # 5000m nc_g_planar = st_transform(nc_g, 2264) # planar coordinates, US foot plot(nc_g_planar[1]) plot(st_simplify(nc_g_planar[1], dTolerance = 1e3)) # 1000 foot plot(st_simplify(nc_g_planar[1], dTolerance = 5e3)) # 5000 foot par(op) if (compareVersion(sf_extSoftVersion()[["GEOS"]], "3.10.0") > -1) { pts = rbind(c(0,0), c(1,0), c(1,1), c(.5,.5), c(0,1), c(0,0)) po = st_polygon(list(pts)) co = st_triangulate_constrained(po) tr = st_triangulate(po) plot(po, col = NA, border = 'grey', lwd = 15) plot(tr, border = 'green', col = NA, lwd = 5, add = TRUE) plot(co, border = 'red', col = 'NA', add = TRUE) } if (compareVersion(sf_extSoftVersion()[["GEOS"]], "3.9.0") > -1) { nc_t = st_transform(nc, 'EPSG:2264') x = st_inscribed_circle(st_geometry(nc_t)) plot(st_geometry(nc_t), asp = 1, col = grey(.9)) plot(x, add = TRUE, col = '#ff9999') } set.seed(1) x = st_multipoint(matrix(runif(10),,2)) box = st_polygon(list(rbind(c(0,0),c(1,0),c(1,1),c(0,1),c(0,0)))) if (compareVersion(sf_extSoftVersion()[["GEOS"]], "3.5.0") > -1) { v = st_sfc(st_voronoi(x, st_sfc(box))) plot(v, col = 0, border = 1, axes = TRUE) plot(box, add = TRUE, col = 0, border = 1) # a larger box is returned, as documented plot(x, add = TRUE, col = 'red', cex=2, pch=16) plot(st_intersection(st_cast(v), box)) # clip to smaller box plot(x, add = TRUE, col = 'red', cex=2, pch=16) # matching Voronoi polygons to data points: # https://github.com/r-spatial/sf/issues/1030 # generate 50 random unif points: n = 100 pts = st_as_sf(data.frame(matrix(runif(n), , 2), id = 1:(n/2)), coords = c("X1", "X2")) # compute Voronoi polygons: pols = st_collection_extract(st_voronoi(do.call(c, st_geometry(pts)))) # match them to points: pts_pol = st_intersects(pts, pols) pts$pols = pols[unlist(pts_pol)] # re-order if (isTRUE(try(compareVersion(sf_extSoftVersion()["GEOS"], "3.12.0") > -1, silent = TRUE))) { pols_po = st_collection_extract(st_voronoi(do.call(c, st_geometry(pts)), point_order = TRUE)) # GEOS >= 3.12 can preserve order of inputs pts_pol_po = st_intersects(pts, pols_po) print(all(unlist(pts_pol_po) == 1:(n/2))) } plot(pts["id"], pch = 16) # ID is color plot(st_set_geometry(pts, "pols")["id"], xlim = c(0,1), ylim = c(0,1), reset = FALSE) plot(st_geometry(pts), add = TRUE) layout(matrix(1)) # reset plot layout } mls = st_multilinestring(list(matrix(c(0,0,0,1,1,1,0,0),,2,byrow=TRUE))) st_polygonize(st_sfc(mls)) mls = st_multilinestring(list(rbind(c(0,0), c(1,1)), rbind(c(2,0), c(1,1)))) st_line_merge(st_sfc(mls)) plot(nc_g, axes = TRUE) plot(st_centroid(nc_g), add = TRUE, pch = 3, col = 'red') mp = st_combine(st_buffer(st_sfc(lapply(1:3, function(x) st_point(c(x,x)))), 0.2 * 1:3)) plot(mp) plot(st_centroid(mp), add = TRUE, col = 'red') # centroid of combined geometry plot(st_centroid(mp, of_largest_polygon = TRUE), add = TRUE, col = 'blue', pch = 3) plot(nc_g, axes = TRUE) plot(st_point_on_surface(nc_g), add = TRUE, pch = 3, col = 'red') if (compareVersion(sf_extSoftVersion()[["GEOS"]], "3.7.0") > -1) { st_reverse(st_linestring(rbind(c(1,1), c(2,2), c(3,3)))) } (l = st_linestring(rbind(c(0,0), c(1,1), c(0,1), c(1,0), c(0,0)))) st_polygonize(st_node(l)) st_node(st_multilinestring(list(rbind(c(0,0), c(1,1), c(0,1), c(1,0), c(0,0))))) sf = st_sf(a=1, geom=st_sfc(st_linestring(rbind(c(0,0),c(1,1)))), crs = 4326) if (require(lwgeom, quietly = TRUE)) { seg = st_segmentize(sf, units::set_units(100, km)) seg = st_segmentize(sf, units::set_units(0.01, rad)) nrow(seg$geom[[1]]) }
## st_buffer, style options (taken from rgeos gBuffer) l1 = st_as_sfc("LINESTRING(0 0,1 5,4 5,5 2,8 2,9 4,4 6.5)") op = par(mfrow=c(2,3)) plot(st_buffer(l1, dist = 1, endCapStyle="ROUND"), reset = FALSE, main = "endCapStyle: ROUND") plot(l1,col='blue',add=TRUE) plot(st_buffer(l1, dist = 1, endCapStyle="FLAT"), reset = FALSE, main = "endCapStyle: FLAT") plot(l1,col='blue',add=TRUE) plot(st_buffer(l1, dist = 1, endCapStyle="SQUARE"), reset = FALSE, main = "endCapStyle: SQUARE") plot(l1,col='blue',add=TRUE) plot(st_buffer(l1, dist = 1, nQuadSegs=1), reset = FALSE, main = "nQuadSegs: 1") plot(l1,col='blue',add=TRUE) plot(st_buffer(l1, dist = 1, nQuadSegs=2), reset = FALSE, main = "nQuadSegs: 2") plot(l1,col='blue',add=TRUE) plot(st_buffer(l1, dist = 1, nQuadSegs= 5), reset = FALSE, main = "nQuadSegs: 5") plot(l1,col='blue',add=TRUE) par(op) l2 = st_as_sfc("LINESTRING(0 0,1 5,3 2)") op = par(mfrow = c(2, 3)) plot(st_buffer(l2, dist = 1, joinStyle="ROUND"), reset = FALSE, main = "joinStyle: ROUND") plot(l2, col = 'blue', add = TRUE) plot(st_buffer(l2, dist = 1, joinStyle="MITRE"), reset = FALSE, main = "joinStyle: MITRE") plot(l2, col= 'blue', add = TRUE) plot(st_buffer(l2, dist = 1, joinStyle="BEVEL"), reset = FALSE, main = "joinStyle: BEVEL") plot(l2, col= 'blue', add=TRUE) plot(st_buffer(l2, dist = 1, joinStyle="MITRE" , mitreLimit=0.5), reset = FALSE, main = "mitreLimit: 0.5") plot(l2, col = 'blue', add = TRUE) plot(st_buffer(l2, dist = 1, joinStyle="MITRE",mitreLimit=1), reset = FALSE, main = "mitreLimit: 1") plot(l2, col = 'blue', add = TRUE) plot(st_buffer(l2, dist = 1, joinStyle="MITRE",mitreLimit=3), reset = FALSE, main = "mitreLimit: 3") plot(l2, col = 'blue', add = TRUE) par(op) nc = st_read(system.file("shape/nc.shp", package="sf")) nc_g = st_geometry(nc) plot(st_convex_hull(nc_g)) plot(nc_g, border = grey(.5), add = TRUE) pt = st_combine(st_sfc(st_point(c(0,80)), st_point(c(120,80)), st_point(c(240,80)))) st_convex_hull(pt) # R2 st_convex_hull(st_set_crs(pt, 'OGC:CRS84')) # S2 set.seed(131) if (compareVersion(sf_extSoftVersion()[["GEOS"]], "3.11.0") > -1) { pts = cbind(runif(100), runif(100)) m = st_multipoint(pts) co = sf:::st_concave_hull(m, 0.3) coh = sf:::st_concave_hull(m, 0.3, allow_holes = TRUE) plot(co, col = 'grey') plot(coh, add = TRUE, border = 'red') plot(m, add = TRUE) } # st_simplify examples: op = par(mfrow = c(2, 3), mar = rep(0, 4)) plot(nc_g[1]) plot(st_simplify(nc_g[1], dTolerance = 1e3)) # 1000m plot(st_simplify(nc_g[1], dTolerance = 5e3)) # 5000m nc_g_planar = st_transform(nc_g, 2264) # planar coordinates, US foot plot(nc_g_planar[1]) plot(st_simplify(nc_g_planar[1], dTolerance = 1e3)) # 1000 foot plot(st_simplify(nc_g_planar[1], dTolerance = 5e3)) # 5000 foot par(op) if (compareVersion(sf_extSoftVersion()[["GEOS"]], "3.10.0") > -1) { pts = rbind(c(0,0), c(1,0), c(1,1), c(.5,.5), c(0,1), c(0,0)) po = st_polygon(list(pts)) co = st_triangulate_constrained(po) tr = st_triangulate(po) plot(po, col = NA, border = 'grey', lwd = 15) plot(tr, border = 'green', col = NA, lwd = 5, add = TRUE) plot(co, border = 'red', col = 'NA', add = TRUE) } if (compareVersion(sf_extSoftVersion()[["GEOS"]], "3.9.0") > -1) { nc_t = st_transform(nc, 'EPSG:2264') x = st_inscribed_circle(st_geometry(nc_t)) plot(st_geometry(nc_t), asp = 1, col = grey(.9)) plot(x, add = TRUE, col = '#ff9999') } set.seed(1) x = st_multipoint(matrix(runif(10),,2)) box = st_polygon(list(rbind(c(0,0),c(1,0),c(1,1),c(0,1),c(0,0)))) if (compareVersion(sf_extSoftVersion()[["GEOS"]], "3.5.0") > -1) { v = st_sfc(st_voronoi(x, st_sfc(box))) plot(v, col = 0, border = 1, axes = TRUE) plot(box, add = TRUE, col = 0, border = 1) # a larger box is returned, as documented plot(x, add = TRUE, col = 'red', cex=2, pch=16) plot(st_intersection(st_cast(v), box)) # clip to smaller box plot(x, add = TRUE, col = 'red', cex=2, pch=16) # matching Voronoi polygons to data points: # https://github.com/r-spatial/sf/issues/1030 # generate 50 random unif points: n = 100 pts = st_as_sf(data.frame(matrix(runif(n), , 2), id = 1:(n/2)), coords = c("X1", "X2")) # compute Voronoi polygons: pols = st_collection_extract(st_voronoi(do.call(c, st_geometry(pts)))) # match them to points: pts_pol = st_intersects(pts, pols) pts$pols = pols[unlist(pts_pol)] # re-order if (isTRUE(try(compareVersion(sf_extSoftVersion()["GEOS"], "3.12.0") > -1, silent = TRUE))) { pols_po = st_collection_extract(st_voronoi(do.call(c, st_geometry(pts)), point_order = TRUE)) # GEOS >= 3.12 can preserve order of inputs pts_pol_po = st_intersects(pts, pols_po) print(all(unlist(pts_pol_po) == 1:(n/2))) } plot(pts["id"], pch = 16) # ID is color plot(st_set_geometry(pts, "pols")["id"], xlim = c(0,1), ylim = c(0,1), reset = FALSE) plot(st_geometry(pts), add = TRUE) layout(matrix(1)) # reset plot layout } mls = st_multilinestring(list(matrix(c(0,0,0,1,1,1,0,0),,2,byrow=TRUE))) st_polygonize(st_sfc(mls)) mls = st_multilinestring(list(rbind(c(0,0), c(1,1)), rbind(c(2,0), c(1,1)))) st_line_merge(st_sfc(mls)) plot(nc_g, axes = TRUE) plot(st_centroid(nc_g), add = TRUE, pch = 3, col = 'red') mp = st_combine(st_buffer(st_sfc(lapply(1:3, function(x) st_point(c(x,x)))), 0.2 * 1:3)) plot(mp) plot(st_centroid(mp), add = TRUE, col = 'red') # centroid of combined geometry plot(st_centroid(mp, of_largest_polygon = TRUE), add = TRUE, col = 'blue', pch = 3) plot(nc_g, axes = TRUE) plot(st_point_on_surface(nc_g), add = TRUE, pch = 3, col = 'red') if (compareVersion(sf_extSoftVersion()[["GEOS"]], "3.7.0") > -1) { st_reverse(st_linestring(rbind(c(1,1), c(2,2), c(3,3)))) } (l = st_linestring(rbind(c(0,0), c(1,1), c(0,1), c(1,0), c(0,0)))) st_polygonize(st_node(l)) st_node(st_multilinestring(list(rbind(c(0,0), c(1,1), c(0,1), c(1,0), c(0,0))))) sf = st_sf(a=1, geom=st_sfc(st_linestring(rbind(c(0,0),c(1,1)))), crs = 4326) if (require(lwgeom, quietly = TRUE)) { seg = st_segmentize(sf, units::set_units(100, km)) seg = st_segmentize(sf, units::set_units(0.01, rad)) nrow(seg$geom[[1]]) }
Areal-weighted interpolation of polygon data
st_interpolate_aw(x, to, extensive, ...) ## S3 method for class 'sf' st_interpolate_aw(x, to, extensive, ..., keep_NA = FALSE, na.rm = FALSE)
st_interpolate_aw(x, to, extensive, ...) ## S3 method for class 'sf' st_interpolate_aw(x, to, extensive, ..., keep_NA = FALSE, na.rm = FALSE)
x |
object of class |
to |
object of class |
extensive |
logical; if TRUE, the attribute variables are assumed to be spatially extensive (like population) and the sum is preserved, otherwise, spatially intensive (like population density) and the mean is preserved. |
... |
ignored |
keep_NA |
logical; if |
na.rm |
logical; if |
if extensive
is TRUE
and na.rm
is set to TRUE
, geometries with NA
are effectively treated as having zero attribute values.
nc = st_read(system.file("shape/nc.shp", package="sf")) g = st_make_grid(nc, n = c(10, 5)) a1 = st_interpolate_aw(nc["BIR74"], g, extensive = FALSE) sum(a1$BIR74) / sum(nc$BIR74) # not close to one: property is assumed spatially intensive a2 = st_interpolate_aw(nc["BIR74"], g, extensive = TRUE) # verify mass preservation (pycnophylactic) property: sum(a2$BIR74) / sum(nc$BIR74) a1$intensive = a1$BIR74 a1$extensive = a2$BIR74 plot(a1[c("intensive", "extensive")], key.pos = 4)
nc = st_read(system.file("shape/nc.shp", package="sf")) g = st_make_grid(nc, n = c(10, 5)) a1 = st_interpolate_aw(nc["BIR74"], g, extensive = FALSE) sum(a1$BIR74) / sum(nc$BIR74) # not close to one: property is assumed spatially intensive a2 = st_interpolate_aw(nc["BIR74"], g, extensive = TRUE) # verify mass preservation (pycnophylactic) property: sum(a2$BIR74) / sum(nc$BIR74) a1$intensive = a1$BIR74 a1$extensive = a2$BIR74 plot(a1[c("intensive", "extensive")], key.pos = 4)
Search through the driver table if driver is listed
is_driver_available(drv, drivers = st_drivers())
is_driver_available(drv, drivers = st_drivers())
drv |
character. Name of driver |
drivers |
data.frame. Table containing driver names and support. Default
is from |
Search through the driver table to match a driver name with
an action (e.g. "write"
) and check if the action is supported.
is_driver_can(drv, drivers = st_drivers(), operation = "write")
is_driver_can(drv, drivers = st_drivers(), operation = "write")
drv |
character. Name of driver |
drivers |
data.frame. Table containing driver names and support. Default
is from |
operation |
character. What action to check |
Check if the columns could be of a coercable type for sf
is_geometry_column(con, x, classes = "")
is_geometry_column(con, x, classes = "")
con |
database connection |
x |
inherits data.frame |
classes |
classes inherited |
merge method for sf and data.frame object
## S3 method for class 'sf' merge(x, y, ...)
## S3 method for class 'sf' merge(x, y, ...)
x |
object of class |
y |
object of class |
... |
arguments passed on to |
a = data.frame(a = 1:3, b = 5:7) st_geometry(a) = st_sfc(st_point(c(0,0)), st_point(c(1,1)), st_point(c(2,2))) b = data.frame(x = c("a", "b", "c"), b = c(2,5,6)) merge(a, b) merge(a, b, all = TRUE)
a = data.frame(a = 1:3, b = 5:7) st_geometry(a) = st_sfc(st_point(c(0,0)), st_point(c(1,1)), st_point(c(2,2))) b = data.frame(x = c("a", "b", "c"), b = c(2,5,6)) merge(a, b) merge(a, b, all = TRUE)
Sudden Infant Death Syndrome (SIDS) sample data for North Carolina counties,
two time periods (1974-78 and 1979-84). The details of the columns can be
found in a spdep package vignette.
Please note that, though this is basically the same as nc.sids
dataset in spData
package, nc
only contains a subset of variables. The differences are
also discussed on the vignette.
A sf
object
https://r-spatial.github.io/spdep/articles/sids.html
nc <- st_read(system.file("shape/nc.shp", package="sf"))
nc <- st_read(system.file("shape/nc.shp", package="sf"))
Arithmetic operators for simple feature geometries
## S3 method for class 'sfg' Ops(e1, e2) ## S3 method for class 'sfc' Ops(e1, e2)
## S3 method for class 'sfg' Ops(e1, e2) ## S3 method for class 'sfc' Ops(e1, e2)
e1 |
object of class |
e2 |
numeric, or object of class |
in case e2
is numeric, +, -, *, /, %% and %/% add, subtract, multiply, divide, modulo, or integer-divide by e2
. In case e2
is an n x n matrix, * matrix-multiplies and / multiplies by its inverse. If e2
is an sfg
object, |, /, & and %/% result in the geometric union, difference, intersection and symmetric difference respectively, and ==
and !=
return geometric (in)equality, using st_equals. If e2
is an sfg
or sfc
object, for operations +
and -
it has to have POINT
geometries.
If e1
is of class sfc
, and e2
is a length 2 numeric, then it is considered a two-dimensional point (and if needed repeated as such) only for operations +
and -
, in other cases the individual numbers are repeated; see commented examples.
It has been reported (https://github.com/r-spatial/sf/issues/2067) that certain ATLAS versions result in invalid polygons, where the final point in a ring is no longer equal to the first point. In that case, setting the precisions with st_set_precision may help.
object of class sfg
st_point(c(1,2,3)) + 4 st_point(c(1,2,3)) * 3 + 4 m = matrix(0, 2, 2) diag(m) = c(1, 3) # affine: st_point(c(1,2)) * m + c(2,5) # world in 0-360 range: if (require(maps, quietly = TRUE)) { w = st_as_sf(map('world', plot = FALSE, fill = TRUE)) w2 = (st_geometry(w) + c(360,90)) %% c(360) - c(0,90) w3 = st_wrap_dateline(st_set_crs(w2 - c(180,0), 4326)) + c(180,0) plot(st_set_crs(w3, 4326), axes = TRUE) } (mp <- st_point(c(1,2)) + st_point(c(3,4))) # MULTIPOINT (1 2, 3 4) mp - st_point(c(3,4)) # POINT (1 2) opar = par(mfrow = c(2,2), mar = c(0, 0, 1, 0)) a = st_buffer(st_point(c(0,0)), 2) b = a + c(2, 0) p = function(m) { plot(c(a,b)); plot(eval(parse(text=m)), col=grey(.9), add = TRUE); title(m) } o = lapply(c('a | b', 'a / b', 'a & b', 'a %/% b'), p) par(opar) sfc = st_sfc(st_point(0:1), st_point(2:3)) sfc + c(2,3) # added to EACH geometry sfc * c(2,3) # first geometry multiplied by 2, second by 3 nc = st_transform(st_read(system.file("gpkg/nc.gpkg", package="sf")), 32119) # nc state plane, m b = st_buffer(st_centroid(st_union(nc)), units::set_units(50, km)) # shoot a hole in nc: plot(st_geometry(nc) / b, col = grey(.9))
st_point(c(1,2,3)) + 4 st_point(c(1,2,3)) * 3 + 4 m = matrix(0, 2, 2) diag(m) = c(1, 3) # affine: st_point(c(1,2)) * m + c(2,5) # world in 0-360 range: if (require(maps, quietly = TRUE)) { w = st_as_sf(map('world', plot = FALSE, fill = TRUE)) w2 = (st_geometry(w) + c(360,90)) %% c(360) - c(0,90) w3 = st_wrap_dateline(st_set_crs(w2 - c(180,0), 4326)) + c(180,0) plot(st_set_crs(w3, 4326), axes = TRUE) } (mp <- st_point(c(1,2)) + st_point(c(3,4))) # MULTIPOINT (1 2, 3 4) mp - st_point(c(3,4)) # POINT (1 2) opar = par(mfrow = c(2,2), mar = c(0, 0, 1, 0)) a = st_buffer(st_point(c(0,0)), 2) b = a + c(2, 0) p = function(m) { plot(c(a,b)); plot(eval(parse(text=m)), col=grey(.9), add = TRUE); title(m) } o = lapply(c('a | b', 'a / b', 'a & b', 'a %/% b'), p) par(opar) sfc = st_sfc(st_point(0:1), st_point(2:3)) sfc + c(2,3) # added to EACH geometry sfc * c(2,3) # first geometry multiplied by 2, second by 3 nc = st_transform(st_read(system.file("gpkg/nc.gpkg", package="sf")), 32119) # nc state plane, m b = st_buffer(st_centroid(st_union(nc)), units::set_units(50, km)) # shoot a hole in nc: plot(st_geometry(nc) / b, col = grey(.9))
plot one or more attributes of an sf object on a map Plot sf object
## S3 method for class 'sf' plot( x, y, ..., main, pal = NULL, nbreaks = 10, breaks = "pretty", max.plot = getOption("sf_max.plot", default = 9), key.pos = get_key_pos(x, ...), key.length = 0.618, key.width = kw_dflt(x, key.pos), reset = TRUE, logz = FALSE, extent = x, xlim = st_bbox(extent)[c(1, 3)], ylim = st_bbox(extent)[c(2, 4)], compact = FALSE ) get_key_pos(x, ...) ## S3 method for class 'sfc_POINT' plot( x, y, ..., pch = 1, cex = 1, col = 1, bg = 0, lwd = 1, lty = 1, type = "p", add = FALSE ) ## S3 method for class 'sfc_MULTIPOINT' plot( x, y, ..., pch = 1, cex = 1, col = 1, bg = 0, lwd = 1, lty = 1, type = "p", add = FALSE ) ## S3 method for class 'sfc_LINESTRING' plot(x, y, ..., lty = 1, lwd = 1, col = 1, pch = 1, type = "l", add = FALSE) ## S3 method for class 'sfc_CIRCULARSTRING' plot(x, y, ...) ## S3 method for class 'sfc_MULTILINESTRING' plot(x, y, ..., lty = 1, lwd = 1, col = 1, pch = 1, type = "l", add = FALSE) ## S3 method for class 'sfc_POLYGON' plot( x, y, ..., lty = 1, lwd = 1, col = NA, cex = 1, pch = NA, border = 1, add = FALSE, rule = "evenodd", xpd = par("xpd") ) ## S3 method for class 'sfc_MULTIPOLYGON' plot( x, y, ..., lty = 1, lwd = 1, col = NA, border = 1, add = FALSE, rule = "evenodd", xpd = par("xpd") ) ## S3 method for class 'sfc_GEOMETRYCOLLECTION' plot( x, y, ..., pch = 1, cex = 1, bg = 0, lty = 1, lwd = 1, col = 1, border = 1, add = FALSE ) ## S3 method for class 'sfc_GEOMETRY' plot( x, y, ..., pch = 1, cex = 1, bg = 0, lty = 1, lwd = 1, col = ifelse(st_dimension(x) == 2, NA, 1), border = 1, add = FALSE ) ## S3 method for class 'sfg' plot(x, ...) plot_sf( x, xlim = NULL, ylim = NULL, asp = NA, axes = FALSE, bgc = par("bg"), ..., xaxs, yaxs, lab, setParUsrBB = FALSE, bgMap = NULL, expandBB = c(0, 0, 0, 0), graticule = NA_crs_, col_graticule = "grey", border, extent = x ) sf.colors(n = 10, cutoff.tails = c(0.35, 0.2), alpha = 1, categorical = FALSE) ## S3 method for class 'sf' text(x, labels = row.names(x), ...) ## S3 method for class 'sfc' text(x, labels = seq_along(x), ..., of_largest_polygon = FALSE) ## S3 method for class 'sf' points(x, ...) ## S3 method for class 'sfc' points(x, ..., of_largest_polygon = FALSE)
## S3 method for class 'sf' plot( x, y, ..., main, pal = NULL, nbreaks = 10, breaks = "pretty", max.plot = getOption("sf_max.plot", default = 9), key.pos = get_key_pos(x, ...), key.length = 0.618, key.width = kw_dflt(x, key.pos), reset = TRUE, logz = FALSE, extent = x, xlim = st_bbox(extent)[c(1, 3)], ylim = st_bbox(extent)[c(2, 4)], compact = FALSE ) get_key_pos(x, ...) ## S3 method for class 'sfc_POINT' plot( x, y, ..., pch = 1, cex = 1, col = 1, bg = 0, lwd = 1, lty = 1, type = "p", add = FALSE ) ## S3 method for class 'sfc_MULTIPOINT' plot( x, y, ..., pch = 1, cex = 1, col = 1, bg = 0, lwd = 1, lty = 1, type = "p", add = FALSE ) ## S3 method for class 'sfc_LINESTRING' plot(x, y, ..., lty = 1, lwd = 1, col = 1, pch = 1, type = "l", add = FALSE) ## S3 method for class 'sfc_CIRCULARSTRING' plot(x, y, ...) ## S3 method for class 'sfc_MULTILINESTRING' plot(x, y, ..., lty = 1, lwd = 1, col = 1, pch = 1, type = "l", add = FALSE) ## S3 method for class 'sfc_POLYGON' plot( x, y, ..., lty = 1, lwd = 1, col = NA, cex = 1, pch = NA, border = 1, add = FALSE, rule = "evenodd", xpd = par("xpd") ) ## S3 method for class 'sfc_MULTIPOLYGON' plot( x, y, ..., lty = 1, lwd = 1, col = NA, border = 1, add = FALSE, rule = "evenodd", xpd = par("xpd") ) ## S3 method for class 'sfc_GEOMETRYCOLLECTION' plot( x, y, ..., pch = 1, cex = 1, bg = 0, lty = 1, lwd = 1, col = 1, border = 1, add = FALSE ) ## S3 method for class 'sfc_GEOMETRY' plot( x, y, ..., pch = 1, cex = 1, bg = 0, lty = 1, lwd = 1, col = ifelse(st_dimension(x) == 2, NA, 1), border = 1, add = FALSE ) ## S3 method for class 'sfg' plot(x, ...) plot_sf( x, xlim = NULL, ylim = NULL, asp = NA, axes = FALSE, bgc = par("bg"), ..., xaxs, yaxs, lab, setParUsrBB = FALSE, bgMap = NULL, expandBB = c(0, 0, 0, 0), graticule = NA_crs_, col_graticule = "grey", border, extent = x ) sf.colors(n = 10, cutoff.tails = c(0.35, 0.2), alpha = 1, categorical = FALSE) ## S3 method for class 'sf' text(x, labels = row.names(x), ...) ## S3 method for class 'sfc' text(x, labels = seq_along(x), ..., of_largest_polygon = FALSE) ## S3 method for class 'sf' points(x, ...) ## S3 method for class 'sfc' points(x, ..., of_largest_polygon = FALSE)
x |
object of class sf |
y |
ignored |
... |
|
main |
title for plot ( |
pal |
palette function, similar to rainbow, or palette values; if omitted, |
nbreaks |
number of colors breaks (ignored for |
breaks |
either a numeric vector with the actual breaks, or a name of a method accepted by the |
max.plot |
integer; lower boundary to maximum number of attributes to plot; the default value (9) can be overridden by setting the global option |
key.pos |
numeric; side to plot a color key: 1 bottom, 2 left, 3 top, 4 right; set to |
key.length |
amount of space reserved for the key along its axis, length of the scale bar |
key.width |
amount of space reserved for the key (incl. labels), thickness/width of the scale bar |
reset |
logical; if |
logz |
logical; if |
extent |
object with an |
xlim |
see plot.window |
ylim |
see plot.window |
compact |
logical; compact sub-plots over plotting space? |
pch |
plotting symbol |
cex |
symbol size |
col |
color for plotting features; if |
bg |
symbol background color |
lwd |
line width |
lty |
line type |
type |
plot type: 'p' for points, 'l' for lines, 'b' for both |
add |
logical; add to current plot? Note that when using |
border |
color of polygon border(s); using |
rule |
see polypath; for |
xpd |
see par; sets polygon clipping strategy; only implemented for POLYGON and MULTIPOLYGON |
asp |
see below, and see par |
axes |
logical; should axes be plotted? (default FALSE) |
bgc |
background color |
xaxs |
see par |
yaxs |
see par |
lab |
see par |
setParUsrBB |
default FALSE; set the |
bgMap |
object of class |
expandBB |
numeric; fractional values to expand the bounding box with, in each direction (bottom, left, top, right) |
graticule |
logical, or object of class |
col_graticule |
color to used for the graticule (if present) |
n |
integer; number of colors |
cutoff.tails |
numeric, in |
alpha |
numeric, in |
categorical |
logical; do we want colors for a categorical variable? (see details) |
labels |
character, text to draw (one per row of input) |
of_largest_polygon |
logical, passed on to st_centroid |
plot.sf
maximally plots max.plot
maps with colors following from attribute columns,
one map per attribute. It uses sf.colors
for default colors. For more control over placement of individual maps,
set parameter mfrow
with par prior to plotting, and plot single maps one by one; note that this only works
in combination with setting parameters key.pos=NULL
(no legend) and reset=FALSE
.
plot.sfc
plots the geometry, additional parameters can be passed on
to control color, lines or symbols.
When setting reset
to FALSE
, the original device parameters are lost, and the device must be reset using dev.off()
in order to reset it.
parameter at
can be set to specify where labels are placed along the key; see examples.
The features are plotted in the order as they apppear in the sf object. See examples for when a different plotting order is wanted.
plot_sf
sets up the plotting area, axes, graticule, or webmap background; it
is called by all plot
methods before anything is drawn.
The argument setParUsrBB
may be used to pass the logical value TRUE
to functions within plot.Spatial
. When set to TRUE
, par(“usr”) will be overwritten with c(xlim, ylim)
, which defaults to the bounding box of the spatial object. This is only needed in the particular context of graphic output to a specified device with given width and height, to be matched to the spatial object, when using par(“xaxs”) and par(“yaxs”) in addition to par(mar=c(0,0,0,0))
.
The default aspect for map plots is 1; if however data are not
projected (coordinates are long/lat), the aspect is by default set to
1/cos(My * pi/180) with My the y coordinate of the middle of the map
(the mean of ylim
, which defaults to the y range of bounding box). This
implies an Equirectangular projection.
non-categorical colors from sf.colors
were taken from bpy.colors, with modified cutoff.tails
defaults
If categorical is TRUE
, default colors are from https://colorbrewer2.org/ (if n < 9, Set2, else Set3).
text.sf
adds text to an existing base graphic. Text is placed at the centroid of
each feature in x
. Provide POINT features for further control of placement.
points.sf
adds point symbols to an existing base graphic. If points of text are not shown
correctly, try setting argument reset
to FALSE
in the plot()
call.
nc = st_read(system.file("gpkg/nc.gpkg", package="sf"), quiet = TRUE) # plot single attribute, auto-legend: plot(nc["SID74"]) # plot multiple: plot(nc[c("SID74", "SID79")]) # better use ggplot2::geom_sf to facet and get a single legend! # adding to a plot of an sf object only works when using reset=FALSE in the first plot: plot(nc["SID74"], reset = FALSE) plot(st_centroid(st_geometry(nc)), add = TRUE) # log10 z-scale: plot(nc["SID74"], logz = TRUE, breaks = c(0,.5,1,1.5,2), at = c(0,.5,1,1.5,2)) # and we need to reset the plotting device after that, e.g. by layout(1) # when plotting only geometries, the reset=FALSE is not needed: plot(st_geometry(nc)) plot(st_geometry(nc)[1], col = 'red', add = TRUE) # add a custom legend to an arbitray plot: layout(matrix(1:2, ncol = 2), widths = c(1, lcm(2))) plot(1) .image_scale(1:10, col = sf.colors(9), key.length = lcm(8), key.pos = 4, at = 1:10) # manipulate plotting order, plot largest polygons first: p = st_polygon(list(rbind(c(0,0), c(1,0), c(1,1), c(0,1), c(0,0)))) x = st_sf(a=1:4, st_sfc(p, p * 2, p * 3, p * 4)) # plot(x, col=2:5) only shows the largest polygon! plot(x[order(st_area(x), decreasing = TRUE),], col = 2:5) # plot largest polygons first sf.colors(10) text(nc, labels = substring(nc$NAME,1,1))
nc = st_read(system.file("gpkg/nc.gpkg", package="sf"), quiet = TRUE) # plot single attribute, auto-legend: plot(nc["SID74"]) # plot multiple: plot(nc[c("SID74", "SID79")]) # better use ggplot2::geom_sf to facet and get a single legend! # adding to a plot of an sf object only works when using reset=FALSE in the first plot: plot(nc["SID74"], reset = FALSE) plot(st_centroid(st_geometry(nc)), add = TRUE) # log10 z-scale: plot(nc["SID74"], logz = TRUE, breaks = c(0,.5,1,1.5,2), at = c(0,.5,1,1.5,2)) # and we need to reset the plotting device after that, e.g. by layout(1) # when plotting only geometries, the reset=FALSE is not needed: plot(st_geometry(nc)) plot(st_geometry(nc)[1], col = 'red', add = TRUE) # add a custom legend to an arbitray plot: layout(matrix(1:2, ncol = 2), widths = c(1, lcm(2))) plot(1) .image_scale(1:10, col = sf.colors(9), key.length = lcm(8), key.pos = 4, at = 1:10) # manipulate plotting order, plot largest polygons first: p = st_polygon(list(rbind(c(0,0), c(1,0), c(1,1), c(0,1), c(0,0)))) x = st_sf(a=1:4, st_sfc(p, p * 2, p * 3, p * 4)) # plot(x, col=2:5) only shows the largest polygon! plot(x[order(st_area(x), decreasing = TRUE),], col = 2:5) # plot largest polygons first sf.colors(10) text(nc, labels = substring(nc$NAME,1,1))
Map prefix to driver
prefix_map
prefix_map
An object of class list
of length 10.
Query or manage PROJ search path and network settings
sf_proj_search_paths(paths = character(0), with_proj = NA) sf_proj_network(enable = FALSE, url = character(0)) sf_proj_pipelines( source_crs, target_crs, authority = character(0), AOI = numeric(0), Use = "NONE", grid_availability = "USED", desired_accuracy = -1, strict_containment = FALSE, axis_order_authority_compliant = st_axis_order() )
sf_proj_search_paths(paths = character(0), with_proj = NA) sf_proj_network(enable = FALSE, url = character(0)) sf_proj_pipelines( source_crs, target_crs, authority = character(0), AOI = numeric(0), Use = "NONE", grid_availability = "USED", desired_accuracy = -1, strict_containment = FALSE, axis_order_authority_compliant = st_axis_order() )
paths |
the search path to be set; omit if paths need to be queried |
with_proj |
logical; if |
enable |
logical; set this to enable ( |
url |
character; use this to specify and override the default proj network CDN |
source_crs , target_crs
|
object of class |
authority |
character; constrain output pipelines to those of authority |
AOI |
length four numeric; desired area of interest for the resulting coordinate transformations (west, south, east, north, in degrees). For an area of interest crossing the anti-meridian, west will be greater than east. |
Use |
one of "NONE", "BOTH", "INTERSECTION", "SMALLEST", indicating how AOI's of source_crs and target_crs are being used |
grid_availability |
character; one of "USED" (Grid availability is only used for sorting results. Operations where some grids are missing will be sorted last), "DISCARD" (Completely discard an operation if a required grid is missing) , "IGNORED" (Ignore grid availability at all. Results will be presented as if all grids were available.), or "AVAILABLE" (Results will be presented as if grids known to PROJ (that is registered in the grid_alternatives table of its database) were available. Used typically when networking is enabled.) |
desired_accuracy |
numeric; only return pipelines with at least this accuracy |
strict_containment |
logical; default |
axis_order_authority_compliant |
logical; if |
sf_proj_search_paths()
returns the search path (possibly after setting it)
sf_proj_network
when called without arguments returns a logical indicating whether
network search of datum grids is enabled, when called with arguments it returns a character
vector with the URL of the CDN used (or specified with url
).
sf_proj_pipelines()
returns a table with candidate coordinate transformation
pipelines along with their accuracy; NA
accuracy indicates ballpark accuracy.
Convert raw vector(s) into hexadecimal character string(s)
rawToHex(x)
rawToHex(x)
x |
raw vector, or list with raw vectors |
functions for spherical geometry, using the s2 package based on the google s2geometry.io library
sf_use_s2(use_s2) st_as_s2(x, ...) ## S3 method for class 'sf' st_as_s2(x, ...) ## S3 method for class 'sfc' st_as_s2( x, ..., oriented = getOption("s2_oriented", FALSE) || isTRUE(attr(x, "oriented")), rebuild = FALSE )
sf_use_s2(use_s2) st_as_s2(x, ...) ## S3 method for class 'sf' st_as_s2(x, ...) ## S3 method for class 'sfc' st_as_s2( x, ..., oriented = getOption("s2_oriented", FALSE) || isTRUE(attr(x, "oriented")), rebuild = FALSE )
use_s2 |
logical; if |
x |
object of class |
... |
passed on |
oriented |
logical; if |
rebuild |
logical; call s2_rebuild on the geometry (think of this as a |
st_as_s2
converts an sf
POLYGON object into a form readable by s2
.
sf_use_s2
returns the value of this variable before (re)setting it,
invisibly if use_s2
is not missing.
m = rbind(c(-1,-1), c(1,-1), c(1,1), c(-1,1), c(-1,-1)) m1 = rbind(c(-1,-1), c(1,-1), c(1,1), c(-1,1), c(-1,0), c(-1,-1)) m0 = m[5:1,] mp = st_multipolygon(list( list(m, 0.8 * m0, 0.01 * m1 + 0.9), list(0.7* m, 0.6*m0), list(0.5 * m0), list(m+2), list(m+4,(.9*m0)+4) )) sf = st_sfc(mp, mp, crs = 'EPSG:4326') s2 = st_as_s2(sf)
m = rbind(c(-1,-1), c(1,-1), c(1,1), c(-1,1), c(-1,-1)) m1 = rbind(c(-1,-1), c(1,-1), c(1,1), c(-1,1), c(-1,0), c(-1,-1)) m0 = m[5:1,] mp = st_multipolygon(list( list(m, 0.8 * m0, 0.01 * m1 + 0.9), list(0.7* m, 0.6*m0), list(0.5 * m0), list(m+2), list(m+4,(.9*m0)+4) )) sf = st_sfc(mp, mp, crs = 'EPSG:4326') s2 = st_as_s2(sf)
Create sf, which extends data.frame-like objects with a simple feature list column.
To convert a data frame object to sf
, use st_as_sf()
st_sf( ..., agr = NA_agr_, row.names, stringsAsFactors = sf_stringsAsFactors(), crs, precision, sf_column_name = NULL, check_ring_dir = FALSE, sfc_last = TRUE ) ## S3 method for class 'sf' x[i, j, ..., drop = FALSE, op = st_intersects] ## S3 method for class 'sf' print(x, ..., n = getOption("sf_max_print", default = 10))
st_sf( ..., agr = NA_agr_, row.names, stringsAsFactors = sf_stringsAsFactors(), crs, precision, sf_column_name = NULL, check_ring_dir = FALSE, sfc_last = TRUE ) ## S3 method for class 'sf' x[i, j, ..., drop = FALSE, op = st_intersects] ## S3 method for class 'sf' print(x, ..., n = getOption("sf_max_print", default = 10))
... |
column elements to be binded into an |
agr |
character vector; see details below. |
row.names |
row.names for the created |
stringsAsFactors |
logical; see st_read |
crs |
coordinate reference system, something suitable as input to st_crs |
precision |
numeric; see st_as_binary |
sf_column_name |
character; name of the active list-column with simple feature geometries; in case
there is more than one and |
check_ring_dir |
see st_read |
sfc_last |
logical; if |
x |
object of class |
i |
record selection, see [.data.frame, or a |
j |
variable selection, see [.data.frame |
drop |
logical, default |
op |
function; geometrical binary predicate function to apply when |
n |
maximum number of features to print; can be set globally by |
agr
, attribute-geometry-relationship, specifies for each non-geometry attribute column how it relates to the geometry, and can have one of following values: "constant", "aggregate", "identity". "constant" is used for attributes that are constant throughout the geometry (e.g. land use), "aggregate" where the attribute is an aggregate value over the geometry (e.g. population density or population count), "identity" when the attributes uniquely identifies the geometry of particular "thing", such as a building ID or a city name. The default value, NA_agr_
, implies we don't know.
When a single value is provided to agr
, it is cascaded across all input columns; otherwise, a named vector like c(feature1='constant', ...)
will set agr
value to 'constant'
for the input column named feature1
. See demo(nc)
for a worked example of this.
When confronted with a data.frame-like object, st_sf
will try to find a geometry column of class sfc
, and otherwise try to convert list-columns when available into a geometry column, using st_as_sfc.
[.sf
will return a data.frame
or vector if the geometry column (of class sfc
) is dropped (drop=TRUE
), an sfc
object if only the geometry column is selected, and otherwise return an sf
object; see also [.data.frame; for [.sf
...
arguments are passed to op
.
g = st_sfc(st_point(1:2)) st_sf(a=3,g) st_sf(g, a=3) st_sf(a=3, st_sfc(st_point(1:2))) # better to name it! # create empty structure with preallocated empty geometries: nrows <- 10 geometry = st_sfc(lapply(1:nrows, function(x) st_geometrycollection())) df <- st_sf(id = 1:nrows, geometry = geometry) g = st_sfc(st_point(1:2), st_point(3:4)) s = st_sf(a=3:4, g) s[1,] class(s[1,]) s[,1] class(s[,1]) s[,2] class(s[,2]) g = st_sf(a=2:3, g) pol = st_sfc(st_polygon(list(cbind(c(0,3,3,0,0),c(0,0,3,3,0))))) h = st_sf(r = 5, pol) g[h,] h[g,]
g = st_sfc(st_point(1:2)) st_sf(a=3,g) st_sf(g, a=3) st_sf(a=3, st_sfc(st_point(1:2))) # better to name it! # create empty structure with preallocated empty geometries: nrows <- 10 geometry = st_sfc(lapply(1:nrows, function(x) st_geometrycollection())) df <- st_sf(id = 1:nrows, geometry = geometry) g = st_sfc(st_point(1:2), st_point(3:4)) s = st_sf(a=3:4, g) s[1,] class(s[1,]) s[,1] class(s[,1]) s[,2] class(s[,2]) g = st_sf(a=2:3, g) pol = st_sfc(st_polygon(list(cbind(c(0,3,3,0,0),c(0,0,3,3,0))))) h = st_sf(r = 5, pol) g[h,] h[g,]
Provide the external dependencies versions of the libraries linked to sf
sf_extSoftVersion()
sf_extSoftVersion()
directly transform a set of coordinates
sf_add_proj_units() sf_project( from = character(0), to = character(0), pts, keep = FALSE, warn = TRUE, authority_compliant = st_axis_order() )
sf_add_proj_units() sf_project( from = character(0), to = character(0), pts, keep = FALSE, warn = TRUE, authority_compliant = st_axis_order() )
from |
character description of source CRS, or object of class |
to |
character description of target CRS, or object of class |
pts |
two-, three- or four-column numeric matrix, or object that can be coerced into a matrix; columns 3 and 4 contain z and t values. |
keep |
logical value controlling the handling of unprojectable points. If
|
warn |
logical; if |
authority_compliant |
logical; |
sf_add_proj_units
loads the PROJ units link
, us_in
, ind_yd
, ind_ft
, and ind_ch
into the udunits database, and returns TRUE
invisibly on success.
two-column numeric matrix with transformed/converted coordinates, returning invalid values as Inf
sf_add_proj_units()
sf_add_proj_units()
Create simple feature geometry list column, set class, and add coordinate reference system and precision.
For data.frame alternatives see st_sf()
. To convert a foreign object to sfc
, see st_as_sfc()
st_sfc( ..., crs = NA_crs_, precision = 0, check_ring_dir = FALSE, dim, recompute_bbox = FALSE, oriented = NA ) ## S3 method for class 'sfc' x[i, j, ..., op = st_intersects]
st_sfc( ..., crs = NA_crs_, precision = 0, check_ring_dir = FALSE, dim, recompute_bbox = FALSE, oriented = NA ) ## S3 method for class 'sfc' x[i, j, ..., op = st_intersects]
... |
zero or more simple feature geometries (objects of class |
crs |
coordinate reference system: integer with the EPSG code, or character with proj4string |
precision |
numeric; see st_as_binary |
check_ring_dir |
see st_read |
dim |
character; if this function is called without valid geometries, this argument may carry the right dimension to set empty geometries |
recompute_bbox |
logical; use |
oriented |
logical; if |
x |
object of class |
i |
record selection. Might also be an |
j |
ignored if |
op |
function, geometrical binary predicate function to apply when
|
A simple feature geometry list-column is a list of class
c("stc_TYPE", "sfc")
which most often contains objects of identical type;
in case of a mix of types or an empty set, TYPE
is set to the
superclass GEOMETRY
.
if x
has a dim
attribute (i.e. is an array
or matrix
) then op
cannot be used.
an object of class sfc
, which is a classed list-column with simple feature geometries.
pt1 = st_point(c(0,1)) pt2 = st_point(c(1,1)) (sfc = st_sfc(pt1, pt2)) sfc[sfc[1], op = st_is_within_distance, dist = 0.5] d = st_sf(data.frame(a=1:2, geom=sfc))
pt1 = st_point(c(0,1)) pt2 = st_point(c(1,1)) (sfc = st_sfc(pt1, pt2)) sfc[sfc[1], op = st_is_within_distance, dist = 0.5] d = st_sf(data.frame(a=1:2, geom=sfc))
Methods for dealing with sparse geometry binary predicate lists
## S3 method for class 'sgbp' print(x, ..., n = 10, max_nb = 10) ## S3 method for class 'sgbp' t(x) ## S3 method for class 'sgbp' as.matrix(x, ...) ## S3 method for class 'sgbp' dim(x) ## S3 method for class 'sgbp' Ops(e1, e2) ## S3 method for class 'sgbp' as.data.frame(x, ...)
## S3 method for class 'sgbp' print(x, ..., n = 10, max_nb = 10) ## S3 method for class 'sgbp' t(x) ## S3 method for class 'sgbp' as.matrix(x, ...) ## S3 method for class 'sgbp' dim(x) ## S3 method for class 'sgbp' Ops(e1, e2) ## S3 method for class 'sgbp' as.data.frame(x, ...)
x |
object of class |
... |
ignored |
n |
integer; maximum number of items to print |
max_nb |
integer; maximum number of neighbours to print for each item |
e1 |
object of class |
e2 |
object of class |
sgbp
are sparse matrices, stored as a list with integer vectors holding the ordered TRUE
indices of each row. This means that for a dense, matrix
Q
and a list L
, if Q[i,j]
is TRUE
then is an element of
L[[i]]
. Reversed: when is the value of
L[[i]][j]
, then Q[i,k]
is TRUE
.
==
compares only the dimension and index values, not the attributes of two sgbp
object; use identical
to check for equality of everything.
Create simple feature from a numeric vector, matrix or list
st_point(x = c(NA_real_, NA_real_), dim = "XYZ") st_multipoint(x = matrix(numeric(0), 0, 2), dim = "XYZ") st_linestring(x = matrix(numeric(0), 0, 2), dim = "XYZ") st_polygon(x = list(), dim = if (length(x)) "XYZ" else "XY") st_multilinestring(x = list(), dim = if (length(x)) "XYZ" else "XY") st_multipolygon(x = list(), dim = if (length(x)) "XYZ" else "XY") st_geometrycollection(x = list(), dims = "XY") ## S3 method for class 'sfg' print(x, ..., width = 0) ## S3 method for class 'sfg' head(x, n = 10L, ...) ## S3 method for class 'sfg' format(x, ..., width = 30) ## S3 method for class 'sfg' c(..., recursive = FALSE, flatten = TRUE) ## S3 method for class 'sfg' as.matrix(x, ...)
st_point(x = c(NA_real_, NA_real_), dim = "XYZ") st_multipoint(x = matrix(numeric(0), 0, 2), dim = "XYZ") st_linestring(x = matrix(numeric(0), 0, 2), dim = "XYZ") st_polygon(x = list(), dim = if (length(x)) "XYZ" else "XY") st_multilinestring(x = list(), dim = if (length(x)) "XYZ" else "XY") st_multipolygon(x = list(), dim = if (length(x)) "XYZ" else "XY") st_geometrycollection(x = list(), dims = "XY") ## S3 method for class 'sfg' print(x, ..., width = 0) ## S3 method for class 'sfg' head(x, n = 10L, ...) ## S3 method for class 'sfg' format(x, ..., width = 30) ## S3 method for class 'sfg' c(..., recursive = FALSE, flatten = TRUE) ## S3 method for class 'sfg' as.matrix(x, ...)
x |
for |
dim |
character, indicating dimensions: "XY", "XYZ", "XYM", or "XYZM"; only really needed for three-dimensional points (which can be either XYZ or XYM) or empty geometries; see details |
dims |
character; specify dimensionality in case of an empty (NULL) geometrycollection, in which case |
... |
objects to be pasted together into a single simple feature |
width |
integer; number of characters to be printed (max 30; 0 means print everything) |
n |
integer; number of elements to be selected |
recursive |
logical; ignored |
flatten |
logical; if |
"XYZ" refers to coordinates where the third dimension represents altitude, "XYM" refers to three-dimensional coordinates where the third dimension refers to something else ("M" for measure); checking of the sanity of x
may be only partial.
When flatten=TRUE
, this method may merge points into a multipoint structure, and may not preserve order, and hence cannot be reverted. When given fish, it returns fish soup.
object of the same nature as x
, but with appropriate class attribute set
as.matrix returns the set of points that form a geometry as a single matrix, where each point is a row; use unlist(x, recursive = FALSE)
to get sets of matrices.
(p1 = st_point(c(1,2))) class(p1) st_bbox(p1) (p2 = st_point(c(1,2,3))) class(p2) (p3 = st_point(c(1,2,3), "XYM")) pts = matrix(1:10, , 2) (mp1 = st_multipoint(pts)) pts = matrix(1:15, , 3) (mp2 = st_multipoint(pts)) (mp3 = st_multipoint(pts, "XYM")) pts = matrix(1:20, , 4) (mp4 = st_multipoint(pts)) pts = matrix(1:10, , 2) (ls1 = st_linestring(pts)) pts = matrix(1:15, , 3) (ls2 = st_linestring(pts)) (ls3 = st_linestring(pts, "XYM")) pts = matrix(1:20, , 4) (ls4 = st_linestring(pts)) outer = matrix(c(0,0,10,0,10,10,0,10,0,0),ncol=2, byrow=TRUE) hole1 = matrix(c(1,1,1,2,2,2,2,1,1,1),ncol=2, byrow=TRUE) hole2 = matrix(c(5,5,5,6,6,6,6,5,5,5),ncol=2, byrow=TRUE) pts = list(outer, hole1, hole2) (ml1 = st_multilinestring(pts)) pts3 = lapply(pts, function(x) cbind(x, 0)) (ml2 = st_multilinestring(pts3)) (ml3 = st_multilinestring(pts3, "XYM")) pts4 = lapply(pts3, function(x) cbind(x, 0)) (ml4 = st_multilinestring(pts4)) outer = matrix(c(0,0,10,0,10,10,0,10,0,0),ncol=2, byrow=TRUE) hole1 = matrix(c(1,1,1,2,2,2,2,1,1,1),ncol=2, byrow=TRUE) hole2 = matrix(c(5,5,5,6,6,6,6,5,5,5),ncol=2, byrow=TRUE) pts = list(outer, hole1, hole2) (pl1 = st_polygon(pts)) pts3 = lapply(pts, function(x) cbind(x, 0)) (pl2 = st_polygon(pts3)) (pl3 = st_polygon(pts3, "XYM")) pts4 = lapply(pts3, function(x) cbind(x, 0)) (pl4 = st_polygon(pts4)) pol1 = list(outer, hole1, hole2) pol2 = list(outer + 12, hole1 + 12) pol3 = list(outer + 24) mp = list(pol1,pol2,pol3) (mp1 = st_multipolygon(mp)) pts3 = lapply(mp, function(x) lapply(x, function(y) cbind(y, 0))) (mp2 = st_multipolygon(pts3)) (mp3 = st_multipolygon(pts3, "XYM")) pts4 = lapply(mp2, function(x) lapply(x, function(y) cbind(y, 0))) (mp4 = st_multipolygon(pts4)) (gc = st_geometrycollection(list(p1, ls1, pl1, mp1))) st_geometrycollection() # empty geometry c(st_point(1:2), st_point(5:6)) c(st_point(1:2), st_multipoint(matrix(5:8,2))) c(st_multipoint(matrix(1:4,2)), st_multipoint(matrix(5:8,2))) c(st_linestring(matrix(1:6,3)), st_linestring(matrix(11:16,3))) c(st_multilinestring(list(matrix(1:6,3))), st_multilinestring(list(matrix(11:16,3)))) pl = list(rbind(c(0,0), c(1,0), c(1,1), c(0,1), c(0,0))) c(st_polygon(pl), st_polygon(pl)) c(st_polygon(pl), st_multipolygon(list(pl))) c(st_linestring(matrix(1:6,3)), st_point(1:2)) c(st_geometrycollection(list(st_point(1:2), st_linestring(matrix(1:6,3)))), st_geometrycollection(list(st_multilinestring(list(matrix(11:16,3)))))) c(st_geometrycollection(list(st_point(1:2), st_linestring(matrix(1:6,3)))), st_multilinestring(list(matrix(11:16,3))), st_point(5:6), st_geometrycollection(list(st_point(10:11))))
(p1 = st_point(c(1,2))) class(p1) st_bbox(p1) (p2 = st_point(c(1,2,3))) class(p2) (p3 = st_point(c(1,2,3), "XYM")) pts = matrix(1:10, , 2) (mp1 = st_multipoint(pts)) pts = matrix(1:15, , 3) (mp2 = st_multipoint(pts)) (mp3 = st_multipoint(pts, "XYM")) pts = matrix(1:20, , 4) (mp4 = st_multipoint(pts)) pts = matrix(1:10, , 2) (ls1 = st_linestring(pts)) pts = matrix(1:15, , 3) (ls2 = st_linestring(pts)) (ls3 = st_linestring(pts, "XYM")) pts = matrix(1:20, , 4) (ls4 = st_linestring(pts)) outer = matrix(c(0,0,10,0,10,10,0,10,0,0),ncol=2, byrow=TRUE) hole1 = matrix(c(1,1,1,2,2,2,2,1,1,1),ncol=2, byrow=TRUE) hole2 = matrix(c(5,5,5,6,6,6,6,5,5,5),ncol=2, byrow=TRUE) pts = list(outer, hole1, hole2) (ml1 = st_multilinestring(pts)) pts3 = lapply(pts, function(x) cbind(x, 0)) (ml2 = st_multilinestring(pts3)) (ml3 = st_multilinestring(pts3, "XYM")) pts4 = lapply(pts3, function(x) cbind(x, 0)) (ml4 = st_multilinestring(pts4)) outer = matrix(c(0,0,10,0,10,10,0,10,0,0),ncol=2, byrow=TRUE) hole1 = matrix(c(1,1,1,2,2,2,2,1,1,1),ncol=2, byrow=TRUE) hole2 = matrix(c(5,5,5,6,6,6,6,5,5,5),ncol=2, byrow=TRUE) pts = list(outer, hole1, hole2) (pl1 = st_polygon(pts)) pts3 = lapply(pts, function(x) cbind(x, 0)) (pl2 = st_polygon(pts3)) (pl3 = st_polygon(pts3, "XYM")) pts4 = lapply(pts3, function(x) cbind(x, 0)) (pl4 = st_polygon(pts4)) pol1 = list(outer, hole1, hole2) pol2 = list(outer + 12, hole1 + 12) pol3 = list(outer + 24) mp = list(pol1,pol2,pol3) (mp1 = st_multipolygon(mp)) pts3 = lapply(mp, function(x) lapply(x, function(y) cbind(y, 0))) (mp2 = st_multipolygon(pts3)) (mp3 = st_multipolygon(pts3, "XYM")) pts4 = lapply(mp2, function(x) lapply(x, function(y) cbind(y, 0))) (mp4 = st_multipolygon(pts4)) (gc = st_geometrycollection(list(p1, ls1, pl1, mp1))) st_geometrycollection() # empty geometry c(st_point(1:2), st_point(5:6)) c(st_point(1:2), st_multipoint(matrix(5:8,2))) c(st_multipoint(matrix(1:4,2)), st_multipoint(matrix(5:8,2))) c(st_linestring(matrix(1:6,3)), st_linestring(matrix(11:16,3))) c(st_multilinestring(list(matrix(1:6,3))), st_multilinestring(list(matrix(11:16,3)))) pl = list(rbind(c(0,0), c(1,0), c(1,1), c(0,1), c(0,0))) c(st_polygon(pl), st_polygon(pl)) c(st_polygon(pl), st_multipolygon(list(pl))) c(st_linestring(matrix(1:6,3)), st_point(1:2)) c(st_geometrycollection(list(st_point(1:2), st_linestring(matrix(1:6,3)))), st_geometrycollection(list(st_multilinestring(list(matrix(11:16,3)))))) c(st_geometrycollection(list(st_point(1:2), st_linestring(matrix(1:6,3)))), st_multilinestring(list(matrix(11:16,3))), st_point(5:6), st_geometrycollection(list(st_point(10:11))))
sf
objectget or set relation_to_geometry attribute of an sf
object
NA_agr_ st_agr(x, ...) st_agr(x) <- value st_set_agr(x, value)
NA_agr_ st_agr(x, ...) st_agr(x) <- value st_set_agr(x, value)
x |
object of class |
... |
ignored |
value |
character, or factor with appropriate levels; if named, names should correspond to the non-geometry list-column columns of |
An object of class factor
of length 1.
NA_agr_
is the agr
object with a missing value.
Convert sfc object to an WKB object
st_as_binary(x, ...) ## S3 method for class 'sfc' st_as_binary( x, ..., EWKB = FALSE, endian = .Platform$endian, pureR = FALSE, precision = attr(x, "precision"), hex = FALSE ) ## S3 method for class 'sfg' st_as_binary( x, ..., endian = .Platform$endian, EWKB = FALSE, pureR = FALSE, hex = FALSE, srid = 0 )
st_as_binary(x, ...) ## S3 method for class 'sfc' st_as_binary( x, ..., EWKB = FALSE, endian = .Platform$endian, pureR = FALSE, precision = attr(x, "precision"), hex = FALSE ) ## S3 method for class 'sfg' st_as_binary( x, ..., endian = .Platform$endian, EWKB = FALSE, pureR = FALSE, hex = FALSE, srid = 0 )
x |
object to convert |
... |
ignored |
EWKB |
logical; use EWKB (PostGIS), or (default) ISO-WKB? |
endian |
character; either "big" or "little"; default: use that of platform |
pureR |
logical; use pure R solution, or C++? |
precision |
numeric; if zero, do not modify; to reduce precision: negative values convert to float (4-byte real); positive values convert to round(x*precision)/precision. See details. |
hex |
logical; return as (unclassed) hexadecimal encoded character vector? |
srid |
integer; override srid (can be used when the srid is unavailable locally). |
st_as_binary
is called on sfc objects on their way to the GDAL or GEOS libraries, and hence does rounding (if requested) on the fly before e.g. computing spatial predicates like st_intersects. The examples show a round-trip of an sfc
to and from binary.
For the precision model used, see also https://locationtech.github.io/jts/javadoc/org/locationtech/jts/geom/PrecisionModel.html. There, it is written that: “... to specify 3 decimal places of precision, use a scale factor of 1000. To specify -3 decimal places of precision (i.e. rounding to the nearest 1000), use a scale factor of 0.001.”. Note that ALL coordinates, so also Z or M values (if present) are affected.
# examples of setting precision: st_point(c(1/3, 1/6)) %>% st_sfc(precision = 1000) %>% st_as_binary %>% st_as_sfc st_point(c(1/3, 1/6)) %>% st_sfc(precision = 100) %>% st_as_binary %>% st_as_sfc st_point(1e6 * c(1/3, 1/6)) %>% st_sfc(precision = 0.01) %>% st_as_binary %>% st_as_sfc st_point(1e6 * c(1/3, 1/6)) %>% st_sfc(precision = 0.001) %>% st_as_binary %>% st_as_sfc
# examples of setting precision: st_point(c(1/3, 1/6)) %>% st_sfc(precision = 1000) %>% st_as_binary %>% st_as_sfc st_point(c(1/3, 1/6)) %>% st_sfc(precision = 100) %>% st_as_binary %>% st_as_sfc st_point(1e6 * c(1/3, 1/6)) %>% st_sfc(precision = 0.01) %>% st_as_binary %>% st_as_sfc st_point(1e6 * c(1/3, 1/6)) %>% st_sfc(precision = 0.001) %>% st_as_binary %>% st_as_sfc
Convert sf* object to an grid graphics object (grob)
st_as_grob(x, ...)
st_as_grob(x, ...)
x |
object to be converted into an object class |
... |
passed on to the xxxGrob function, e.g. |
Convert foreign object to an sf object
st_as_sf(x, ...) ## S3 method for class 'data.frame' st_as_sf( x, ..., agr = NA_agr_, coords, wkt, dim = "XYZ", remove = TRUE, na.fail = TRUE, sf_column_name = NULL ) ## S3 method for class 'sf' st_as_sf(x, ...) ## S3 method for class 'sfc' st_as_sf(x, ...) ## S3 method for class 'Spatial' st_as_sf(x, ...) ## S3 method for class 'map' st_as_sf(x, ..., fill = TRUE, group = TRUE) ## S3 method for class 'ppp' st_as_sf(x, ...) ## S3 method for class 'psp' st_as_sf(x, ...) ## S3 method for class 'lpp' st_as_sf(x, ...) ## S3 method for class 's2_geography' st_as_sf(x, ..., crs = st_crs(4326))
st_as_sf(x, ...) ## S3 method for class 'data.frame' st_as_sf( x, ..., agr = NA_agr_, coords, wkt, dim = "XYZ", remove = TRUE, na.fail = TRUE, sf_column_name = NULL ) ## S3 method for class 'sf' st_as_sf(x, ...) ## S3 method for class 'sfc' st_as_sf(x, ...) ## S3 method for class 'Spatial' st_as_sf(x, ...) ## S3 method for class 'map' st_as_sf(x, ..., fill = TRUE, group = TRUE) ## S3 method for class 'ppp' st_as_sf(x, ...) ## S3 method for class 'psp' st_as_sf(x, ...) ## S3 method for class 'lpp' st_as_sf(x, ...) ## S3 method for class 's2_geography' st_as_sf(x, ..., crs = st_crs(4326))
x |
object to be converted into an object class |
... |
passed on to st_sf, might included named arguments |
agr |
character vector; see details section of st_sf |
coords |
in case of point data: names or numbers of the numeric columns holding coordinates |
wkt |
name or number of the character column that holds WKT encoded geometries |
dim |
specify what 3- or 4-dimensional points reflect: passed on to st_point (only when argument coords is given) |
remove |
logical; when coords or wkt is given, remove these columns from data.frame? |
na.fail |
logical; if |
sf_column_name |
character; name of the active list-column with simple feature geometries; in case
there is more than one and |
fill |
logical; the value for |
group |
logical; if |
crs |
coordinate reference system to be assigned; object of class |
setting argument wkt
annihilates the use of argument coords
. If x
contains a column called "geometry", coords
will result in overwriting of this column by the sfc geometry list-column. Setting wkt
will replace this column with the geometry list-column, unless remove
is FALSE
.
If coords
has length 4, and dim
is not XYZM
, the four columns are taken as the xmin, ymin, xmax, ymax corner coordinates of a rectangle, and polygons are returned.
pt1 = st_point(c(0,1)) pt2 = st_point(c(1,1)) st_sfc(pt1, pt2) d = data.frame(a = 1:2) d$geom = st_sfc(pt1, pt2) df = st_as_sf(d) d$geom = c("POINT(0 0)", "POINT(0 1)") df = st_as_sf(d, wkt = "geom") d$geom2 = st_sfc(pt1, pt2) st_as_sf(d) # should warn if (require(sp, quietly = TRUE)) { data(meuse, package = "sp") meuse_sf = st_as_sf(meuse, coords = c("x", "y"), crs = 28992, agr = "constant") meuse_sf[1:3,] summary(meuse_sf) } if (require(sp, quietly = TRUE)) { x = rbind(c(-1,-1), c(1,-1), c(1,1), c(-1,1), c(-1,-1)) x1 = 0.1 * x + 0.1 x2 = 0.1 * x + 0.4 x3 = 0.1 * x + 0.7 y = x + 3 y1 = x1 + 3 y3 = x3 + 3 m = matrix(c(3, 0), 5, 2, byrow = TRUE) z = x + m z1 = x1 + m z2 = x2 + m z3 = x3 + m p1 = Polygons(list( Polygon(x[5:1,]), Polygon(x2), Polygon(x3), Polygon(y[5:1,]), Polygon(y1), Polygon(x1), Polygon(y3)), "ID1") p2 = Polygons(list( Polygon(z[5:1,]), Polygon(z2), Polygon(z3), Polygon(z1)), "ID2") r = SpatialPolygons(list(p1,p2)) a = suppressWarnings(st_as_sf(r)) summary(a) demo(meuse, ask = FALSE, echo = FALSE) summary(st_as_sf(meuse)) summary(st_as_sf(meuse.grid)) summary(st_as_sf(meuse.area)) summary(st_as_sf(meuse.riv)) summary(st_as_sf(as(meuse.riv, "SpatialLines"))) pol.grd = as(meuse.grid, "SpatialPolygonsDataFrame") # summary(st_as_sf(pol.grd)) # summary(st_as_sf(as(pol.grd, "SpatialLinesDataFrame"))) } if (require(spatstat.geom)) { g = st_as_sf(gorillas) # select only the points: g[st_is(g, "POINT"),] } if (require(spatstat.linnet)) { data(chicago) plot(st_as_sf(chicago)["label"]) plot(st_as_sf(chicago)[-1,"label"]) }
pt1 = st_point(c(0,1)) pt2 = st_point(c(1,1)) st_sfc(pt1, pt2) d = data.frame(a = 1:2) d$geom = st_sfc(pt1, pt2) df = st_as_sf(d) d$geom = c("POINT(0 0)", "POINT(0 1)") df = st_as_sf(d, wkt = "geom") d$geom2 = st_sfc(pt1, pt2) st_as_sf(d) # should warn if (require(sp, quietly = TRUE)) { data(meuse, package = "sp") meuse_sf = st_as_sf(meuse, coords = c("x", "y"), crs = 28992, agr = "constant") meuse_sf[1:3,] summary(meuse_sf) } if (require(sp, quietly = TRUE)) { x = rbind(c(-1,-1), c(1,-1), c(1,1), c(-1,1), c(-1,-1)) x1 = 0.1 * x + 0.1 x2 = 0.1 * x + 0.4 x3 = 0.1 * x + 0.7 y = x + 3 y1 = x1 + 3 y3 = x3 + 3 m = matrix(c(3, 0), 5, 2, byrow = TRUE) z = x + m z1 = x1 + m z2 = x2 + m z3 = x3 + m p1 = Polygons(list( Polygon(x[5:1,]), Polygon(x2), Polygon(x3), Polygon(y[5:1,]), Polygon(y1), Polygon(x1), Polygon(y3)), "ID1") p2 = Polygons(list( Polygon(z[5:1,]), Polygon(z2), Polygon(z3), Polygon(z1)), "ID2") r = SpatialPolygons(list(p1,p2)) a = suppressWarnings(st_as_sf(r)) summary(a) demo(meuse, ask = FALSE, echo = FALSE) summary(st_as_sf(meuse)) summary(st_as_sf(meuse.grid)) summary(st_as_sf(meuse.area)) summary(st_as_sf(meuse.riv)) summary(st_as_sf(as(meuse.riv, "SpatialLines"))) pol.grd = as(meuse.grid, "SpatialPolygonsDataFrame") # summary(st_as_sf(pol.grd)) # summary(st_as_sf(as(pol.grd, "SpatialLinesDataFrame"))) } if (require(spatstat.geom)) { g = st_as_sf(gorillas) # select only the points: g[st_is(g, "POINT"),] } if (require(spatstat.linnet)) { data(chicago) plot(st_as_sf(chicago)["label"]) plot(st_as_sf(chicago)[-1,"label"]) }
Convert foreign geometry object to an sfc object
## S3 method for class 'pq_geometry' st_as_sfc( x, ..., EWKB = TRUE, spatialite = FALSE, pureR = FALSE, crs = NA_crs_ ) ## S3 method for class 'list' st_as_sfc(x, ..., crs = NA_crs_) ## S3 method for class 'blob' st_as_sfc(x, ...) ## S3 method for class 'bbox' st_as_sfc(x, ...) ## S3 method for class 'WKB' st_as_sfc( x, ..., EWKB = FALSE, spatialite = FALSE, pureR = FALSE, crs = NA_crs_ ) ## S3 method for class 'raw' st_as_sfc(x, ...) ## S3 method for class 'character' st_as_sfc(x, crs = NA_integer_, ..., GeoJSON = FALSE) ## S3 method for class 'factor' st_as_sfc(x, ...) st_as_sfc(x, ...) ## S3 method for class 'SpatialPoints' st_as_sfc(x, ..., precision = 0) ## S3 method for class 'SpatialPixels' st_as_sfc(x, ..., precision = 0) ## S3 method for class 'SpatialMultiPoints' st_as_sfc(x, ..., precision = 0) ## S3 method for class 'SpatialLines' st_as_sfc(x, ..., precision = 0, forceMulti = FALSE) ## S3 method for class 'SpatialPolygons' st_as_sfc(x, ..., precision = 0, forceMulti = FALSE) ## S3 method for class 'map' st_as_sfc(x, ...) ## S3 method for class 's2_geography' st_as_sfc( x, ..., crs = st_crs(4326), endian = match(.Platform$endian, c("big", "little")) - 1L )
## S3 method for class 'pq_geometry' st_as_sfc( x, ..., EWKB = TRUE, spatialite = FALSE, pureR = FALSE, crs = NA_crs_ ) ## S3 method for class 'list' st_as_sfc(x, ..., crs = NA_crs_) ## S3 method for class 'blob' st_as_sfc(x, ...) ## S3 method for class 'bbox' st_as_sfc(x, ...) ## S3 method for class 'WKB' st_as_sfc( x, ..., EWKB = FALSE, spatialite = FALSE, pureR = FALSE, crs = NA_crs_ ) ## S3 method for class 'raw' st_as_sfc(x, ...) ## S3 method for class 'character' st_as_sfc(x, crs = NA_integer_, ..., GeoJSON = FALSE) ## S3 method for class 'factor' st_as_sfc(x, ...) st_as_sfc(x, ...) ## S3 method for class 'SpatialPoints' st_as_sfc(x, ..., precision = 0) ## S3 method for class 'SpatialPixels' st_as_sfc(x, ..., precision = 0) ## S3 method for class 'SpatialMultiPoints' st_as_sfc(x, ..., precision = 0) ## S3 method for class 'SpatialLines' st_as_sfc(x, ..., precision = 0, forceMulti = FALSE) ## S3 method for class 'SpatialPolygons' st_as_sfc(x, ..., precision = 0, forceMulti = FALSE) ## S3 method for class 'map' st_as_sfc(x, ...) ## S3 method for class 's2_geography' st_as_sfc( x, ..., crs = st_crs(4326), endian = match(.Platform$endian, c("big", "little")) - 1L )
x |
object to convert |
... |
further arguments |
EWKB |
logical; if |
spatialite |
logical; if |
pureR |
logical; if |
crs |
coordinate reference system to be assigned; object of class |
GeoJSON |
logical; if |
precision |
precision value; see st_as_binary |
forceMulti |
logical; if |
endian |
integer; 0 or 1: defaults to the endian of the native machine |
When converting from WKB, the object x
is either a character vector such as typically obtained from PostGIS (either with leading "0x" or without), or a list with raw vectors representing the features in binary (raw) form.
If x
is a character vector, it should be a vector containing
well-known-text, or
Postgis EWKT or GeoJSON representations of a single geometry for each vector element.
If x
is a factor
, it is converted to character
.
wkb = structure(list("01010000204071000000000000801A064100000000AC5C1441"), class = "WKB") st_as_sfc(wkb, EWKB = TRUE) wkb = structure(list("0x01010000204071000000000000801A064100000000AC5C1441"), class = "WKB") st_as_sfc(wkb, EWKB = TRUE) st_as_sfc(st_as_binary(st_sfc(st_point(0:1)))[[1]], crs = 4326) st_as_sfc("SRID=3978;LINESTRING(1663106 -105415,1664320 -104617)")
wkb = structure(list("01010000204071000000000000801A064100000000AC5C1441"), class = "WKB") st_as_sfc(wkb, EWKB = TRUE) wkb = structure(list("0x01010000204071000000000000801A064100000000AC5C1441"), class = "WKB") st_as_sfc(wkb, EWKB = TRUE) st_as_sfc(st_as_binary(st_sfc(st_point(0:1)))[[1]], crs = 4326) st_as_sfc("SRID=3978;LINESTRING(1663106 -105415,1664320 -104617)")
Return Well-known Text representation of simple feature geometry or coordinate reference system
## S3 method for class 'crs' st_as_text(x, ..., projjson = FALSE, pretty = FALSE) st_as_text(x, ...) ## S3 method for class 'sfg' st_as_text(x, ...) ## S3 method for class 'sfc' st_as_text(x, ..., EWKT = FALSE)
## S3 method for class 'crs' st_as_text(x, ..., projjson = FALSE, pretty = FALSE) st_as_text(x, ...) ## S3 method for class 'sfg' st_as_text(x, ...) ## S3 method for class 'sfc' st_as_text(x, ..., EWKT = FALSE)
x |
object of class |
... |
modifiers; in particular |
projjson |
logical; if TRUE, return projjson form (requires GDAL 3.1 and PROJ 6.2), else return well-known-text form |
pretty |
logical; if TRUE, print human-readable well-known-text representation of a coordinate reference system |
EWKT |
logical; if TRUE, print SRID=xxx; before the WKT string if |
The returned WKT representation of simple feature geometry conforms to the simple features access specification and extensions (known as EWKT, supported by PostGIS and other simple features implementations for addition of a SRID to a WKT string).
To improve conversion performance, the lwgeom package can be used (it must be installed
beforehand) and set the Sys.setenv("LWGEOM_WKT" = "true")
environment variable. This
will also result in faster printing of complex geometries. Note that the representation as WKT is
different from the sf package and may cause reproducibility problems. An alternative solution is
to use the lwgeom::st_astext()
or wk::as_wkt()
functions.
st_as_text(st_point(1:2)) st_as_text(st_sfc(st_point(c(-90,40)), crs = 4326), EWKT = TRUE)
st_as_text(st_point(1:2)) st_as_text(st_sfc(st_point(c(-90,40)), crs = 4326), EWKT = TRUE)
Return bounding of a simple feature or simple feature set
## S3 method for class 'bbox' is.na(x) st_bbox(obj, ...) ## S3 method for class 'POINT' st_bbox(obj, ...) ## S3 method for class 'MULTIPOINT' st_bbox(obj, ...) ## S3 method for class 'LINESTRING' st_bbox(obj, ...) ## S3 method for class 'POLYGON' st_bbox(obj, ...) ## S3 method for class 'MULTILINESTRING' st_bbox(obj, ...) ## S3 method for class 'MULTIPOLYGON' st_bbox(obj, ...) ## S3 method for class 'GEOMETRYCOLLECTION' st_bbox(obj, ...) ## S3 method for class 'MULTISURFACE' st_bbox(obj, ...) ## S3 method for class 'MULTICURVE' st_bbox(obj, ...) ## S3 method for class 'CURVEPOLYGON' st_bbox(obj, ...) ## S3 method for class 'COMPOUNDCURVE' st_bbox(obj, ...) ## S3 method for class 'POLYHEDRALSURFACE' st_bbox(obj, ...) ## S3 method for class 'TIN' st_bbox(obj, ...) ## S3 method for class 'TRIANGLE' st_bbox(obj, ...) ## S3 method for class 'CIRCULARSTRING' st_bbox(obj, ...) ## S3 method for class 'sfc' st_bbox(obj, ...) ## S3 method for class 'sf' st_bbox(obj, ...) ## S3 method for class 'Spatial' st_bbox(obj, ...) ## S3 method for class 'Raster' st_bbox(obj, ...) ## S3 method for class 'Extent' st_bbox(obj, ..., crs = NA_crs_) ## S3 method for class 'numeric' st_bbox(obj, ..., crs = NA_crs_) NA_bbox_ FULL_bbox_ ## S3 method for class 'bbox' format(x, ...)
## S3 method for class 'bbox' is.na(x) st_bbox(obj, ...) ## S3 method for class 'POINT' st_bbox(obj, ...) ## S3 method for class 'MULTIPOINT' st_bbox(obj, ...) ## S3 method for class 'LINESTRING' st_bbox(obj, ...) ## S3 method for class 'POLYGON' st_bbox(obj, ...) ## S3 method for class 'MULTILINESTRING' st_bbox(obj, ...) ## S3 method for class 'MULTIPOLYGON' st_bbox(obj, ...) ## S3 method for class 'GEOMETRYCOLLECTION' st_bbox(obj, ...) ## S3 method for class 'MULTISURFACE' st_bbox(obj, ...) ## S3 method for class 'MULTICURVE' st_bbox(obj, ...) ## S3 method for class 'CURVEPOLYGON' st_bbox(obj, ...) ## S3 method for class 'COMPOUNDCURVE' st_bbox(obj, ...) ## S3 method for class 'POLYHEDRALSURFACE' st_bbox(obj, ...) ## S3 method for class 'TIN' st_bbox(obj, ...) ## S3 method for class 'TRIANGLE' st_bbox(obj, ...) ## S3 method for class 'CIRCULARSTRING' st_bbox(obj, ...) ## S3 method for class 'sfc' st_bbox(obj, ...) ## S3 method for class 'sf' st_bbox(obj, ...) ## S3 method for class 'Spatial' st_bbox(obj, ...) ## S3 method for class 'Raster' st_bbox(obj, ...) ## S3 method for class 'Extent' st_bbox(obj, ..., crs = NA_crs_) ## S3 method for class 'numeric' st_bbox(obj, ..., crs = NA_crs_) NA_bbox_ FULL_bbox_ ## S3 method for class 'bbox' format(x, ...)
x |
object of class |
obj |
object to compute the bounding box from |
... |
for format.bbox, passed on to format to format individual numbers |
crs |
object of class |
An object of class bbox
of length 4.
An object of class bbox
of length 4.
NA_bbox_
represents the missing value for a bbox
object
NA_bbox_
represents the missing value for a bbox
object
a numeric vector of length four, with xmin
, ymin
, xmax
and ymax
values; if obj
is of class sf
, sfc
, Spatial
or Raster
, the object
returned has a class bbox
, an attribute crs
and a method to print the
bbox and an st_crs
method to retrieve the coordinate reference system
corresponding to obj
(and hence the bounding box). st_as_sfc has a
methods for bbox
objects to generate a polygon around the four bounding box points.
a = st_sf(a = 1:2, geom = st_sfc(st_point(0:1), st_point(1:2)), crs = 4326) st_bbox(a) st_as_sfc(st_bbox(a)) st_bbox(c(xmin = 16.1, xmax = 16.6, ymax = 48.6, ymin = 47.9), crs = st_crs(4326))
a = st_sf(a = 1:2, geom = st_sfc(st_point(0:1), st_point(1:2)), crs = 4326) st_bbox(a) st_as_sfc(st_bbox(a)) st_bbox(c(xmin = 16.1, xmax = 16.6, ymax = 48.6, ymin = 47.9), crs = st_crs(4326))
Longitudes can be broken at the antimeridian of a target central longitude
to permit plotting of (usually world) line or polygon objects centred
on the chosen central longitude. The method may only be used with
non-projected, geographical coordinates and linestring or polygon objects.
s2 is turned off internally to permit the use of a rectangular bounding
box. If the input geometries go outside [-180, 180]
degrees longitude,
the protruding geometries will also be split using the same tol=
values; in this case empty geometries will be dropped first.
st_break_antimeridian(x, lon_0 = 0, tol = 1e-04, ...) ## S3 method for class 'sf' st_break_antimeridian(x, lon_0 = 0, tol = 1e-04, ...) ## S3 method for class 'sfc' st_break_antimeridian(x, lon_0 = 0, tol = 1e-04, ...)
st_break_antimeridian(x, lon_0 = 0, tol = 1e-04, ...) ## S3 method for class 'sf' st_break_antimeridian(x, lon_0 = 0, tol = 1e-04, ...) ## S3 method for class 'sfc' st_break_antimeridian(x, lon_0 = 0, tol = 1e-04, ...)
x |
object of class |
lon_0 |
target central longitude (degrees) |
tol |
half of break width (degrees, default 0.0001) |
... |
ignored here |
if (require("maps", quietly=TRUE)) { opar = par(mfrow=c(3, 2)) wld = st_as_sf(map(fill=FALSE, interior=FALSE, plot=FALSE), fill=FALSE) for (lon_0 in c(-170, -90, -10, 10, 90, 170)) { wld |> st_break_antimeridian(lon_0=lon_0) |> st_transform(paste0("+proj=natearth +lon_0=", lon_0)) |> st_geometry() |> plot(main=lon_0) } par(opar) }
if (require("maps", quietly=TRUE)) { opar = par(mfrow=c(3, 2)) wld = st_as_sf(map(fill=FALSE, interior=FALSE, plot=FALSE), fill=FALSE) for (lon_0 in c(-170, -90, -10, 10, 90, 170)) { wld |> st_break_antimeridian(lon_0=lon_0) |> st_transform(paste0("+proj=natearth +lon_0=", lon_0)) |> st_geometry() |> plot(main=lon_0) } par(opar) }
Cast geometry to another type: either simplify, or cast explicitly
## S3 method for class 'MULTIPOLYGON' st_cast(x, to, ...) ## S3 method for class 'MULTILINESTRING' st_cast(x, to, ...) ## S3 method for class 'MULTIPOINT' st_cast(x, to, ...) ## S3 method for class 'POLYGON' st_cast(x, to, ...) ## S3 method for class 'LINESTRING' st_cast(x, to, ...) ## S3 method for class 'POINT' st_cast(x, to, ...) ## S3 method for class 'GEOMETRYCOLLECTION' st_cast(x, to, ...) ## S3 method for class 'CIRCULARSTRING' st_cast(x, to, ...) ## S3 method for class 'MULTISURFACE' st_cast(x, to, ...) ## S3 method for class 'COMPOUNDCURVE' st_cast(x, to, ...) ## S3 method for class 'MULTICURVE' st_cast(x, to, ...) ## S3 method for class 'CURVE' st_cast(x, to, ...) st_cast(x, to, ...) ## S3 method for class 'sfc' st_cast(x, to, ..., ids = seq_along(x), group_or_split = TRUE) ## S3 method for class 'sf' st_cast(x, to, ..., warn = TRUE, do_split = TRUE) ## S3 method for class 'sfc_CIRCULARSTRING' st_cast(x, to, ...)
## S3 method for class 'MULTIPOLYGON' st_cast(x, to, ...) ## S3 method for class 'MULTILINESTRING' st_cast(x, to, ...) ## S3 method for class 'MULTIPOINT' st_cast(x, to, ...) ## S3 method for class 'POLYGON' st_cast(x, to, ...) ## S3 method for class 'LINESTRING' st_cast(x, to, ...) ## S3 method for class 'POINT' st_cast(x, to, ...) ## S3 method for class 'GEOMETRYCOLLECTION' st_cast(x, to, ...) ## S3 method for class 'CIRCULARSTRING' st_cast(x, to, ...) ## S3 method for class 'MULTISURFACE' st_cast(x, to, ...) ## S3 method for class 'COMPOUNDCURVE' st_cast(x, to, ...) ## S3 method for class 'MULTICURVE' st_cast(x, to, ...) ## S3 method for class 'CURVE' st_cast(x, to, ...) st_cast(x, to, ...) ## S3 method for class 'sfc' st_cast(x, to, ..., ids = seq_along(x), group_or_split = TRUE) ## S3 method for class 'sf' st_cast(x, to, ..., warn = TRUE, do_split = TRUE) ## S3 method for class 'sfc_CIRCULARSTRING' st_cast(x, to, ...)
x |
object of class |
to |
character; target type, if missing, simplification is tried; when |
... |
ignored |
ids |
integer vector, denoting how geometries should be grouped (default: no grouping) |
group_or_split |
logical; if TRUE, group or split geometries; if FALSE, carry out a 1-1 per-geometry conversion. |
warn |
logical; if |
do_split |
logical; if |
When converting a GEOMETRYCOLLECTION to COMPOUNDCURVE, MULTISURFACE or CURVEPOLYGON, the user is responsible for the validity of the resulting object: no checks are being carried out by the software.
When converting mixed, GEOMETRY sets, it may help to first convert to the MULTI-type, see examples
the st_cast
method for sf
objects can only split geometries, e.g. cast MULTIPOINT
into multiple POINT
features. In case of splitting, attributes are repeated and a warning is issued when non-constant attributes are assigned to sub-geometries. To merge feature geometries and attribute values, use aggregate or summarise.
object of class to
if successful, or unmodified object if unsuccessful. If information gets lost while type casting, a warning is raised.
In case to
is missing, st_cast.sfc
will coerce combinations of "POINT" and "MULTIPOINT", "LINESTRING" and "MULTILINESTRING", "POLYGON" and "MULTIPOLYGON" into their "MULTI..." form, or in case all geometries are "GEOMETRYCOLLECTION" will return a list of all the contents of the "GEOMETRYCOLLECTION" objects, or else do nothing. In case to
is specified, if to
is "GEOMETRY", geometries are not converted, else, st_cast
will try to coerce all elements into to
; ids
may be specified to group e.g. "POINT" objects into a "MULTIPOINT", if not specified no grouping takes place. If e.g. a "sfc_MULTIPOINT" is cast to a "sfc_POINT", the objects are split, so no information gets lost, unless group_or_split
is FALSE
.
# example(st_read) nc = st_read(system.file("shape/nc.shp", package="sf")) mpl <- st_geometry(nc)[[4]] #st_cast(x) ## error 'argument "to" is missing, with no default' cast_all <- function(xg) { lapply(c("MULTIPOLYGON", "MULTILINESTRING", "MULTIPOINT", "POLYGON", "LINESTRING", "POINT"), function(x) st_cast(xg, x)) } st_sfc(cast_all(mpl)) ## no closing coordinates should remain for multipoint any(duplicated(unclass(st_cast(mpl, "MULTIPOINT")))) ## should be FALSE ## number of duplicated coordinates in the linestrings should equal the number of polygon rings ## (... in this case, won't always be true) sum(duplicated(do.call(rbind, unclass(st_cast(mpl, "MULTILINESTRING")))) ) == sum(unlist(lapply(mpl, length))) ## should be TRUE p1 <- structure(c(0, 1, 3, 2, 1, 0, 0, 0, 2, 4, 4, 0), .Dim = c(6L, 2L)) p2 <- structure(c(1, 1, 2, 1, 1, 2, 2, 1), .Dim = c(4L, 2L)) st_polygon(list(p1, p2)) mls <- st_cast(st_geometry(nc)[[4]], "MULTILINESTRING") st_sfc(cast_all(mls)) mpt <- st_cast(st_geometry(nc)[[4]], "MULTIPOINT") st_sfc(cast_all(mpt)) pl <- st_cast(st_geometry(nc)[[4]], "POLYGON") st_sfc(cast_all(pl)) ls <- st_cast(st_geometry(nc)[[4]], "LINESTRING") st_sfc(cast_all(ls)) pt <- st_cast(st_geometry(nc)[[4]], "POINT") ## st_sfc(cast_all(pt)) ## Error: cannot create MULTIPOLYGON from POINT st_sfc(lapply(c("POINT", "MULTIPOINT"), function(x) st_cast(pt, x))) s = st_multipoint(rbind(c(1,0))) st_cast(s, "POINT") # https://github.com/r-spatial/sf/issues/1930: pt1 <- st_point(c(0,1)) pt23 <- st_multipoint(matrix(c(1,2,3,4), ncol = 2, byrow = TRUE)) d <- st_sf(geom = st_sfc(pt1, pt23)) st_cast(d, "POINT") # will not convert the entire MULTIPOINT, and warns st_cast(d, "MULTIPOINT") %>% st_cast("POINT")
# example(st_read) nc = st_read(system.file("shape/nc.shp", package="sf")) mpl <- st_geometry(nc)[[4]] #st_cast(x) ## error 'argument "to" is missing, with no default' cast_all <- function(xg) { lapply(c("MULTIPOLYGON", "MULTILINESTRING", "MULTIPOINT", "POLYGON", "LINESTRING", "POINT"), function(x) st_cast(xg, x)) } st_sfc(cast_all(mpl)) ## no closing coordinates should remain for multipoint any(duplicated(unclass(st_cast(mpl, "MULTIPOINT")))) ## should be FALSE ## number of duplicated coordinates in the linestrings should equal the number of polygon rings ## (... in this case, won't always be true) sum(duplicated(do.call(rbind, unclass(st_cast(mpl, "MULTILINESTRING")))) ) == sum(unlist(lapply(mpl, length))) ## should be TRUE p1 <- structure(c(0, 1, 3, 2, 1, 0, 0, 0, 2, 4, 4, 0), .Dim = c(6L, 2L)) p2 <- structure(c(1, 1, 2, 1, 1, 2, 2, 1), .Dim = c(4L, 2L)) st_polygon(list(p1, p2)) mls <- st_cast(st_geometry(nc)[[4]], "MULTILINESTRING") st_sfc(cast_all(mls)) mpt <- st_cast(st_geometry(nc)[[4]], "MULTIPOINT") st_sfc(cast_all(mpt)) pl <- st_cast(st_geometry(nc)[[4]], "POLYGON") st_sfc(cast_all(pl)) ls <- st_cast(st_geometry(nc)[[4]], "LINESTRING") st_sfc(cast_all(ls)) pt <- st_cast(st_geometry(nc)[[4]], "POINT") ## st_sfc(cast_all(pt)) ## Error: cannot create MULTIPOLYGON from POINT st_sfc(lapply(c("POINT", "MULTIPOINT"), function(x) st_cast(pt, x))) s = st_multipoint(rbind(c(1,0))) st_cast(s, "POINT") # https://github.com/r-spatial/sf/issues/1930: pt1 <- st_point(c(0,1)) pt23 <- st_multipoint(matrix(c(1,2,3,4), ncol = 2, byrow = TRUE)) d <- st_sf(geom = st_sfc(pt1, pt23)) st_cast(d, "POINT") # will not convert the entire MULTIPOINT, and warns st_cast(d, "MULTIPOINT") %>% st_cast("POINT")
Mixes of POINTS and MULTIPOINTS, LINESTRING and MULTILINESTRING, POLYGON and MULTIPOLYGON are returned as MULTIPOINTS, MULTILINESTRING and MULTIPOLYGONS respectively
st_cast_sfc_default(x)
st_cast_sfc_default(x)
x |
list of geometries or simple features |
Geometries that are already MULTI* are left unchanged. Features that can't be cast to a single MULTI* geometry are return as a GEOMETRYCOLLECTION
GEOMETRY
or GEOMETRYCOLLECTION
,
return an object consisting only of elements of the specified type.Similar to ST_CollectionExtract in PostGIS. If there are no sub-geometries of the specified type, an empty geometry is returned.
st_collection_extract( x, type = c("POLYGON", "POINT", "LINESTRING"), warn = FALSE ) ## S3 method for class 'sfg' st_collection_extract( x, type = c("POLYGON", "POINT", "LINESTRING"), warn = FALSE ) ## S3 method for class 'sfc' st_collection_extract( x, type = c("POLYGON", "POINT", "LINESTRING"), warn = FALSE ) ## S3 method for class 'sf' st_collection_extract( x, type = c("POLYGON", "POINT", "LINESTRING"), warn = FALSE )
st_collection_extract( x, type = c("POLYGON", "POINT", "LINESTRING"), warn = FALSE ) ## S3 method for class 'sfg' st_collection_extract( x, type = c("POLYGON", "POINT", "LINESTRING"), warn = FALSE ) ## S3 method for class 'sfc' st_collection_extract( x, type = c("POLYGON", "POINT", "LINESTRING"), warn = FALSE ) ## S3 method for class 'sf' st_collection_extract( x, type = c("POLYGON", "POINT", "LINESTRING"), warn = FALSE )
x |
an object of class |
type |
character; one of "POLYGON", "POINT", "LINESTRING" |
warn |
logical; if |
An object having the same class as x
, with geometries
consisting only of elements of the specified type.
For sfg
objects, an sfg
object is returned if there is only
one geometry of the specified type, otherwise the geometries are combined
into an sfc
object of the relevant type. If any subgeometries in the
input are MULTI, then all of the subgeometries in the output will be MULTI.
pt <- st_point(c(1, 0)) ls <- st_linestring(matrix(c(4, 3, 0, 0), ncol = 2)) poly1 <- st_polygon(list(matrix(c(5.5, 7, 7, 6, 5.5, 0, 0, -0.5, -0.5, 0), ncol = 2))) poly2 <- st_polygon(list(matrix(c(6.6, 8, 8, 7, 6.6, 1, 1, 1.5, 1.5, 1), ncol = 2))) multipoly <- st_multipolygon(list(poly1, poly2)) i <- st_geometrycollection(list(pt, ls, poly1, poly2)) j <- st_geometrycollection(list(pt, ls, poly1, poly2, multipoly)) st_collection_extract(i, "POLYGON") st_collection_extract(i, "POINT") st_collection_extract(i, "LINESTRING") ## A GEOMETRYCOLLECTION aa <- rbind(st_sf(a=1, geom = st_sfc(i)), st_sf(a=2, geom = st_sfc(j))) ## With sf objects st_collection_extract(aa, "POLYGON") st_collection_extract(aa, "LINESTRING") st_collection_extract(aa, "POINT") ## With sfc objects st_collection_extract(st_geometry(aa), "POLYGON") st_collection_extract(st_geometry(aa), "LINESTRING") st_collection_extract(st_geometry(aa), "POINT") ## A GEOMETRY of single types bb <- rbind( st_sf(a = 1, geom = st_sfc(pt)), st_sf(a = 2, geom = st_sfc(ls)), st_sf(a = 3, geom = st_sfc(poly1)), st_sf(a = 4, geom = st_sfc(multipoly)) ) st_collection_extract(bb, "POLYGON") ## A GEOMETRY of mixed single types and GEOMETRYCOLLECTIONS cc <- rbind(aa, bb) st_collection_extract(cc, "POLYGON")
pt <- st_point(c(1, 0)) ls <- st_linestring(matrix(c(4, 3, 0, 0), ncol = 2)) poly1 <- st_polygon(list(matrix(c(5.5, 7, 7, 6, 5.5, 0, 0, -0.5, -0.5, 0), ncol = 2))) poly2 <- st_polygon(list(matrix(c(6.6, 8, 8, 7, 6.6, 1, 1, 1.5, 1.5, 1), ncol = 2))) multipoly <- st_multipolygon(list(poly1, poly2)) i <- st_geometrycollection(list(pt, ls, poly1, poly2)) j <- st_geometrycollection(list(pt, ls, poly1, poly2, multipoly)) st_collection_extract(i, "POLYGON") st_collection_extract(i, "POINT") st_collection_extract(i, "LINESTRING") ## A GEOMETRYCOLLECTION aa <- rbind(st_sf(a=1, geom = st_sfc(i)), st_sf(a=2, geom = st_sfc(j))) ## With sf objects st_collection_extract(aa, "POLYGON") st_collection_extract(aa, "LINESTRING") st_collection_extract(aa, "POINT") ## With sfc objects st_collection_extract(st_geometry(aa), "POLYGON") st_collection_extract(st_geometry(aa), "LINESTRING") st_collection_extract(st_geometry(aa), "POINT") ## A GEOMETRY of single types bb <- rbind( st_sf(a = 1, geom = st_sfc(pt)), st_sf(a = 2, geom = st_sfc(ls)), st_sf(a = 3, geom = st_sfc(poly1)), st_sf(a = 4, geom = st_sfc(multipoly)) ) st_collection_extract(bb, "POLYGON") ## A GEOMETRY of mixed single types and GEOMETRYCOLLECTIONS cc <- rbind(aa, bb) st_collection_extract(cc, "POLYGON")
retrieve coordinates in matrix form
st_coordinates(x, ...)
st_coordinates(x, ...)
x |
object of class sf, sfc or sfg |
... |
ignored |
matrix with coordinates (X, Y, possibly Z and/or M) in rows, possibly followed by integer indicators L1
,...,L3
that point out to which structure the coordinate belongs; for POINT
this is absent (each coordinate is a feature), for LINESTRING
L1
refers to the feature, for MULTILINESTRING
L1
refers to the part and L2
to the simple feature, for POLYGON
L1
refers to the main ring or holes and L2
to the simple feature, for MULTIPOLYGON
L1
refers to the main ring or holes, L2
to the ring id in the MULTIPOLYGON
, and L3
to the simple feature.
For POLYGONS
, L1
can be used to identify exterior rings and inner holes.
The exterior ring is when L1
is equal to 1. Interior rings are identified
when L1
is greater than 1. L2
can be used to differentiate between the
feature. Whereas for MULTIPOLYGON
, L3
refers to the MULTIPOLYGON
feature and L2
refers to the component POLYGON
.
crop an sf object to a specific rectangle
st_crop(x, y, ...) ## S3 method for class 'sfc' st_crop(x, y, ..., xmin, ymin, xmax, ymax) ## S3 method for class 'sf' st_crop(x, y, ...)
st_crop(x, y, ...) ## S3 method for class 'sfc' st_crop(x, y, ..., xmin, ymin, xmax, ymax) ## S3 method for class 'sf' st_crop(x, y, ...)
x |
object of class |
y |
numeric vector with named elements |
... |
ignored |
xmin |
minimum x extent of cropping area |
ymin |
minimum y extent of cropping area |
xmax |
maximum x extent of cropping area |
ymax |
maximum y extent of cropping area |
setting arguments xmin
, ymin
, xmax
and ymax
implies that argument y
gets ignored.
box = c(xmin = 0, ymin = 0, xmax = 1, ymax = 1) pol = st_sfc(st_buffer(st_point(c(.5, .5)), .6)) pol_sf = st_sf(a=1, geom=pol) plot(st_crop(pol, box)) plot(st_crop(pol_sf, st_bbox(box))) # alternative: plot(st_crop(pol, xmin = 0, ymin = 0, xmax = 1, ymax = 1))
box = c(xmin = 0, ymin = 0, xmax = 1, ymax = 1) pol = st_sfc(st_buffer(st_point(c(.5, .5)), .6)) pol_sf = st_sf(a=1, geom=pol) plot(st_crop(pol, box)) plot(st_crop(pol_sf, st_bbox(box))) # alternative: plot(st_crop(pol, xmin = 0, ymin = 0, xmax = 1, ymax = 1))
Retrieve coordinate reference system from sf or sfc object
Set or replace retrieve coordinate reference system from object
st_crs(x, ...) ## S3 method for class 'sf' st_crs(x, ...) ## S3 method for class 'numeric' st_crs(x, ...) ## S3 method for class 'character' st_crs(x, ...) ## S3 method for class 'sfc' st_crs(x, ..., parameters = FALSE) ## S3 method for class 'bbox' st_crs(x, ...) ## S3 method for class 'CRS' st_crs(x, ...) ## S3 method for class 'crs' st_crs(x, ...) st_crs(x) <- value ## S3 replacement method for class 'sf' st_crs(x) <- value ## S3 replacement method for class 'sfc' st_crs(x) <- value st_set_crs(x, value) NA_crs_ ## S3 method for class 'crs' is.na(x) ## S3 method for class 'crs' x$name ## S3 method for class 'crs' format(x, ...) st_axis_order(authority_compliant = logical(0))
st_crs(x, ...) ## S3 method for class 'sf' st_crs(x, ...) ## S3 method for class 'numeric' st_crs(x, ...) ## S3 method for class 'character' st_crs(x, ...) ## S3 method for class 'sfc' st_crs(x, ..., parameters = FALSE) ## S3 method for class 'bbox' st_crs(x, ...) ## S3 method for class 'CRS' st_crs(x, ...) ## S3 method for class 'crs' st_crs(x, ...) st_crs(x) <- value ## S3 replacement method for class 'sf' st_crs(x) <- value ## S3 replacement method for class 'sfc' st_crs(x) <- value st_set_crs(x, value) NA_crs_ ## S3 method for class 'crs' is.na(x) ## S3 method for class 'crs' x$name ## S3 method for class 'crs' format(x, ...) st_axis_order(authority_compliant = logical(0))
x |
|
... |
ignored |
parameters |
logical; |
value |
one of (i) character: a string accepted by GDAL, (ii) integer, a valid EPSG value (numeric), or (iii) an object of class |
name |
element name |
authority_compliant |
logical; specify whether axis order should be handled compliant to the authority; if omitted, the current value is printed. |
An object of class crs
of length 2.
The *crs functions create, get, set or replace the crs
attribute
of a simple feature geometry list-column. This attribute is of class crs
,
and is a list consisting of input
(user input, e.g. "EPSG:4326" or "WGS84"
or a proj4string), and wkt
, an automatically generated wkt2 representation of the crs.
If x
is identical to the wkt2 representation, and the CRS has a name, this name
is used for the input
field.
Comparison of two objects of class crs
uses the GDAL function
OGRSpatialReference::IsSame
.
In case a coordinate reference system is replaced, no transformation takes place and a warning is raised to stress this.
NA_crs_
is the crs
object with missing values for input
and wkt
.
the $
method for crs
objects retrieves named elements
using the GDAL interface; named elements include
SemiMajor
, SemiMinor
, InvFlattening
, IsGeographic
,
units_gdal
, IsVertical
, WktPretty
, Wkt
,
Name
, proj4string
, epsg
, yx
,
ud_unit
, and axes
(this may be subject to changes in future GDAL versions).
Note that not all valid CRS have a corresponding proj4string
.
ud_unit
returns a valid units object or NULL
if units are missing.
format.crs returns NA if the crs is missing valued, or else the name of a crs if it is different from "unknown", or else the user input if it was set, or else its "proj4string" representation;
st_axis_order
can be used to get and set the axis order: TRUE
indicates axes order according to the authority
(e.g. EPSG:4326 defining coordinates to be latitude,longitude pairs), FALSE
indicates the usual GIS (display) order (longitude,latitude). This can be useful
when data are read, or have to be written, with coordinates in authority compliant order.
The return value is the current state of this (FALSE
, by default).
If x
is numeric, return crs
object for EPSG:x
;
if x
is character, return crs
object for x
;
if x
is of class sf
or sfc
, return its crs
object.
Object of class crs
, which is a list with elements input
(length-1 character)
and wkt
(length-1 character).
Elements may be NA
valued; if all elements are NA
the CRS is missing valued, and coordinates are
assumed to relate to an arbitrary Cartesian coordinate system.
st_axis_order
returns the (logical) current value if called without
argument, or (invisibly) the previous value if it is being set.
sfc = st_sfc(st_point(c(0,0)), st_point(c(1,1))) sf = st_sf(a = 1:2, geom = sfc) st_crs(sf) = 4326 st_geometry(sf) sfc = st_sfc(st_point(c(0,0)), st_point(c(1,1))) st_crs(sfc) = 4326 sfc sfc = st_sfc(st_point(c(0,0)), st_point(c(1,1))) sfc %>% st_set_crs(4326) %>% st_transform(3857) st_crs("EPSG:3857")$input st_crs(3857)$proj4string pt = st_sfc(st_point(c(0, 60)), crs = 4326) # st_axis_order() only has effect in GDAL >= 2.5.0: st_axis_order() # query default: FALSE means interpret pt as (longitude latitude) st_transform(pt, 3857)[[1]] old_value = FALSE if (sf_extSoftVersion()["GDAL"] >= "2.5.0") (old_value = st_axis_order(TRUE)) # now interpret pt as (latitude longitude), as EPSG:4326 prescribes: st_axis_order() # query current value st_transform(pt, 3857)[[1]] st_axis_order(old_value) # set back to old value
sfc = st_sfc(st_point(c(0,0)), st_point(c(1,1))) sf = st_sf(a = 1:2, geom = sfc) st_crs(sf) = 4326 st_geometry(sf) sfc = st_sfc(st_point(c(0,0)), st_point(c(1,1))) st_crs(sfc) = 4326 sfc sfc = st_sfc(st_point(c(0,0)), st_point(c(1,1))) sfc %>% st_set_crs(4326) %>% st_transform(3857) st_crs("EPSG:3857")$input st_crs(3857)$proj4string pt = st_sfc(st_point(c(0, 60)), crs = 4326) # st_axis_order() only has effect in GDAL >= 2.5.0: st_axis_order() # query default: FALSE means interpret pt as (longitude latitude) st_transform(pt, 3857)[[1]] old_value = FALSE if (sf_extSoftVersion()["GDAL"] >= "2.5.0") (old_value = st_axis_order(TRUE)) # now interpret pt as (latitude longitude), as EPSG:4326 prescribes: st_axis_order() # query current value st_transform(pt, 3857)[[1]] st_axis_order(old_value) # set back to old value
Get a list of the available GDAL drivers
st_drivers(what = "vector", regex)
st_drivers(what = "vector", regex)
what |
character: |
regex |
character; regular expression to filter the |
The drivers available will depend on the installation of GDAL/OGR,
and can vary; the st_drivers()
function shows all the drivers that are
readable, and which may be written. The field vsi
refers to the driver's
capability to read/create datasets through the VSI*L API. See GDAL website for additional details on driver support
A data.frame
with driver metadata.
# The following driver lists depend on the GDAL setup and platform used: st_drivers() st_drivers("raster", "GeoT")
# The following driver lists depend on the GDAL setup and platform used: st_drivers() st_drivers("raster", "GeoT")
Get, set, replace or rename geometry from an sf object
## S3 method for class 'sfc' st_geometry(obj, ...) st_geometry(obj, ...) ## S3 method for class 'sf' st_geometry(obj, ...) ## S3 method for class 'sfc' st_geometry(obj, ...) ## S3 method for class 'sfg' st_geometry(obj, ...) st_geometry(x) <- value st_set_geometry(x, value) st_drop_geometry(x, ...) ## S3 method for class 'sf' st_drop_geometry(x, ...) ## Default S3 method: st_drop_geometry(x, ...)
## S3 method for class 'sfc' st_geometry(obj, ...) st_geometry(obj, ...) ## S3 method for class 'sf' st_geometry(obj, ...) ## S3 method for class 'sfc' st_geometry(obj, ...) ## S3 method for class 'sfg' st_geometry(obj, ...) st_geometry(x) <- value st_set_geometry(x, value) st_drop_geometry(x, ...) ## S3 method for class 'sf' st_drop_geometry(x, ...) ## Default S3 method: st_drop_geometry(x, ...)
obj |
object of class |
... |
ignored |
x |
object of class |
value |
object of class |
when applied to a data.frame
and when value
is an object of class sfc
, st_set_geometry
and st_geometry<-
will first check for the existence of an attribute sf_column
and overwrite that, or else look for list-columns of class sfc
and overwrite the first of that, or else write the geometry list-column to a column named geometry
. In case value
is character and x
is of class sf
, the "active" geometry column is set to x[[value]]
.
the replacement function applied to sf
objects will overwrite the geometry list-column, if value
is NULL
, it will remove it and coerce x
to a data.frame
.
if x
is of class sf
, st_drop_geometry
drops the geometry of its argument, and reclasses it accordingly; otherwise it returns x
unmodified.
st_geometry returns an object of class sfc, a list-column with geometries
st_geometry
returns an object of class sfc. Assigning geometry to a data.frame
creates an sf object, assigning it to an sf object replaces the geometry list-column.
df = data.frame(a = 1:2) sfc = st_sfc(st_point(c(3,4)), st_point(c(10,11))) st_geometry(sfc) st_geometry(df) <- sfc class(df) st_geometry(df) st_geometry(df) <- sfc # replaces st_geometry(df) <- NULL # remove geometry, coerce to data.frame sf <- st_set_geometry(df, sfc) # set geometry, return sf st_set_geometry(sf, NULL) # remove geometry, coerce to data.frame
df = data.frame(a = 1:2) sfc = st_sfc(st_point(c(3,4)), st_point(c(10,11))) st_geometry(sfc) st_geometry(df) <- sfc class(df) st_geometry(df) st_geometry(df) <- sfc # replaces st_geometry(df) <- NULL # remove geometry, coerce to data.frame sf <- st_set_geometry(df, sfc) # set geometry, return sf st_set_geometry(sf, NULL) # remove geometry, coerce to data.frame
Return geometry type of an object, as a factor
st_geometry_type(x, by_geometry = TRUE)
st_geometry_type(x, by_geometry = TRUE)
x |
|
by_geometry |
logical; if |
a factor with the geometry type of each simple feature geometry
in x
, or that of the whole set
Compute graticules and their parameters
st_graticule( x = c(-180, -90, 180, 90), crs = st_crs(x), datum = st_crs(4326), ..., lon = NULL, lat = NULL, ndiscr = 100, margin = 0.001 )
st_graticule( x = c(-180, -90, 180, 90), crs = st_crs(x), datum = st_crs(4326), ..., lon = NULL, lat = NULL, ndiscr = 100, margin = 0.001 )
x |
object of class |
crs |
object of class |
datum |
either an object of class |
... |
ignored |
lon |
numeric; values in degrees East for the meridians, associated with |
lat |
numeric; values in degrees North for the parallels, associated with |
ndiscr |
integer; number of points to discretize a parallel or meridian |
margin |
numeric; small number to trim a longlat bounding box that touches or crosses +/-180 long or +/-90 latitude. |
an object of class sf
with additional attributes describing the type
(E: meridian, N: parallel) degree value, label, start and end coordinates and angle;
see example.
In cartographic visualization, the use of graticules is not advised, unless the graphical output will be used for measurement or navigation, or the direction of North is important for the interpretation of the content, or the content is intended to display distortions and artifacts created by projection. Unnecessary use of graticules only adds visual clutter but little relevant information. Use of coastlines, administrative boundaries or place names permits most viewers of the output to orient themselves better than a graticule.
library(sf) if (require(maps, quietly = TRUE)) { usa = st_as_sf(map('usa', plot = FALSE, fill = TRUE)) laea = st_crs("+proj=laea +lat_0=30 +lon_0=-95") # Lambert equal area usa <- st_transform(usa, laea) bb = st_bbox(usa) bbox = st_linestring(rbind(c( bb[1],bb[2]),c( bb[3],bb[2]), c( bb[3],bb[4]),c( bb[1],bb[4]),c( bb[1],bb[2]))) g = st_graticule(usa) plot(usa, xlim = 1.2 * c(-2450853.4, 2186391.9), reset = FALSE) plot(g[1], add = TRUE, col = 'grey') plot(bbox, add = TRUE) points(g$x_start, g$y_start, col = 'red') points(g$x_end, g$y_end, col = 'blue') invisible(lapply(seq_len(nrow(g)), function(i) { if (g$type[i] == "N" && g$x_start[i] - min(g$x_start) < 1000) text(g$x_start[i], g$y_start[i], labels = parse(text = g$degree_label[i]), srt = g$angle_start[i], pos = 2, cex = .7) if (g$type[i] == "E" && g$y_start[i] - min(g$y_start) < 1000) text(g$x_start[i], g$y_start[i], labels = parse(text = g$degree_label[i]), srt = g$angle_start[i] - 90, pos = 1, cex = .7) if (g$type[i] == "N" && g$x_end[i] - max(g$x_end) > -1000) text(g$x_end[i], g$y_end[i], labels = parse(text = g$degree_label[i]), srt = g$angle_end[i], pos = 4, cex = .7) if (g$type[i] == "E" && g$y_end[i] - max(g$y_end) > -1000) text(g$x_end[i], g$y_end[i], labels = parse(text = g$degree_label[i]), srt = g$angle_end[i] - 90, pos = 3, cex = .7) })) plot(usa, graticule = st_crs(4326), axes = TRUE, lon = seq(-60,-130,by=-10)) }
library(sf) if (require(maps, quietly = TRUE)) { usa = st_as_sf(map('usa', plot = FALSE, fill = TRUE)) laea = st_crs("+proj=laea +lat_0=30 +lon_0=-95") # Lambert equal area usa <- st_transform(usa, laea) bb = st_bbox(usa) bbox = st_linestring(rbind(c( bb[1],bb[2]),c( bb[3],bb[2]), c( bb[3],bb[4]),c( bb[1],bb[4]),c( bb[1],bb[2]))) g = st_graticule(usa) plot(usa, xlim = 1.2 * c(-2450853.4, 2186391.9), reset = FALSE) plot(g[1], add = TRUE, col = 'grey') plot(bbox, add = TRUE) points(g$x_start, g$y_start, col = 'red') points(g$x_end, g$y_end, col = 'blue') invisible(lapply(seq_len(nrow(g)), function(i) { if (g$type[i] == "N" && g$x_start[i] - min(g$x_start) < 1000) text(g$x_start[i], g$y_start[i], labels = parse(text = g$degree_label[i]), srt = g$angle_start[i], pos = 2, cex = .7) if (g$type[i] == "E" && g$y_start[i] - min(g$y_start) < 1000) text(g$x_start[i], g$y_start[i], labels = parse(text = g$degree_label[i]), srt = g$angle_start[i] - 90, pos = 1, cex = .7) if (g$type[i] == "N" && g$x_end[i] - max(g$x_end) > -1000) text(g$x_end[i], g$y_end[i], labels = parse(text = g$degree_label[i]), srt = g$angle_end[i], pos = 4, cex = .7) if (g$type[i] == "E" && g$y_end[i] - max(g$y_end) > -1000) text(g$x_end[i], g$y_end[i], labels = parse(text = g$degree_label[i]), srt = g$angle_end[i] - 90, pos = 3, cex = .7) })) plot(usa, graticule = st_crs(4326), axes = TRUE, lon = seq(-60,-130,by=-10)) }
test equality between the geometry type and a class or set of classes
st_is(x, type)
st_is(x, type)
x |
object of class |
type |
character; class, or set of classes, to test against |
st_is(st_point(0:1), "POINT") sfc = st_sfc(st_point(0:1), st_linestring(matrix(1:6,,2))) st_is(sfc, "POINT") st_is(sfc, "POLYGON") st_is(sfc, "LINESTRING") st_is(st_sf(a = 1:2, sfc), "LINESTRING") st_is(sfc, c("POINT", "LINESTRING"))
st_is(st_point(0:1), "POINT") sfc = st_sfc(st_point(0:1), st_linestring(matrix(1:6,,2))) st_is(sfc, "POINT") st_is(sfc, "POLYGON") st_is(sfc, "LINESTRING") st_is(st_sf(a = 1:2, sfc), "LINESTRING") st_is(sfc, c("POINT", "LINESTRING"))
predicate whether a geometry is equal to a POLYGON FULL
st_is_full(x, ...) ## S3 method for class 'sfg' st_is_full(x, ..., is_longlat = NULL) ## S3 method for class 'sfc' st_is_full(x, ...) ## S3 method for class 'sf' st_is_full(x, ...) ## S3 method for class 'bbox' st_is_full(x, ...)
st_is_full(x, ...) ## S3 method for class 'sfg' st_is_full(x, ..., is_longlat = NULL) ## S3 method for class 'sfc' st_is_full(x, ...) ## S3 method for class 'sf' st_is_full(x, ...) ## S3 method for class 'bbox' st_is_full(x, ...)
x |
object of class |
... |
ignored, except when it contains a |
is_longlat |
logical; output of st_is_longlat of the parent |
logical, indicating whether geometries are POLYGON FULL (a spherical polygon covering the entire sphere)
Assert whether simple feature coordinates are longlat degrees
st_is_longlat(x)
st_is_longlat(x)
x |
object of class sf or sfc, or otherwise an object of a class that has an st_crs method returning a |
TRUE
if x
has geographic coordinates, FALSE
if it has projected coordinates, or NA
if is.na(st_crs(x))
.
jitter geometries
st_jitter(x, amount, factor = 0.002)
st_jitter(x, amount, factor = 0.002)
x |
object of class |
amount |
numeric; amount of jittering applied; if missing, the amount is set to factor * the bounding box diagonal; units of coordinates. |
factor |
numeric; fractional amount of jittering to be applied |
jitters coordinates with an amount such that runif(1, -amount, amount)
is added to the coordinates. x- and y-coordinates are jittered independently but all coordinates of a single geometry are jittered with the same amount, meaning that the geometry shape does not change. For longlat data, a latitude correction is made such that jittering in East and North directions are identical in distance in the center of the bounding box of x
.
nc = st_read(system.file("gpkg/nc.gpkg", package="sf")) pts = st_centroid(st_geometry(nc)) plot(pts) plot(st_jitter(pts, .05), add = TRUE, col = 'red') plot(st_geometry(nc)) plot(st_jitter(st_geometry(nc), factor = .01), add = TRUE, col = '#ff8888')
nc = st_read(system.file("gpkg/nc.gpkg", package="sf")) pts = st_centroid(st_geometry(nc)) plot(pts) plot(st_jitter(pts, .05), add = TRUE, col = 'red') plot(st_geometry(nc)) plot(st_jitter(st_geometry(nc), factor = .01), add = TRUE, col = '#ff8888')
spatial join, spatial filter
st_join(x, y, join, ...) ## S3 method for class 'sf' st_join( x, y, join = st_intersects, ..., suffix = c(".x", ".y"), left = TRUE, largest = FALSE ) st_filter(x, y, ...) ## S3 method for class 'sf' st_filter(x, y, ..., .predicate = st_intersects)
st_join(x, y, join, ...) ## S3 method for class 'sf' st_join( x, y, join = st_intersects, ..., suffix = c(".x", ".y"), left = TRUE, largest = FALSE ) st_filter(x, y, ...) ## S3 method for class 'sf' st_filter(x, y, ..., .predicate = st_intersects)
x |
object of class |
y |
object of class |
join |
geometry predicate function with the same profile as st_intersects; see details |
... |
for |
suffix |
length 2 character vector; see merge |
left |
logical; if |
largest |
logical; if |
.predicate |
geometry predicate function with the same profile as st_intersects; see details |
alternative values for argument join
are:
st_relate (which will require pattern
to be set),
or any user-defined function of the same profile as the above
A left join returns all records of the x
object with y
fields for non-matched records filled with NA
values; an inner join returns only records that spatially match.
To replicate the results of st_within(x, y)
you will need to use st_join(x, y, join = "st_within", left = FALSE)
.
an object of class sf
, joined based on geometry
a = st_sf(a = 1:3, geom = st_sfc(st_point(c(1,1)), st_point(c(2,2)), st_point(c(3,3)))) b = st_sf(a = 11:14, geom = st_sfc(st_point(c(10,10)), st_point(c(2,2)), st_point(c(2,2)), st_point(c(3,3)))) st_join(a, b) st_join(a, b, left = FALSE) # two ways to aggregate y's attribute values outcome over x's geometries: st_join(a, b) %>% aggregate(list(.$a.x), mean) if (require(dplyr, quietly = TRUE)) { st_join(a, b) %>% group_by(a.x) %>% summarise(mean(a.y)) } # example of largest = TRUE: nc <- st_transform(st_read(system.file("shape/nc.shp", package="sf")), 2264) gr = st_sf( label = apply(expand.grid(1:10, LETTERS[10:1])[,2:1], 1, paste0, collapse = " "), geom = st_make_grid(st_as_sfc(st_bbox(nc)))) gr$col = sf.colors(10, categorical = TRUE, alpha = .3) # cut, to check, NA's work out: gr = gr[-(1:30),] nc_j <- st_join(nc, gr, largest = TRUE) # the two datasets: opar = par(mfrow = c(2,1), mar = rep(0,4)) plot(st_geometry(nc_j)) plot(st_geometry(gr), add = TRUE, col = gr$col) text(st_coordinates(st_centroid(gr)), labels = gr$label) # the joined dataset: plot(st_geometry(nc_j), border = 'black', col = nc_j$col) text(st_coordinates(st_centroid(nc_j)), labels = nc_j$label, cex = .8) plot(st_geometry(gr), border = 'green', add = TRUE) par(opar) # st_filter keeps the geometries in x where .predicate(x,y) returns any match in y for x st_filter(a, b) # for an anti-join, use the union of y st_filter(a, st_union(b), .predicate = st_disjoint)
a = st_sf(a = 1:3, geom = st_sfc(st_point(c(1,1)), st_point(c(2,2)), st_point(c(3,3)))) b = st_sf(a = 11:14, geom = st_sfc(st_point(c(10,10)), st_point(c(2,2)), st_point(c(2,2)), st_point(c(3,3)))) st_join(a, b) st_join(a, b, left = FALSE) # two ways to aggregate y's attribute values outcome over x's geometries: st_join(a, b) %>% aggregate(list(.$a.x), mean) if (require(dplyr, quietly = TRUE)) { st_join(a, b) %>% group_by(a.x) %>% summarise(mean(a.y)) } # example of largest = TRUE: nc <- st_transform(st_read(system.file("shape/nc.shp", package="sf")), 2264) gr = st_sf( label = apply(expand.grid(1:10, LETTERS[10:1])[,2:1], 1, paste0, collapse = " "), geom = st_make_grid(st_as_sfc(st_bbox(nc)))) gr$col = sf.colors(10, categorical = TRUE, alpha = .3) # cut, to check, NA's work out: gr = gr[-(1:30),] nc_j <- st_join(nc, gr, largest = TRUE) # the two datasets: opar = par(mfrow = c(2,1), mar = rep(0,4)) plot(st_geometry(nc_j)) plot(st_geometry(gr), add = TRUE, col = gr$col) text(st_coordinates(st_centroid(gr)), labels = gr$label) # the joined dataset: plot(st_geometry(nc_j), border = 'black', col = nc_j$col) text(st_coordinates(st_centroid(nc_j)), labels = nc_j$label, cex = .8) plot(st_geometry(gr), border = 'green', add = TRUE) par(opar) # st_filter keeps the geometries in x where .predicate(x,y) returns any match in y for x st_filter(a, b) # for an anti-join, use the union of y st_filter(a, st_union(b), .predicate = st_disjoint)
Return properties of layers in a datasource
st_layers(dsn, options = character(0), do_count = FALSE)
st_layers(dsn, options = character(0), do_count = FALSE)
dsn |
data source name (interpretation varies by driver - for some drivers, |
options |
character; driver dependent dataset open options, multiple options supported. |
do_count |
logical; if TRUE, count the features by reading them, even if their count is not reported by the driver |
list object of class sf_layers
with elements
name of the layer
list with for each layer the geometry types
number of features (if reported; see do_count
)
number of fields
list with for each layer the crs
object
Project point on linestring, interpolate along a linestring
st_line_project(line, point, normalized = FALSE) st_line_interpolate(line, dist, normalized = FALSE)
st_line_project(line, point, normalized = FALSE) st_line_interpolate(line, dist, normalized = FALSE)
line |
object of class |
point |
object of class |
normalized |
logical; if |
dist |
numeric, vector with distance value(s) |
arguments line
, point
and dist
are recycled to common length when needed
st_line_project
returns the distance(s) of point(s) along line(s), when projected on the line(s)
st_line_interpolate
returns the point(s) at dist(s), when measured along (interpolated on) the line(s)
st_line_project(st_as_sfc("LINESTRING (0 0, 10 10)"), st_as_sfc(c("POINT (0 0)", "POINT (5 5)"))) st_line_project(st_as_sfc("LINESTRING (0 0, 10 10)"), st_as_sfc("POINT (5 5)"), TRUE) st_line_interpolate(st_as_sfc("LINESTRING (0 0, 1 1)"), 1) st_line_interpolate(st_as_sfc("LINESTRING (0 0, 1 1)"), 1, TRUE)
st_line_project(st_as_sfc("LINESTRING (0 0, 10 10)"), st_as_sfc(c("POINT (0 0)", "POINT (5 5)"))) st_line_project(st_as_sfc("LINESTRING (0 0, 10 10)"), st_as_sfc("POINT (5 5)"), TRUE) st_line_interpolate(st_as_sfc("LINESTRING (0 0, 1 1)"), 1) st_line_interpolate(st_as_sfc("LINESTRING (0 0, 1 1)"), 1, TRUE)
Sample points on a linear geometry
st_line_sample(x, n, density, type = "regular", sample = NULL)
st_line_sample(x, n, density, type = "regular", sample = NULL)
x |
object of class |
n |
integer; number of points to choose per geometry; if missing, n will be computed as |
density |
numeric; density (points per distance unit) of the sampling, possibly a vector of length equal to the number of features (otherwise recycled); |
type |
character; indicate the sampling type, either "regular" or "random" |
sample |
numeric; a vector of numbers between 0 and 1 indicating the points to sample - if defined sample overrules n, density and type. |
ls = st_sfc(st_linestring(rbind(c(0,0),c(0,1))), st_linestring(rbind(c(0,0),c(10,0)))) st_line_sample(ls, density = 1) ls = st_sfc(st_linestring(rbind(c(0,0),c(0,1))), st_linestring(rbind(c(0,0),c(.1,0))), crs = 4326) try(st_line_sample(ls, density = 1/1000)) # error st_line_sample(st_transform(ls, 3857), n = 5) # five points for each line st_line_sample(st_transform(ls, 3857), n = c(1, 3)) # one and three points st_line_sample(st_transform(ls, 3857), density = 1/1000) # one per km st_line_sample(st_transform(ls, 3857), density = c(1/1000, 1/10000)) # one per km, one per 10 km st_line_sample(st_transform(ls, 3857), density = units::set_units(1, 1/km)) # one per km # five equidistant points including start and end: st_line_sample(st_transform(ls, 3857), sample = c(0, 0.25, 0.5, 0.75, 1))
ls = st_sfc(st_linestring(rbind(c(0,0),c(0,1))), st_linestring(rbind(c(0,0),c(10,0)))) st_line_sample(ls, density = 1) ls = st_sfc(st_linestring(rbind(c(0,0),c(0,1))), st_linestring(rbind(c(0,0),c(.1,0))), crs = 4326) try(st_line_sample(ls, density = 1/1000)) # error st_line_sample(st_transform(ls, 3857), n = 5) # five points for each line st_line_sample(st_transform(ls, 3857), n = c(1, 3)) # one and three points st_line_sample(st_transform(ls, 3857), density = 1/1000) # one per km st_line_sample(st_transform(ls, 3857), density = c(1/1000, 1/10000)) # one per km, one per 10 km st_line_sample(st_transform(ls, 3857), density = units::set_units(1, 1/km)) # one per km # five equidistant points including start and end: st_line_sample(st_transform(ls, 3857), sample = c(0, 0.25, 0.5, 0.75, 1))
Return 'm' range of a simple feature or simple feature set
## S3 method for class 'm_range' is.na(x) st_m_range(obj, ...) ## S3 method for class 'POINT' st_m_range(obj, ...) ## S3 method for class 'MULTIPOINT' st_m_range(obj, ...) ## S3 method for class 'LINESTRING' st_m_range(obj, ...) ## S3 method for class 'POLYGON' st_m_range(obj, ...) ## S3 method for class 'MULTILINESTRING' st_m_range(obj, ...) ## S3 method for class 'MULTIPOLYGON' st_m_range(obj, ...) ## S3 method for class 'GEOMETRYCOLLECTION' st_m_range(obj, ...) ## S3 method for class 'MULTISURFACE' st_m_range(obj, ...) ## S3 method for class 'MULTICURVE' st_m_range(obj, ...) ## S3 method for class 'CURVEPOLYGON' st_m_range(obj, ...) ## S3 method for class 'COMPOUNDCURVE' st_m_range(obj, ...) ## S3 method for class 'POLYHEDRALSURFACE' st_m_range(obj, ...) ## S3 method for class 'TIN' st_m_range(obj, ...) ## S3 method for class 'TRIANGLE' st_m_range(obj, ...) ## S3 method for class 'CIRCULARSTRING' st_m_range(obj, ...) ## S3 method for class 'sfc' st_m_range(obj, ...) ## S3 method for class 'sf' st_m_range(obj, ...) ## S3 method for class 'numeric' st_m_range(obj, ..., crs = NA_crs_) NA_m_range_
## S3 method for class 'm_range' is.na(x) st_m_range(obj, ...) ## S3 method for class 'POINT' st_m_range(obj, ...) ## S3 method for class 'MULTIPOINT' st_m_range(obj, ...) ## S3 method for class 'LINESTRING' st_m_range(obj, ...) ## S3 method for class 'POLYGON' st_m_range(obj, ...) ## S3 method for class 'MULTILINESTRING' st_m_range(obj, ...) ## S3 method for class 'MULTIPOLYGON' st_m_range(obj, ...) ## S3 method for class 'GEOMETRYCOLLECTION' st_m_range(obj, ...) ## S3 method for class 'MULTISURFACE' st_m_range(obj, ...) ## S3 method for class 'MULTICURVE' st_m_range(obj, ...) ## S3 method for class 'CURVEPOLYGON' st_m_range(obj, ...) ## S3 method for class 'COMPOUNDCURVE' st_m_range(obj, ...) ## S3 method for class 'POLYHEDRALSURFACE' st_m_range(obj, ...) ## S3 method for class 'TIN' st_m_range(obj, ...) ## S3 method for class 'TRIANGLE' st_m_range(obj, ...) ## S3 method for class 'CIRCULARSTRING' st_m_range(obj, ...) ## S3 method for class 'sfc' st_m_range(obj, ...) ## S3 method for class 'sf' st_m_range(obj, ...) ## S3 method for class 'numeric' st_m_range(obj, ..., crs = NA_crs_) NA_m_range_
x |
object of class |
obj |
object to compute the m range from |
... |
ignored |
crs |
object of class |
An object of class m_range
of length 2.
NA_m_range_
represents the missing value for a m_range
object
a numeric vector of length two, with mmin
and mmax
values;
if obj
is of class sf
or sfc
the object
if obj
is of class sf
or sfc
the object
returned has a class m_range
a = st_sf(a = 1:2, geom = st_sfc(st_point(0:3), st_point(1:4)), crs = 4326) st_m_range(a) st_m_range(c(mmin = 16.1, mmax = 16.6), crs = st_crs(4326))
a = st_sf(a = 1:2, geom = st_sfc(st_point(0:3), st_point(1:4)), crs = 4326) st_m_range(a) st_m_range(c(mmin = 16.1, mmax = 16.6), crs = st_crs(4326))
Create a square or hexagonal grid covering the bounding box of the geometry of an sf or sfc object
st_make_grid( x, cellsize = c(diff(st_bbox(x)[c(1, 3)]), diff(st_bbox(x)[c(2, 4)]))/n, offset = st_bbox(x)[c("xmin", "ymin")], n = c(10, 10), crs = if (missing(x)) NA_crs_ else st_crs(x), what = "polygons", square = TRUE, flat_topped = FALSE )
st_make_grid( x, cellsize = c(diff(st_bbox(x)[c(1, 3)]), diff(st_bbox(x)[c(2, 4)]))/n, offset = st_bbox(x)[c("xmin", "ymin")], n = c(10, 10), crs = if (missing(x)) NA_crs_ else st_crs(x), what = "polygons", square = TRUE, flat_topped = FALSE )
x |
|
cellsize |
numeric of length 1 or 2 with target cellsize: for square or rectangular cells the width and height, for hexagonal cells the distance between opposite edges (edge length is cellsize/sqrt(3)). A length units object can be passed, or an area unit object with area size of the square or hexagonal cell. |
offset |
numeric of length 2; lower left corner coordinates (x, y) of the grid |
n |
integer of length 1 or 2, number of grid cells in x and y direction (columns, rows) |
crs |
object of class |
what |
character; one of: |
square |
logical; if |
flat_topped |
logical; if |
Object of class sfc
(simple feature geometry list column) with, depending on what
and square
,
square or hexagonal polygons, corner points of these polygons, or center points of these polygons.
plot(st_make_grid(what = "centers"), axes = TRUE) plot(st_make_grid(what = "corners"), add = TRUE, col = 'green', pch=3) sfc = st_sfc(st_polygon(list(rbind(c(0,0), c(1,0), c(1,1), c(0,0))))) plot(st_make_grid(sfc, cellsize = .1, square = FALSE)) plot(sfc, add = TRUE) # non-default offset: plot(st_make_grid(sfc, cellsize = .1, square = FALSE, offset = c(0, .05 / (sqrt(3)/2)))) plot(sfc, add = TRUE) nc = st_read(system.file("shape/nc.shp", package="sf")) g = st_make_grid(nc) plot(g) plot(st_geometry(nc), add = TRUE) # g[nc] selects cells that intersect with nc: plot(g[nc], col = '#ff000088', add = TRUE)
plot(st_make_grid(what = "centers"), axes = TRUE) plot(st_make_grid(what = "corners"), add = TRUE, col = 'green', pch=3) sfc = st_sfc(st_polygon(list(rbind(c(0,0), c(1,0), c(1,1), c(0,0))))) plot(st_make_grid(sfc, cellsize = .1, square = FALSE)) plot(sfc, add = TRUE) # non-default offset: plot(st_make_grid(sfc, cellsize = .1, square = FALSE, offset = c(0, .05 / (sqrt(3)/2)))) plot(sfc, add = TRUE) nc = st_read(system.file("shape/nc.shp", package="sf")) g = st_make_grid(nc) plot(g) plot(st_geometry(nc), add = TRUE) # g[nc] selects cells that intersect with nc: plot(g[nc], col = '#ff000088', add = TRUE)
get index of nearest feature
st_nearest_feature( x, y, ..., check_crs = TRUE, longlat = isTRUE(st_is_longlat(x)) )
st_nearest_feature( x, y, ..., check_crs = TRUE, longlat = isTRUE(st_is_longlat(x)) )
x |
object of class |
y |
object of class |
... |
ignored |
check_crs |
logical; should |
longlat |
logical; does |
for each feature (geometry) in x
the index of the nearest feature (geometry) in
set y
, or in the remaining set of x
if y
is missing;
empty geometries result in NA
indexes
st_nearest_points for finding the nearest points for pairs of feature geometries
ls1 = st_linestring(rbind(c(0,0), c(1,0))) ls2 = st_linestring(rbind(c(0,0.1), c(1,0.1))) ls3 = st_linestring(rbind(c(0,1), c(1,1))) (l = st_sfc(ls1, ls2, ls3)) p1 = st_point(c(0.1, -0.1)) p2 = st_point(c(0.1, 0.11)) p3 = st_point(c(0.1, 0.09)) p4 = st_point(c(0.1, 0.9)) (p = st_sfc(p1, p2, p3, p4)) try(st_nearest_feature(p, l)) try(st_nearest_points(p, l[st_nearest_feature(p,l)], pairwise = TRUE)) r = sqrt(2)/10 b1 = st_buffer(st_point(c(.1,.1)), r) b2 = st_buffer(st_point(c(.9,.9)), r) b3 = st_buffer(st_point(c(.9,.1)), r) circles = st_sfc(b1, b2, b3) plot(circles, col = NA, border = 2:4) pts = st_sfc(st_point(c(.3,.1)), st_point(c(.6,.2)), st_point(c(.6,.6)), st_point(c(.4,.8))) plot(pts, add = TRUE, col = 1) # draw points to nearest circle: nearest = try(st_nearest_feature(pts, circles)) if (inherits(nearest, "try-error")) # GEOS 3.6.1 not available nearest = c(1, 3, 2, 2) ls = st_nearest_points(pts, circles[nearest], pairwise = TRUE) plot(ls, col = 5:8, add = TRUE) # compute distance between pairs of nearest features: st_distance(pts, circles[nearest], by_element = TRUE)
ls1 = st_linestring(rbind(c(0,0), c(1,0))) ls2 = st_linestring(rbind(c(0,0.1), c(1,0.1))) ls3 = st_linestring(rbind(c(0,1), c(1,1))) (l = st_sfc(ls1, ls2, ls3)) p1 = st_point(c(0.1, -0.1)) p2 = st_point(c(0.1, 0.11)) p3 = st_point(c(0.1, 0.09)) p4 = st_point(c(0.1, 0.9)) (p = st_sfc(p1, p2, p3, p4)) try(st_nearest_feature(p, l)) try(st_nearest_points(p, l[st_nearest_feature(p,l)], pairwise = TRUE)) r = sqrt(2)/10 b1 = st_buffer(st_point(c(.1,.1)), r) b2 = st_buffer(st_point(c(.9,.9)), r) b3 = st_buffer(st_point(c(.9,.1)), r) circles = st_sfc(b1, b2, b3) plot(circles, col = NA, border = 2:4) pts = st_sfc(st_point(c(.3,.1)), st_point(c(.6,.2)), st_point(c(.6,.6)), st_point(c(.4,.8))) plot(pts, add = TRUE, col = 1) # draw points to nearest circle: nearest = try(st_nearest_feature(pts, circles)) if (inherits(nearest, "try-error")) # GEOS 3.6.1 not available nearest = c(1, 3, 2, 2) ls = st_nearest_points(pts, circles[nearest], pairwise = TRUE) plot(ls, col = 5:8, add = TRUE) # compute distance between pairs of nearest features: st_distance(pts, circles[nearest], by_element = TRUE)
get nearest points between pairs of geometries
st_nearest_points(x, y, ...) ## S3 method for class 'sfc' st_nearest_points(x, y, ..., pairwise = FALSE) ## S3 method for class 'sfg' st_nearest_points(x, y, ...) ## S3 method for class 'sf' st_nearest_points(x, y, ...)
st_nearest_points(x, y, ...) ## S3 method for class 'sfc' st_nearest_points(x, y, ..., pairwise = FALSE) ## S3 method for class 'sfg' st_nearest_points(x, y, ...) ## S3 method for class 'sf' st_nearest_points(x, y, ...)
x |
object of class |
y |
object of class |
... |
ignored |
pairwise |
logical; if |
in case x
lies inside y
, when using S2, the end points
are on polygon boundaries, when using GEOS the end point are identical to x
.
an sfc object with all two-point LINESTRING
geometries of point pairs from the first to the second geometry, of length x * y, with y cycling fastest. See examples for ideas how to convert these to POINT
geometries.
st_nearest_feature for finding the nearest feature
r = sqrt(2)/10 pt1 = st_point(c(.1,.1)) pt2 = st_point(c(.9,.9)) pt3 = st_point(c(.9,.1)) b1 = st_buffer(pt1, r) b2 = st_buffer(pt2, r) b3 = st_buffer(pt3, r) (ls0 = st_nearest_points(b1, b2)) # sfg (ls = st_nearest_points(st_sfc(b1), st_sfc(b2, b3))) # sfc plot(b1, xlim = c(-.2,1.2), ylim = c(-.2,1.2), col = NA, border = 'green') plot(st_sfc(b2, b3), add = TRUE, col = NA, border = 'blue') plot(ls, add = TRUE, col = 'red') nc = st_read(system.file("gpkg/nc.gpkg", package="sf")) plot(st_geometry(nc)) ls = st_nearest_points(nc[1,], nc) plot(ls, col = 'red', add = TRUE) pts = st_cast(ls, "POINT") # gives all start & end points # starting, "from" points, corresponding to x: plot(pts[seq(1, 200, 2)], add = TRUE, col = 'blue') # ending, "to" points, corresponding to y: plot(pts[seq(2, 200, 2)], add = TRUE, col = 'green')
r = sqrt(2)/10 pt1 = st_point(c(.1,.1)) pt2 = st_point(c(.9,.9)) pt3 = st_point(c(.9,.1)) b1 = st_buffer(pt1, r) b2 = st_buffer(pt2, r) b3 = st_buffer(pt3, r) (ls0 = st_nearest_points(b1, b2)) # sfg (ls = st_nearest_points(st_sfc(b1), st_sfc(b2, b3))) # sfc plot(b1, xlim = c(-.2,1.2), ylim = c(-.2,1.2), col = NA, border = 'green') plot(st_sfc(b2, b3), add = TRUE, col = NA, border = 'blue') plot(ls, add = TRUE, col = 'red') nc = st_read(system.file("gpkg/nc.gpkg", package="sf")) plot(st_geometry(nc)) ls = st_nearest_points(nc[1,], nc) plot(ls, col = 'red', add = TRUE) pts = st_cast(ls, "POINT") # gives all start & end points # starting, "from" points, corresponding to x: plot(pts[seq(1, 200, 2)], add = TRUE, col = 'blue') # ending, "to" points, corresponding to y: plot(pts[seq(2, 200, 2)], add = TRUE, col = 'green')
st_normalize
transforms the coordinates in the input feature to fall
between 0 and 1. By default the current domain is set to the bounding box of
the input, but other domains can be used as well
st_normalize(x, domain = st_bbox(x), ...)
st_normalize(x, domain = st_bbox(x), ...)
x |
object of class sf, sfc or sfg |
domain |
The domain |
... |
ignored |
p1 = st_point(c(7,52)) st_normalize(p1, domain = c(0, 0, 10, 100)) p2 = st_point(c(-30,20)) sfc = st_sfc(p1, p2, crs = 4326) sfc sfc_norm <- st_normalize(sfc) st_bbox(sfc_norm)
p1 = st_point(c(7,52)) st_normalize(p1, domain = c(0, 0, 10, 100)) p2 = st_point(c(-30,20)) sfc = st_sfc(p1, p2, crs = 4326) sfc sfc_norm <- st_normalize(sfc) st_bbox(sfc_norm)
Get precision
Set precision
st_precision(x) st_set_precision(x, precision) st_precision(x) <- value
st_precision(x) st_set_precision(x, precision) st_precision(x) <- value
x |
object of class |
precision |
numeric, or object of class |
value |
precision value |
If precision
is a units
object, the object on which we set precision must have a coordinate reference system with compatible distance units.
Setting a precision
has no direct effect on coordinates of geometries, but merely set an attribute tag to an sfc
object.
The effect takes place in st_as_binary or, more precise, in the C++ function CPL_write_wkb
, where simple feature geometries are being serialized to well-known-binary (WKB).
This happens always when routines are called in GEOS library (geometrical operations or predicates), for writing geometries using st_write or write_sf, st_make_valid
in package lwgeom
; also aggregate and summarise by default union geometries, which calls a GEOS library function.
Routines in these libraries receive rounded coordinates, and possibly return results based on them. st_as_binary contains an example of a roundtrip of sfc
geometries through WKB, in order to see the rounding happening to R data.
The reason to support precision is that geometrical operations in GEOS or liblwgeom may work better at reduced precision. For writing data from R to external resources it is harder to think of a good reason to limiting precision.
st_as_binary for an explanation of what setting precision does, and the examples therein.
x <- st_sfc(st_point(c(pi, pi))) st_precision(x) st_precision(x) <- 0.01 st_precision(x)
x <- st_sfc(st_point(c(pi, pi))) st_precision(x) st_precision(x) <- 0.01 st_precision(x)
Read simple features from file or database, or retrieve layer names and their geometry type(s)
Read PostGIS table directly through DBI and RPostgreSQL interface, converting Well-Know Binary geometries to sfc
st_read(dsn, layer, ...) ## S3 method for class 'character' st_read( dsn, layer, ..., query = NA, options = NULL, quiet = FALSE, geometry_column = 1L, type = 0, promote_to_multi = TRUE, stringsAsFactors = sf_stringsAsFactors(), int64_as_string = FALSE, check_ring_dir = FALSE, fid_column_name = character(0), drivers = character(0), wkt_filter = character(0), optional = FALSE, use_stream = default_st_read_use_stream() ) read_sf(..., quiet = TRUE, stringsAsFactors = FALSE, as_tibble = TRUE) ## S3 method for class 'DBIObject' st_read( dsn = NULL, layer = NULL, query = NULL, EWKB = TRUE, quiet = TRUE, as_tibble = FALSE, geometry_column = NULL, ... )
st_read(dsn, layer, ...) ## S3 method for class 'character' st_read( dsn, layer, ..., query = NA, options = NULL, quiet = FALSE, geometry_column = 1L, type = 0, promote_to_multi = TRUE, stringsAsFactors = sf_stringsAsFactors(), int64_as_string = FALSE, check_ring_dir = FALSE, fid_column_name = character(0), drivers = character(0), wkt_filter = character(0), optional = FALSE, use_stream = default_st_read_use_stream() ) read_sf(..., quiet = TRUE, stringsAsFactors = FALSE, as_tibble = TRUE) ## S3 method for class 'DBIObject' st_read( dsn = NULL, layer = NULL, query = NULL, EWKB = TRUE, quiet = TRUE, as_tibble = FALSE, geometry_column = NULL, ... )
dsn |
data source name (interpretation varies by driver - for some
drivers, |
layer |
layer name (varies by driver, may be a file name without
extension); in case |
... |
parameter(s) passed on to st_as_sf |
query |
SQL query to select records; see details |
options |
character; driver dependent dataset open options, multiple options supported. For possible values, see the "Open options" section of the GDAL documentation of the corresponding driver, and https://github.com/r-spatial/sf/issues/1157 for an example. |
quiet |
logical; suppress info on name, driver, size and spatial reference, or signaling no or multiple layers |
geometry_column |
integer or character; in case of multiple geometry fields, which one to take? |
type |
integer; ISO number of desired simple feature type; see details.
If left zero, and |
promote_to_multi |
logical; in case of a mix of Point and MultiPoint, or
of LineString and MultiLineString, or of Polygon and MultiPolygon, convert
all to the Multi variety; defaults to |
stringsAsFactors |
logical; logical: should character vectors be
converted to factors? Default for |
int64_as_string |
logical; if |
check_ring_dir |
logical; if |
fid_column_name |
character; name of column to write feature IDs to; defaults to not doing this |
drivers |
character; limited set of driver short names to be tried (default: try all) |
wkt_filter |
character; WKT representation of a spatial filter (may be used as bounding box, selecting overlapping geometries); see examples |
optional |
logical; passed to as.data.frame; always |
use_stream |
Use |
as_tibble |
logical; should the returned table be of class tibble or data.frame? |
EWKB |
logical; is the WKB of type EWKB? if missing, defaults to
|
for geometry_column
, see also
https://gdal.org/en/latest/development/rfc/rfc41_multiple_geometry_fields.html
for values for type
see
https://en.wikipedia.org/wiki/Well-known_text#Well-known_binary, but
note that not every target value may lead to successful conversion. The
typical conversion from POLYGON (3) to MULTIPOLYGON (6) should work; the
other way around (type=3), secondary rings from MULTIPOLYGONS may be dropped
without warnings. promote_to_multi
is handled on a per-geometry column
basis; type
may be specified for each geometry column.
Note that stray files in data source directories (such as *.dbf
) may
lead to spurious errors that accompanying *.shp
are missing.
In case of problems reading shapefiles from USB drives on OSX, please see
https://github.com/r-spatial/sf/issues/252. Reading shapefiles (or other
data sources) directly from zip files can be done by prepending the path
with /vsizip/
. This is part of the GDAL Virtual File Systems interface
that also supports .gz, curl, and other operations, including chaining; see
https://gdal.org/en/latest/user/virtual_file_systems.html for a complete
description and examples.
For query
with a character dsn
the query text is handed to
'ExecuteSQL' on the GDAL/OGR data set and will result in the creation of a
new layer (and layer
is ignored). See 'OGRSQL'
https://gdal.org/en/latest/user/ogr_sql_dialect.html for details. Please note that the
'FID' special field is driver-dependent, and may be either 0-based (e.g. ESRI
Shapefile), 1-based (e.g. MapInfo) or arbitrary (e.g. OSM). Other features of
OGRSQL are also likely to be driver dependent. The available layer names may
be obtained with
st_layers. Care will be required to properly escape the use of some layer names.
read_sf
and write_sf
are aliases for st_read
and st_write
, respectively, with some
modified default arguments.
read_sf
and write_sf
are quiet by default: they do not print information
about the data source. read_sf
returns an sf-tibble rather than an sf-data.frame.
write_sf
delete layers by default: it overwrites existing files without asking or warning.
if table
is not given but query
is, the spatial
reference system (crs) of the table queried is only available in case it
has been stored into each geometry record (e.g., by PostGIS, when using
EWKB)
The function will automatically find the geometry
type columns for
drivers that support it. For the other drivers, it will try to cast all the
character columns, which can be slow for very wide tables.
object of class sf when a layer was successfully read; in case
argument layer
is missing and data source dsn
does not
contain a single layer, an object of class sf_layers
is returned
with the layer names, each with their geometry type(s). Note that the
number of layers may also be zero.
The use of system.file
in examples make sure that examples run regardless where R is installed:
typical users will not use system.file
but give the file name directly, either with full path or relative
to the current working directory (see getwd). "Shapefiles" consist of several files with the same basename
that reside in the same directory, only one of them having extension .shp
.
nc = st_read(system.file("shape/nc.shp", package="sf")) summary(nc) # note that AREA was computed using Euclidian area on lon/lat degrees ## only three fields by select clause ## only two features by where clause nc_sql = st_read(system.file("shape/nc.shp", package="sf"), query = "SELECT NAME, SID74, FIPS FROM \"nc\" WHERE BIR74 > 20000") ## Not run: library(sp) example(meuse, ask = FALSE, echo = FALSE) try(st_write(st_as_sf(meuse), "PG:dbname=postgis", "meuse", layer_options = "OVERWRITE=true")) try(st_meuse <- st_read("PG:dbname=postgis", "meuse")) if (exists("st_meuse")) summary(st_meuse) ## End(Not run) ## Not run: ## note that we need special escaping of layer within single quotes (nc.gpkg) ## and that geom needs to be included in the select, otherwise we don't detect it layer <- st_layers(system.file("gpkg/nc.gpkg", package = "sf"))$name[1] nc_gpkg_sql = st_read(system.file("gpkg/nc.gpkg", package = "sf"), query = sprintf("SELECT NAME, SID74, FIPS, geom FROM \"%s\" WHERE BIR74 > 20000", layer)) ## End(Not run) # spatial filter, as wkt: wkt = st_as_text(st_geometry(nc[1,])) # filter by (bbox overlaps of) first feature geometry: st_read(system.file("gpkg/nc.gpkg", package="sf"), wkt_filter = wkt) # read geojson from string: geojson_txt <- paste("{\"type\":\"MultiPoint\",\"coordinates\":", "[[3.2,4],[3,4.6],[3.8,4.4],[3.5,3.8],[3.4,3.6],[3.9,4.5]]}") x = st_read(geojson_txt) x ## Not run: library(RPostgreSQL) try(conn <- dbConnect(PostgreSQL(), dbname = "postgis")) if (exists("conn") && !inherits(conn, "try-error")) { x = st_read(conn, "meuse", query = "select * from meuse limit 3;") x = st_read(conn, table = "public.meuse") print(st_crs(x)) # SRID resolved by the database, not by GDAL! dbDisconnect(conn) } ## End(Not run)
nc = st_read(system.file("shape/nc.shp", package="sf")) summary(nc) # note that AREA was computed using Euclidian area on lon/lat degrees ## only three fields by select clause ## only two features by where clause nc_sql = st_read(system.file("shape/nc.shp", package="sf"), query = "SELECT NAME, SID74, FIPS FROM \"nc\" WHERE BIR74 > 20000") ## Not run: library(sp) example(meuse, ask = FALSE, echo = FALSE) try(st_write(st_as_sf(meuse), "PG:dbname=postgis", "meuse", layer_options = "OVERWRITE=true")) try(st_meuse <- st_read("PG:dbname=postgis", "meuse")) if (exists("st_meuse")) summary(st_meuse) ## End(Not run) ## Not run: ## note that we need special escaping of layer within single quotes (nc.gpkg) ## and that geom needs to be included in the select, otherwise we don't detect it layer <- st_layers(system.file("gpkg/nc.gpkg", package = "sf"))$name[1] nc_gpkg_sql = st_read(system.file("gpkg/nc.gpkg", package = "sf"), query = sprintf("SELECT NAME, SID74, FIPS, geom FROM \"%s\" WHERE BIR74 > 20000", layer)) ## End(Not run) # spatial filter, as wkt: wkt = st_as_text(st_geometry(nc[1,])) # filter by (bbox overlaps of) first feature geometry: st_read(system.file("gpkg/nc.gpkg", package="sf"), wkt_filter = wkt) # read geojson from string: geojson_txt <- paste("{\"type\":\"MultiPoint\",\"coordinates\":", "[[3.2,4],[3,4.6],[3.8,4.4],[3.5,3.8],[3.4,3.6],[3.9,4.5]]}") x = st_read(geojson_txt) x ## Not run: library(RPostgreSQL) try(conn <- dbConnect(PostgreSQL(), dbname = "postgis")) if (exists("conn") && !inherits(conn, "try-error")) { x = st_read(conn, "meuse", query = "select * from meuse limit 3;") x = st_read(conn, table = "public.meuse") print(st_crs(x)) # SRID resolved by the database, not by GDAL! dbDisconnect(conn) } ## End(Not run)
Compute DE9-IM relation between pairs of geometries, or match it to a given pattern
st_relate(x, y, pattern = NA_character_, sparse = !is.na(pattern))
st_relate(x, y, pattern = NA_character_, sparse = !is.na(pattern))
x |
object of class |
y |
object of class |
pattern |
character; define the pattern to match to, see details. |
sparse |
logical; should a sparse matrix be returned ( |
In case pattern
is not given, st_relate
returns a dense character
matrix; element [i,j]
has nine characters, referring to the DE9-IM relationship between x[i]
and y[j]
, encoded as IxIy,IxBy,IxEy,BxIy,BxBy,BxEy,ExIy,ExBy,ExEy where I refers to interior, B to boundary, and E to exterior, and e.g. BxIy the dimensionality of the intersection of the the boundary of x[i]
and the interior of y[j]
, which is one of: 0, 1, 2, or F; digits denoting dimensionality of intersection, F denoting no intersection. When pattern
is given, a dense logical matrix or sparse index list returned with matches to the given pattern; see st_intersects for a description of the returned matrix or list. See also https://en.wikipedia.org/wiki/DE-9IM for further explanation.
p1 = st_point(c(0,0)) p2 = st_point(c(2,2)) pol1 = st_polygon(list(rbind(c(0,0),c(1,0),c(1,1),c(0,1),c(0,0)))) - 0.5 pol2 = pol1 + 1 pol3 = pol1 + 2 st_relate(st_sfc(p1, p2), st_sfc(pol1, pol2, pol3)) sfc = st_sfc(st_point(c(0,0)), st_point(c(3,3))) grd = st_make_grid(sfc, n = c(3,3)) st_intersects(grd) st_relate(grd, pattern = "****1****") # sides, not corners, internals st_relate(grd, pattern = "****0****") # only corners touch st_rook = function(a, b = a) st_relate(a, b, pattern = "F***1****") st_rook(grd) # queen neighbours, see \url{https://github.com/r-spatial/sf/issues/234#issuecomment-300511129} st_queen <- function(a, b = a) st_relate(a, b, pattern = "F***T****")
p1 = st_point(c(0,0)) p2 = st_point(c(2,2)) pol1 = st_polygon(list(rbind(c(0,0),c(1,0),c(1,1),c(0,1),c(0,0)))) - 0.5 pol2 = pol1 + 1 pol3 = pol1 + 2 st_relate(st_sfc(p1, p2), st_sfc(pol1, pol2, pol3)) sfc = st_sfc(st_point(c(0,0)), st_point(c(3,3))) grd = st_make_grid(sfc, n = c(3,3)) st_intersects(grd) st_relate(grd, pattern = "****1****") # sides, not corners, internals st_relate(grd, pattern = "****0****") # only corners touch st_rook = function(a, b = a) st_relate(a, b, pattern = "F***1****") st_rook(grd) # queen neighbours, see \url{https://github.com/r-spatial/sf/issues/234#issuecomment-300511129} st_queen <- function(a, b = a) st_relate(a, b, pattern = "F***T****")
Sample points on or in (sets of) spatial features.
By default, returns a pre-specified number of points that is equal to
size
(if type = "random"
and exact = TRUE
) or an approximation of
size
otherwise. spatstat
methods are
interfaced and do not use the size
argument, see examples.
st_sample(x, size, ...) ## S3 method for class 'sf' st_sample(x, size, ...) ## S3 method for class 'sfc' st_sample( x, size, ..., type = "random", exact = TRUE, warn_if_not_integer = TRUE, by_polygon = FALSE, progress = FALSE, force = FALSE ) ## S3 method for class 'sfg' st_sample(x, size, ...) ## S3 method for class 'bbox' st_sample( x, size, ..., great_circles = FALSE, segments = units::set_units(2, "degree", mode = "standard") )
st_sample(x, size, ...) ## S3 method for class 'sf' st_sample(x, size, ...) ## S3 method for class 'sfc' st_sample( x, size, ..., type = "random", exact = TRUE, warn_if_not_integer = TRUE, by_polygon = FALSE, progress = FALSE, force = FALSE ) ## S3 method for class 'sfg' st_sample(x, size, ...) ## S3 method for class 'bbox' st_sample( x, size, ..., great_circles = FALSE, segments = units::set_units(2, "degree", mode = "standard") )
x |
object of class |
size |
sample size(s) requested; either total size, or a numeric vector with sample sizes for each feature geometry. When sampling polygons, the returned sampling size may differ from the requested size, as the bounding box is sampled, and sampled points intersecting the polygon are returned. |
... |
passed on to sample for |
type |
character; indicates the spatial sampling type; one of |
exact |
logical; should the length of output be exactly |
warn_if_not_integer |
logical; if |
by_polygon |
logical; for |
progress |
logical; if |
force |
logical; if |
great_circles |
logical; if |
segments |
units, or numeric (degrees); segment sizes for segmenting a bounding box polygon if |
The function is vectorised: it samples size
points across all geometries in
the object if size
is a single number, or the specified number of points
in each feature if size
is a vector of integers equal in length to the geometry
of x
.
if x
has dimension 2 (polygons) and geographical coordinates (long/lat), uniform random sampling on the sphere is applied, see e.g. https://mathworld.wolfram.com/SpherePointPicking.html.
For regular
or hexagonal
sampling of polygons, the resulting size is only an approximation.
As parameter called offset
can be passed to control ("fix") regular or hexagonal sampling: for polygons a length 2 numeric vector (by default: a random point from st_bbox(x)
); for lines use a number like runif(1)
.
Fibonacci sampling see: Alvaro Gonzalez, 2010. Measurement of Areas on a Sphere Using Fibonacci and Latitude-Longitude Lattices. Mathematical Geosciences 42(1), p. 49-64
For regular sampling on the sphere, see also geosphere::regularCoordinates
.
Sampling methods from package spatstat
are interfaced (see examples), and need their own parameters to be set.
For instance, to use spatstat.random::rThomas()
, set type = "Thomas"
.
For sampling polygons one can specify oriented=TRUE
to make sure that polygons larger than half the globe are not reverted, e.g. when specifying a polygon from a bounding box of a global dataset. The st_sample
method for bbox
does this by default.
an sfc
object containing the sampled POINT
geometries
nc = st_read(system.file("shape/nc.shp", package="sf")) p1 = st_sample(nc[1:3, ], 6) p2 = st_sample(nc[1:3, ], 1:3) plot(st_geometry(nc)[1:3]) plot(p1, add = TRUE) plot(p2, add = TRUE, pch = 2) x = st_sfc(st_polygon(list(rbind(c(0,0),c(90,0),c(90,90),c(0,90),c(0,0)))), crs = st_crs(4326)) plot(x, axes = TRUE, graticule = TRUE) if (sf_extSoftVersion()["proj.4"] >= "4.9.0") plot(p <- st_sample(x, 1000), add = TRUE) if (require(lwgeom, quietly = TRUE)) { # for st_segmentize() x2 = st_transform(st_segmentize(x, 1e4), st_crs("+proj=ortho +lat_0=30 +lon_0=45")) g = st_transform(st_graticule(), st_crs("+proj=ortho +lat_0=30 +lon_0=45")) plot(x2, graticule = g) if (sf_extSoftVersion()["proj.4"] >= "4.9.0") { p2 = st_transform(p, st_crs("+proj=ortho +lat_0=30 +lon_0=45")) plot(p2, add = TRUE) } } x = st_sfc(st_polygon(list(rbind(c(0,0),c(90,0),c(90,10),c(0,90),c(0,0))))) # NOT long/lat: plot(x) p_exact = st_sample(x, 1000, exact = TRUE) p_not_exact = st_sample(x, 1000, exact = FALSE) length(p_exact); length(p_not_exact) plot(st_sample(x, 1000), add = TRUE) x = st_sfc(st_polygon(list(rbind(c(-180,-90),c(180,-90),c(180,90),c(-180,90),c(-180,-90)))), crs=st_crs(4326)) # FIXME: #if (sf_extSoftVersion()["proj.4"] >= "4.9.0") { # p = st_sample(x, 1000) # st_sample(p, 3) #} # hexagonal: sfc = st_sfc(st_polygon(list(rbind(c(0,0), c(1,0), c(1,1), c(0,0))))) plot(sfc) h = st_sample(sfc, 100, type = "hexagonal") h1 = st_sample(sfc, 100, type = "hexagonal") plot(h, add = TRUE) plot(h1, col = 'red', add = TRUE) c(length(h), length(h1)) # approximate! pt = st_multipoint(matrix(1:20,,2)) ls = st_sfc(st_linestring(rbind(c(0,0),c(0,1))), st_linestring(rbind(c(0,0),c(.1,0))), st_linestring(rbind(c(0,1),c(.1,1))), st_linestring(rbind(c(2,2),c(2,2.00001)))) st_sample(ls, 80) plot(st_sample(ls, 80)) # spatstat example: if (require(spatstat.random)) { x <- sf::st_sfc(sf::st_polygon(list(rbind(c(0, 0), c(10, 0), c(10, 10), c(0, 0))))) # for spatstat.random::rThomas(), set type = "Thomas": pts <- st_sample(x, kappa = 1, mu = 10, scale = 0.1, type = "Thomas") } bbox = st_bbox( c(xmin = 0, xmax = 40, ymax = 70, ymin = 60), crs = st_crs('OGC:CRS84') ) set.seed(13531) s1 = st_sample(bbox, 400) st_bbox(s1) # within bbox s2 = st_sample(bbox, 400, great_circles = TRUE) st_bbox(s2) # outside bbox
nc = st_read(system.file("shape/nc.shp", package="sf")) p1 = st_sample(nc[1:3, ], 6) p2 = st_sample(nc[1:3, ], 1:3) plot(st_geometry(nc)[1:3]) plot(p1, add = TRUE) plot(p2, add = TRUE, pch = 2) x = st_sfc(st_polygon(list(rbind(c(0,0),c(90,0),c(90,90),c(0,90),c(0,0)))), crs = st_crs(4326)) plot(x, axes = TRUE, graticule = TRUE) if (sf_extSoftVersion()["proj.4"] >= "4.9.0") plot(p <- st_sample(x, 1000), add = TRUE) if (require(lwgeom, quietly = TRUE)) { # for st_segmentize() x2 = st_transform(st_segmentize(x, 1e4), st_crs("+proj=ortho +lat_0=30 +lon_0=45")) g = st_transform(st_graticule(), st_crs("+proj=ortho +lat_0=30 +lon_0=45")) plot(x2, graticule = g) if (sf_extSoftVersion()["proj.4"] >= "4.9.0") { p2 = st_transform(p, st_crs("+proj=ortho +lat_0=30 +lon_0=45")) plot(p2, add = TRUE) } } x = st_sfc(st_polygon(list(rbind(c(0,0),c(90,0),c(90,10),c(0,90),c(0,0))))) # NOT long/lat: plot(x) p_exact = st_sample(x, 1000, exact = TRUE) p_not_exact = st_sample(x, 1000, exact = FALSE) length(p_exact); length(p_not_exact) plot(st_sample(x, 1000), add = TRUE) x = st_sfc(st_polygon(list(rbind(c(-180,-90),c(180,-90),c(180,90),c(-180,90),c(-180,-90)))), crs=st_crs(4326)) # FIXME: #if (sf_extSoftVersion()["proj.4"] >= "4.9.0") { # p = st_sample(x, 1000) # st_sample(p, 3) #} # hexagonal: sfc = st_sfc(st_polygon(list(rbind(c(0,0), c(1,0), c(1,1), c(0,0))))) plot(sfc) h = st_sample(sfc, 100, type = "hexagonal") h1 = st_sample(sfc, 100, type = "hexagonal") plot(h, add = TRUE) plot(h1, col = 'red', add = TRUE) c(length(h), length(h1)) # approximate! pt = st_multipoint(matrix(1:20,,2)) ls = st_sfc(st_linestring(rbind(c(0,0),c(0,1))), st_linestring(rbind(c(0,0),c(.1,0))), st_linestring(rbind(c(0,1),c(.1,1))), st_linestring(rbind(c(2,2),c(2,2.00001)))) st_sample(ls, 80) plot(st_sample(ls, 80)) # spatstat example: if (require(spatstat.random)) { x <- sf::st_sfc(sf::st_polygon(list(rbind(c(0, 0), c(10, 0), c(10, 10), c(0, 0))))) # for spatstat.random::rThomas(), set type = "Thomas": pts <- st_sample(x, kappa = 1, mu = 10, scale = 0.1, type = "Thomas") } bbox = st_bbox( c(xmin = 0, xmax = 40, ymax = 70, ymin = 60), crs = st_crs('OGC:CRS84') ) set.seed(13531) s1 = st_sample(bbox, 400) st_bbox(s1) # within bbox s2 = st_sample(bbox, 400, great_circles = TRUE) st_bbox(s2) # outside bbox
All longitudes < 0 are added to 360, to avoid for instance parts of Alaska
being represented on the far left and right of a plot because they have
values straddling 180 degrees. In general, using a projected
coordinate reference system is to be preferred, but this method permits a
geographical coordinate reference system to be used. This is the sf
equivalent of recenter in the sp package and
ST_ShiftLongitude
in PostGIS.
st_shift_longitude(x) ## S3 method for class 'sfc' st_shift_longitude(x, ...) ## S3 method for class 'sf' st_shift_longitude(x, ...)
st_shift_longitude(x) ## S3 method for class 'sfc' st_shift_longitude(x, ...) ## S3 method for class 'sf' st_shift_longitude(x, ...)
x |
object of class |
... |
ignored |
## sfc pt1 = st_point(c(-170, 50)) pt2 = st_point(c(170, 50)) (sfc = st_sfc(pt1, pt2)) sfc = st_set_crs(sfc, 4326) st_shift_longitude(sfc) ## sf d = st_as_sf(data.frame(id = 1:2, geometry = sfc)) st_shift_longitude(d)
## sfc pt1 = st_point(c(-170, 50)) pt2 = st_point(c(170, 50)) (sfc = st_sfc(pt1, pt2)) sfc = st_set_crs(sfc, 4326) st_shift_longitude(sfc) ## sf d = st_as_sf(data.frame(id = 1:2, geometry = sfc)) st_shift_longitude(d)
Transform or convert coordinates of simple feature
st_can_transform(src, dst) st_transform(x, crs, ...) ## S3 method for class 'sfc' st_transform( x, crs = st_crs(x), ..., aoi = numeric(0), pipeline = character(0), reverse = FALSE, desired_accuracy = -1, allow_ballpark = TRUE, partial = TRUE, check = FALSE ) ## S3 method for class 'sf' st_transform(x, crs = st_crs(x), ...) ## S3 method for class 'sfg' st_transform(x, crs = st_crs(x), ...) ## S3 method for class 'bbox' st_transform(x, crs, ..., densify = 21) st_wrap_dateline(x, options, quiet) ## S3 method for class 'sfc' st_wrap_dateline(x, options = "WRAPDATELINE=YES", quiet = TRUE) ## S3 method for class 'sf' st_wrap_dateline(x, options = "WRAPDATELINE=YES", quiet = TRUE) ## S3 method for class 'sfg' st_wrap_dateline(x, options = "WRAPDATELINE=YES", quiet = TRUE) sf_proj_info(type = "proj", path)
st_can_transform(src, dst) st_transform(x, crs, ...) ## S3 method for class 'sfc' st_transform( x, crs = st_crs(x), ..., aoi = numeric(0), pipeline = character(0), reverse = FALSE, desired_accuracy = -1, allow_ballpark = TRUE, partial = TRUE, check = FALSE ) ## S3 method for class 'sf' st_transform(x, crs = st_crs(x), ...) ## S3 method for class 'sfg' st_transform(x, crs = st_crs(x), ...) ## S3 method for class 'bbox' st_transform(x, crs, ..., densify = 21) st_wrap_dateline(x, options, quiet) ## S3 method for class 'sfc' st_wrap_dateline(x, options = "WRAPDATELINE=YES", quiet = TRUE) ## S3 method for class 'sf' st_wrap_dateline(x, options = "WRAPDATELINE=YES", quiet = TRUE) ## S3 method for class 'sfg' st_wrap_dateline(x, options = "WRAPDATELINE=YES", quiet = TRUE) sf_proj_info(type = "proj", path)
src |
source crs |
dst |
destination crs |
x |
object of class sf, sfc or sfg |
crs |
target coordinate reference system: object of class |
... |
ignored |
aoi |
area of interest, in degrees: WestLongitude, SouthLatitude, EastLongitude, NorthLatitude |
pipeline |
character; coordinate operation pipeline, for overriding the default operation |
reverse |
boolean; has only an effect when |
desired_accuracy |
numeric; Only coordinate operations that offer an accuracy of at least the one specified will be considered; a negative value disables this feature (requires GDAL >= 3.3) |
allow_ballpark |
logical; are ballpark (low accuracy) transformations allowed? (requires GDAL >= 3.3) |
partial |
logical; allow for partial projection, if not all points of a geometry can be projected (corresponds to setting environment variable |
check |
logical; if |
densify |
integer, number of points for discretizing lines between bounding box corner points; see Details |
options |
character; should have "WRAPDATELINE=YES" to function; another parameter that is used is "DATELINEOFFSET=10" (where 10 is the default value) |
quiet |
logical; print options after they have been parsed? |
type |
character; one of |
path |
character; PROJ search path to be set |
st_can_transform
returns a boolean indicating whether
coordinates with CRS src can be transformed into CRS dst
Transforms coordinates of object to new projection.
Features that cannot be transformed are returned as empty geometries.
Transforms using the pipeline=
argument may fail if there is
ambiguity in the axis order of the specified coordinate reference system;
if you need the traditional GIS order, use "OGC:CRS84"
, not
"EPSG:4326"
. Extra care is needed with the ESRI Shapefile format,
because WKT1 does not store axis order unambiguously.
The st_transform
method for sfg
objects assumes that the CRS of the object is available as an attribute of that name.
the method for bbox
objects densifies lines for geographic coordinates along Cartesian lines, not great circle arcs
For a discussion of using options
, see https://github.com/r-spatial/sf/issues/280 and https://github.com/r-spatial/sf/issues/1983
sf_proj_info
lists the available projections, ellipses, datums, units, or data search path of the PROJ library when type
is equal to proj, ellps, datum, units or path; when type
equals have_datum_files
a boolean is returned indicating whether datum files are installed and accessible (checking for conus
). path
returns the PROJ_INFO.searchpath
field directly, as a single string with path separaters (:
or ;
).
for PROJ >= 6, sf_proj_info
does not provide option type = "datums"
.
PROJ < 6 does not provide the option type = "prime_meridians"
.
for PROJ >= 7.1.0, the "units" query of sf_proj_info
returns the to_meter
variable as numeric, previous versions return a character vector containing a numeric expression.
st_transform_proj, part of package lwgeom.
sf_project projects a matrix of coordinates, bypassing GDAL altogether
p1 = st_point(c(7,52)) p2 = st_point(c(-30,20)) sfc = st_sfc(p1, p2, crs = 4326) sfc st_transform(sfc, 3857) st_transform(st_sf(a=2:1, geom=sfc), "EPSG:3857") if (sf_extSoftVersion()["GDAL"] >= "3.0.0") { st_transform(sfc, pipeline = "+proj=pipeline +step +proj=axisswap +order=2,1") # reverse axes st_transform(sfc, pipeline = "+proj=pipeline +step +proj=axisswap +order=2,1", reverse = TRUE) # also reverse axes } nc = st_read(system.file("shape/nc.shp", package="sf")) st_area(nc[1,]) # area from long/lat st_area(st_transform(nc[1,], 32119)) # NC state plane, m st_area(st_transform(nc[1,], 2264)) # NC state plane, US foot library(units) set_units(st_area(st_transform(nc[1,], 2264)), m^2) st_transform(structure(p1, proj4string = "EPSG:4326"), "EPSG:3857") st_wrap_dateline(st_sfc(st_linestring(rbind(c(-179,0),c(179,0))), crs = 4326)) sf_proj_info("datum")
p1 = st_point(c(7,52)) p2 = st_point(c(-30,20)) sfc = st_sfc(p1, p2, crs = 4326) sfc st_transform(sfc, 3857) st_transform(st_sf(a=2:1, geom=sfc), "EPSG:3857") if (sf_extSoftVersion()["GDAL"] >= "3.0.0") { st_transform(sfc, pipeline = "+proj=pipeline +step +proj=axisswap +order=2,1") # reverse axes st_transform(sfc, pipeline = "+proj=pipeline +step +proj=axisswap +order=2,1", reverse = TRUE) # also reverse axes } nc = st_read(system.file("shape/nc.shp", package="sf")) st_area(nc[1,]) # area from long/lat st_area(st_transform(nc[1,], 32119)) # NC state plane, m st_area(st_transform(nc[1,], 2264)) # NC state plane, US foot library(units) set_units(st_area(st_transform(nc[1,], 2264)), m^2) st_transform(structure(p1, proj4string = "EPSG:4326"), "EPSG:3857") st_wrap_dateline(st_sfc(st_linestring(rbind(c(-179,0),c(179,0))), crs = 4326)) sf_proj_info("datum")
Create viewport from sf, sfc or sfg object
st_viewport(x, ..., bbox = st_bbox(x), asp)
st_viewport(x, ..., bbox = st_bbox(x), asp)
x |
object of class sf, sfc or sfg object |
... |
parameters passed on to viewport |
bbox |
the bounding box used for aspect ratio |
asp |
numeric; target aspect ratio (y/x), see Details |
parameters width
, height
, xscale
and yscale
are set such that aspect ratio is honoured and plot size is maximized in the current viewport; others can be passed as ...
If asp
is missing, it is taken as 1, except when isTRUE(st_is_longlat(x))
, in which case it is set to 1.0 /cos(y)
, with y
the middle of the latitude bounding box.
The output of the call to viewport
library(grid) nc = st_read(system.file("shape/nc.shp", package="sf")) grid.newpage() pushViewport(viewport(width = 0.8, height = 0.8)) pushViewport(st_viewport(nc)) invisible(lapply(st_geometry(nc), function(x) grid.draw(st_as_grob(x, gp = gpar(fill = 'red')))))
library(grid) nc = st_read(system.file("shape/nc.shp", package="sf")) grid.newpage() pushViewport(viewport(width = 0.8, height = 0.8)) pushViewport(st_viewport(nc)) invisible(lapply(st_geometry(nc), function(x) grid.draw(st_as_grob(x, gp = gpar(fill = 'red')))))
Write simple features object to file or database
st_write(obj, dsn, layer, ...) ## S3 method for class 'sfc' st_write(obj, dsn, layer, ...) ## S3 method for class 'sf' st_write( obj, dsn, layer = NULL, ..., driver = guess_driver_can_write(dsn), dataset_options = NULL, layer_options = NULL, quiet = FALSE, factorsAsCharacter = TRUE, append = NA, delete_dsn = FALSE, delete_layer = !is.na(append) && !append, fid_column_name = NULL, config_options = character(0) ) ## S3 method for class 'data.frame' st_write(obj, dsn, layer = NULL, ...) write_sf(..., quiet = TRUE, append = FALSE, delete_layer = !append) st_delete( dsn, layer = character(0), driver = guess_driver_can_write(dsn), quiet = FALSE )
st_write(obj, dsn, layer, ...) ## S3 method for class 'sfc' st_write(obj, dsn, layer, ...) ## S3 method for class 'sf' st_write( obj, dsn, layer = NULL, ..., driver = guess_driver_can_write(dsn), dataset_options = NULL, layer_options = NULL, quiet = FALSE, factorsAsCharacter = TRUE, append = NA, delete_dsn = FALSE, delete_layer = !is.na(append) && !append, fid_column_name = NULL, config_options = character(0) ) ## S3 method for class 'data.frame' st_write(obj, dsn, layer = NULL, ...) write_sf(..., quiet = TRUE, append = FALSE, delete_layer = !append) st_delete( dsn, layer = character(0), driver = guess_driver_can_write(dsn), quiet = FALSE )
obj |
object of class |
dsn |
data source name. Interpretation varies by driver: can be
a filename, a folder, a database name, or a Database Connection
(we officially test support for
|
layer |
layer name. Varies by driver, may be a file name without
extension; for database connection, it is the name of the table. If layer
is missing, the |
... |
other arguments passed to dbWriteTable when |
driver |
character; name of driver to be used; if missing and |
dataset_options |
character; driver dependent dataset creation options; multiple options supported. |
layer_options |
character; driver dependent layer creation options; multiple options supported. |
quiet |
logical; suppress info on name, driver, size and spatial reference |
factorsAsCharacter |
logical; convert |
append |
logical; should we append to an existing layer, or replace it?
if |
delete_dsn |
logical; delete data source |
delete_layer |
logical; delete layer |
fid_column_name |
character, name of column with feature IDs; if specified, this column is no longer written as feature attribute. |
config_options |
character, named vector with GDAL config options |
Columns (variables) of a class not supported are dropped with a warning.
When updating an existing layer, records are appended to it if the updating object has the right variable names and types. If names don't match an error is raised. If types don't match, behaviour is undefined: GDAL may raise warnings or errors or fail silently.
When deleting layers or data sources is not successful, no error is emitted.
delete_dsn
and delete_layer
should be
handled with care; the former may erase complete directories or databases.
st_delete()
deletes layer(s) in a data source, or a data source if layers are
omitted; it returns TRUE
on success, FALSE
on failure, invisibly.
obj
, invisibly
nc = st_read(system.file("shape/nc.shp", package="sf")) st_write(nc, paste0(tempdir(), "/", "nc.shp")) st_write(nc, paste0(tempdir(), "/", "nc.shp"), delete_layer = TRUE) # overwrites if (require(sp, quietly = TRUE)) { data(meuse, package = "sp") # loads data.frame from sp meuse_sf = st_as_sf(meuse, coords = c("x", "y"), crs = 28992) # writes X and Y as columns: st_write(meuse_sf, paste0(tempdir(), "/", "meuse.csv"), layer_options = "GEOMETRY=AS_XY") st_write(meuse_sf, paste0(tempdir(), "/", "meuse.csv"), layer_options = "GEOMETRY=AS_WKT", delete_dsn=TRUE) # overwrites ## Not run: library(sp) example(meuse, ask = FALSE, echo = FALSE) try(st_write(st_as_sf(meuse), "PG:dbname=postgis", "meuse_sf", layer_options = c("OVERWRITE=yes", "LAUNDER=true"))) demo(nc, ask = FALSE) try(st_write(nc, "PG:dbname=postgis", "sids", layer_options = "OVERWRITE=true")) ## End(Not run) }
nc = st_read(system.file("shape/nc.shp", package="sf")) st_write(nc, paste0(tempdir(), "/", "nc.shp")) st_write(nc, paste0(tempdir(), "/", "nc.shp"), delete_layer = TRUE) # overwrites if (require(sp, quietly = TRUE)) { data(meuse, package = "sp") # loads data.frame from sp meuse_sf = st_as_sf(meuse, coords = c("x", "y"), crs = 28992) # writes X and Y as columns: st_write(meuse_sf, paste0(tempdir(), "/", "meuse.csv"), layer_options = "GEOMETRY=AS_XY") st_write(meuse_sf, paste0(tempdir(), "/", "meuse.csv"), layer_options = "GEOMETRY=AS_WKT", delete_dsn=TRUE) # overwrites ## Not run: library(sp) example(meuse, ask = FALSE, echo = FALSE) try(st_write(st_as_sf(meuse), "PG:dbname=postgis", "meuse_sf", layer_options = c("OVERWRITE=yes", "LAUNDER=true"))) demo(nc, ask = FALSE) try(st_write(nc, "PG:dbname=postgis", "sids", layer_options = "OVERWRITE=true")) ## End(Not run) }
Return 'z' range of a simple feature or simple feature set
## S3 method for class 'z_range' is.na(x) st_z_range(obj, ...) ## S3 method for class 'POINT' st_z_range(obj, ...) ## S3 method for class 'MULTIPOINT' st_z_range(obj, ...) ## S3 method for class 'LINESTRING' st_z_range(obj, ...) ## S3 method for class 'POLYGON' st_z_range(obj, ...) ## S3 method for class 'MULTILINESTRING' st_z_range(obj, ...) ## S3 method for class 'MULTIPOLYGON' st_z_range(obj, ...) ## S3 method for class 'GEOMETRYCOLLECTION' st_z_range(obj, ...) ## S3 method for class 'MULTISURFACE' st_z_range(obj, ...) ## S3 method for class 'MULTICURVE' st_z_range(obj, ...) ## S3 method for class 'CURVEPOLYGON' st_z_range(obj, ...) ## S3 method for class 'COMPOUNDCURVE' st_z_range(obj, ...) ## S3 method for class 'POLYHEDRALSURFACE' st_z_range(obj, ...) ## S3 method for class 'TIN' st_z_range(obj, ...) ## S3 method for class 'TRIANGLE' st_z_range(obj, ...) ## S3 method for class 'CIRCULARSTRING' st_z_range(obj, ...) ## S3 method for class 'sfc' st_z_range(obj, ...) ## S3 method for class 'sf' st_z_range(obj, ...) ## S3 method for class 'numeric' st_z_range(obj, ..., crs = NA_crs_) NA_z_range_
## S3 method for class 'z_range' is.na(x) st_z_range(obj, ...) ## S3 method for class 'POINT' st_z_range(obj, ...) ## S3 method for class 'MULTIPOINT' st_z_range(obj, ...) ## S3 method for class 'LINESTRING' st_z_range(obj, ...) ## S3 method for class 'POLYGON' st_z_range(obj, ...) ## S3 method for class 'MULTILINESTRING' st_z_range(obj, ...) ## S3 method for class 'MULTIPOLYGON' st_z_range(obj, ...) ## S3 method for class 'GEOMETRYCOLLECTION' st_z_range(obj, ...) ## S3 method for class 'MULTISURFACE' st_z_range(obj, ...) ## S3 method for class 'MULTICURVE' st_z_range(obj, ...) ## S3 method for class 'CURVEPOLYGON' st_z_range(obj, ...) ## S3 method for class 'COMPOUNDCURVE' st_z_range(obj, ...) ## S3 method for class 'POLYHEDRALSURFACE' st_z_range(obj, ...) ## S3 method for class 'TIN' st_z_range(obj, ...) ## S3 method for class 'TRIANGLE' st_z_range(obj, ...) ## S3 method for class 'CIRCULARSTRING' st_z_range(obj, ...) ## S3 method for class 'sfc' st_z_range(obj, ...) ## S3 method for class 'sf' st_z_range(obj, ...) ## S3 method for class 'numeric' st_z_range(obj, ..., crs = NA_crs_) NA_z_range_
x |
object of class |
obj |
object to compute the z range from |
... |
ignored |
crs |
object of class |
An object of class z_range
of length 2.
NA_z_range_
represents the missing value for a z_range
object
a numeric vector of length two, with zmin
and zmax
values;
if obj
is of class sf
or sfc
the object
returned has a class z_range
a = st_sf(a = 1:2, geom = st_sfc(st_point(0:2), st_point(1:3)), crs = 4326) st_z_range(a) st_z_range(c(zmin = 16.1, zmax = 16.6), crs = st_crs(4326))
a = st_sf(a = 1:2, geom = st_sfc(st_point(0:2), st_point(1:3)), crs = 4326) st_z_range(a) st_z_range(c(zmin = 16.1, zmax = 16.6), crs = st_crs(4326))
Drop Z and/or M dimensions from feature geometries, resetting classes appropriately
st_zm(x, ..., drop = TRUE, what = "ZM")
st_zm(x, ..., drop = TRUE, what = "ZM")
x |
object of class |
... |
ignored |
drop |
logical; drop, or ( |
what |
character which dimensions to drop or add |
Only combinations drop=TRUE
, what = "ZM"
, and drop=FALSE
, what="Z"
are supported so far.
In the latter case, x
should have XY
geometry, and zero values are added for the Z
dimension.
st_zm(st_linestring(matrix(1:32,8))) x = st_sfc(st_linestring(matrix(1:32,8)), st_linestring(matrix(1:8,2))) st_zm(x) a = st_sf(a = 1:2, geom=x) st_zm(a)
st_zm(st_linestring(matrix(1:32,8))) x = st_sfc(st_linestring(matrix(1:32,8)), st_linestring(matrix(1:8,2))) st_zm(x) a = st_sf(a = 1:2, geom=x) st_zm(a)
Summarize simple feature column
## S3 method for class 'sfc' summary(object, ..., maxsum = 7L, maxp4s = 10L)
## S3 method for class 'sfc' summary(object, ..., maxsum = 7L, maxp4s = 10L)
object |
object of class |
... |
ignored |
maxsum |
maximum number of classes to summarize the simple feature column to |
maxp4s |
maximum number of characters to print from the PROJ string |
Summarize simple feature type / item for tibble
type_sum.sfc(x, ...) obj_sum.sfc(x) pillar_shaft.sfc(x, ...)
type_sum.sfc(x, ...) obj_sum.sfc(x) pillar_shaft.sfc(x, ...)
x |
object of class |
... |
ignored |
see type_sum
Tidyverse methods for sf objects. Geometries are sticky, use as.data.frame to let dplyr
's own methods drop them.
Use these methods after loading the tidyverse package with the generic (or after loading package tidyverse).
filter.sf(.data, ..., .dots) arrange.sf(.data, ..., .dots) group_by.sf(.data, ..., add = FALSE) ungroup.sf(x, ...) rowwise.sf(x, ...) mutate.sf(.data, ..., .dots) transmute.sf(.data, ..., .dots) select.sf(.data, ...) rename.sf(.data, ...) rename_with.sf(.data, .fn, .cols, ...) slice.sf(.data, ..., .dots) summarise.sf(.data, ..., .dots, do_union = TRUE, is_coverage = FALSE) distinct.sf(.data, ..., .keep_all = FALSE, exact = FALSE, par = 0) gather.sf( data, key, value, ..., na.rm = FALSE, convert = FALSE, factor_key = FALSE ) pivot_longer.sf( data, cols, names_to = "name", names_prefix = NULL, names_sep = NULL, names_pattern = NULL, names_ptypes = NULL, names_transform = NULL, names_repair = "check_unique", values_to = "value", values_drop_na = FALSE, values_ptypes = NULL, values_transform = NULL, ... ) pivot_wider.sf( data, ..., id_cols = NULL, id_expand = FALSE, names_from = name, names_prefix = "", names_sep = "_", names_glue = NULL, names_sort = FALSE, names_vary = "fastest", names_expand = FALSE, names_repair = "check_unique", values_from = value, values_fill = NULL, values_fn = NULL, unused_fn = NULL ) spread.sf( data, key, value, fill = NA, convert = FALSE, drop = TRUE, sep = NULL ) sample_n.sf(tbl, size, replace = FALSE, weight = NULL, .env = parent.frame()) sample_frac.sf( tbl, size = 1, replace = FALSE, weight = NULL, .env = parent.frame() ) group_split.sf(.tbl, ..., .keep = TRUE) nest.sf(.data, ...) separate.sf( data, col, into, sep = "[^[:alnum:]]+", remove = TRUE, convert = FALSE, extra = "warn", fill = "warn", ... ) separate_rows.sf(data, ..., sep = "[^[:alnum:]]+", convert = FALSE) unite.sf(data, col, ..., sep = "_", remove = TRUE) unnest.sf(data, ..., .preserve = NULL) drop_na.sf(x, ...) inner_join.sf(x, y, by = NULL, copy = FALSE, suffix = c(".x", ".y"), ...) left_join.sf(x, y, by = NULL, copy = FALSE, suffix = c(".x", ".y"), ...) right_join.sf(x, y, by = NULL, copy = FALSE, suffix = c(".x", ".y"), ...) full_join.sf(x, y, by = NULL, copy = FALSE, suffix = c(".x", ".y"), ...) semi_join.sf(x, y, by = NULL, copy = FALSE, suffix = c(".x", ".y"), ...) anti_join.sf(x, y, by = NULL, copy = FALSE, suffix = c(".x", ".y"), ...)
filter.sf(.data, ..., .dots) arrange.sf(.data, ..., .dots) group_by.sf(.data, ..., add = FALSE) ungroup.sf(x, ...) rowwise.sf(x, ...) mutate.sf(.data, ..., .dots) transmute.sf(.data, ..., .dots) select.sf(.data, ...) rename.sf(.data, ...) rename_with.sf(.data, .fn, .cols, ...) slice.sf(.data, ..., .dots) summarise.sf(.data, ..., .dots, do_union = TRUE, is_coverage = FALSE) distinct.sf(.data, ..., .keep_all = FALSE, exact = FALSE, par = 0) gather.sf( data, key, value, ..., na.rm = FALSE, convert = FALSE, factor_key = FALSE ) pivot_longer.sf( data, cols, names_to = "name", names_prefix = NULL, names_sep = NULL, names_pattern = NULL, names_ptypes = NULL, names_transform = NULL, names_repair = "check_unique", values_to = "value", values_drop_na = FALSE, values_ptypes = NULL, values_transform = NULL, ... ) pivot_wider.sf( data, ..., id_cols = NULL, id_expand = FALSE, names_from = name, names_prefix = "", names_sep = "_", names_glue = NULL, names_sort = FALSE, names_vary = "fastest", names_expand = FALSE, names_repair = "check_unique", values_from = value, values_fill = NULL, values_fn = NULL, unused_fn = NULL ) spread.sf( data, key, value, fill = NA, convert = FALSE, drop = TRUE, sep = NULL ) sample_n.sf(tbl, size, replace = FALSE, weight = NULL, .env = parent.frame()) sample_frac.sf( tbl, size = 1, replace = FALSE, weight = NULL, .env = parent.frame() ) group_split.sf(.tbl, ..., .keep = TRUE) nest.sf(.data, ...) separate.sf( data, col, into, sep = "[^[:alnum:]]+", remove = TRUE, convert = FALSE, extra = "warn", fill = "warn", ... ) separate_rows.sf(data, ..., sep = "[^[:alnum:]]+", convert = FALSE) unite.sf(data, col, ..., sep = "_", remove = TRUE) unnest.sf(data, ..., .preserve = NULL) drop_na.sf(x, ...) inner_join.sf(x, y, by = NULL, copy = FALSE, suffix = c(".x", ".y"), ...) left_join.sf(x, y, by = NULL, copy = FALSE, suffix = c(".x", ".y"), ...) right_join.sf(x, y, by = NULL, copy = FALSE, suffix = c(".x", ".y"), ...) full_join.sf(x, y, by = NULL, copy = FALSE, suffix = c(".x", ".y"), ...) semi_join.sf(x, y, by = NULL, copy = FALSE, suffix = c(".x", ".y"), ...) anti_join.sf(x, y, by = NULL, copy = FALSE, suffix = c(".x", ".y"), ...)
.data |
data object of class sf |
... |
other arguments |
.dots |
see corresponding function in package |
add |
see corresponding function in dplyr |
x , y
|
A pair of data frames, data frame extensions (e.g. a tibble), or lazy data frames (e.g. from dbplyr or dtplyr). See Methods, below, for more details. |
.fn , .cols
|
see original docs |
do_union |
logical; in case |
is_coverage |
logical; if |
.keep_all |
see corresponding function in dplyr |
exact |
logical; if |
par |
numeric; passed on to st_equals_exact |
data |
see original function docs |
key |
see original function docs |
value |
see original function docs |
na.rm |
see original function docs |
convert |
see separate_rows |
factor_key |
see original function docs |
cols |
see original function docs |
names_to , names_pattern , names_ptypes , names_transform
|
|
names_prefix , names_sep , names_repair
|
see original function docs. |
values_to , values_drop_na , values_ptypes , values_transform
|
|
id_cols , id_expand , names_from , names_sort , names_glue , names_vary , names_expand
|
|
values_from , values_fill , values_fn , unused_fn
|
|
fill |
see original function docs |
drop |
see original function docs |
sep |
see separate_rows |
tbl |
see original function docs |
size |
see original function docs |
replace |
see original function docs |
weight |
see original function docs |
.env |
see original function docs |
.tbl |
see original function docs |
.keep |
see original function docs |
col |
see separate |
into |
see separate |
remove |
see separate |
extra |
see separate |
.preserve |
see unnest |
by |
A join specification created with If To join on different variables between To join by multiple variables, use a
For simple equality joins, you can alternatively specify a character vector
of variable names to join by. For example, To perform a cross-join, generating all combinations of |
copy |
If |
suffix |
If there are non-joined duplicate variables in |
select
keeps the geometry regardless whether it is selected or not; to deselect it, first pipe through as.data.frame
to let dplyr's own select
drop it.
In case one or more of the arguments (expressions) in the summarise
call creates a geometry list-column, the first of these will be the (active) geometry of the returned object. If this is not the case, a geometry column is created, depending on the value of do_union
.
In case do_union
is FALSE
, summarise
will simply combine geometries using c.sfg. When polygons sharing a boundary are combined, this leads to geometries that are invalid; see for instance https://github.com/r-spatial/sf/issues/681.
distinct
gives distinct records for which all attributes and geometries are distinct; st_equals is used to find out which geometries are distinct.
nest
assumes that a simple feature geometry list-column was among the columns that were nested.
an object of class sf
if (require(dplyr, quietly = TRUE)) { nc = read_sf(system.file("shape/nc.shp", package="sf")) nc %>% filter(AREA > .1) %>% plot() # plot 10 smallest counties in grey: st_geometry(nc) %>% plot() nc %>% select(AREA) %>% arrange(AREA) %>% slice(1:10) %>% plot(add = TRUE, col = 'grey') title("the ten counties with smallest area") nc2 <- nc %>% mutate(area10 = AREA/10) nc %>% slice(1:2) } # plot 10 smallest counties in grey: if (require(dplyr, quietly = TRUE)) { st_geometry(nc) %>% plot() nc %>% select(AREA) %>% arrange(AREA) %>% slice(1:10) %>% plot(add = TRUE, col = 'grey') title("the ten counties with smallest area") } if (require(dplyr, quietly = TRUE)) { nc$area_cl = cut(nc$AREA, c(0, .1, .12, .15, .25)) nc %>% group_by(area_cl) %>% class() } if (require(dplyr, quietly = TRUE)) { nc2 <- nc %>% mutate(area10 = AREA/10) } if (require(dplyr, quietly = TRUE)) { nc %>% transmute(AREA = AREA/10) %>% class() } if (require(dplyr, quietly = TRUE)) { nc %>% select(SID74, SID79) %>% names() nc %>% select(SID74, SID79) %>% class() } if (require(dplyr, quietly = TRUE)) { nc2 <- nc %>% rename(area = AREA) } if (require(dplyr, quietly = TRUE)) { nc %>% slice(1:2) } if (require(dplyr, quietly = TRUE)) { nc$area_cl = cut(nc$AREA, c(0, .1, .12, .15, .25)) nc.g <- nc %>% group_by(area_cl) nc.g %>% summarise(mean(AREA)) nc.g %>% summarise(mean(AREA)) %>% plot(col = grey(3:6 / 7)) nc %>% as.data.frame %>% summarise(mean(AREA)) } if (require(dplyr, quietly = TRUE)) { nc[c(1:100, 1:10), ] %>% distinct() %>% nrow() } if (require(tidyr, quietly = TRUE) && require(dplyr, quietly = TRUE) && "geometry" %in% names(nc)) { nc %>% select(SID74, SID79) %>% gather("VAR", "SID", -geometry) %>% summary() } if (require(tidyr, quietly = TRUE) && require(dplyr, quietly = TRUE) && "geometry" %in% names(nc)) { nc$row = 1:100 # needed for spread to work nc %>% select(SID74, SID79, geometry, row) %>% gather("VAR", "SID", -geometry, -row) %>% spread(VAR, SID) %>% head() } if (require(tidyr, quietly = TRUE) && require(dplyr, quietly = TRUE)) { storms.sf = st_as_sf(storms, coords = c("long", "lat"), crs = 4326) x <- storms.sf %>% group_by(name, year) %>% nest trs = lapply(x$data, function(tr) st_cast(st_combine(tr), "LINESTRING")[[1]]) %>% st_sfc(crs = 4326) trs.sf = st_sf(x[,1:2], trs) plot(trs.sf["year"], axes = TRUE) }
if (require(dplyr, quietly = TRUE)) { nc = read_sf(system.file("shape/nc.shp", package="sf")) nc %>% filter(AREA > .1) %>% plot() # plot 10 smallest counties in grey: st_geometry(nc) %>% plot() nc %>% select(AREA) %>% arrange(AREA) %>% slice(1:10) %>% plot(add = TRUE, col = 'grey') title("the ten counties with smallest area") nc2 <- nc %>% mutate(area10 = AREA/10) nc %>% slice(1:2) } # plot 10 smallest counties in grey: if (require(dplyr, quietly = TRUE)) { st_geometry(nc) %>% plot() nc %>% select(AREA) %>% arrange(AREA) %>% slice(1:10) %>% plot(add = TRUE, col = 'grey') title("the ten counties with smallest area") } if (require(dplyr, quietly = TRUE)) { nc$area_cl = cut(nc$AREA, c(0, .1, .12, .15, .25)) nc %>% group_by(area_cl) %>% class() } if (require(dplyr, quietly = TRUE)) { nc2 <- nc %>% mutate(area10 = AREA/10) } if (require(dplyr, quietly = TRUE)) { nc %>% transmute(AREA = AREA/10) %>% class() } if (require(dplyr, quietly = TRUE)) { nc %>% select(SID74, SID79) %>% names() nc %>% select(SID74, SID79) %>% class() } if (require(dplyr, quietly = TRUE)) { nc2 <- nc %>% rename(area = AREA) } if (require(dplyr, quietly = TRUE)) { nc %>% slice(1:2) } if (require(dplyr, quietly = TRUE)) { nc$area_cl = cut(nc$AREA, c(0, .1, .12, .15, .25)) nc.g <- nc %>% group_by(area_cl) nc.g %>% summarise(mean(AREA)) nc.g %>% summarise(mean(AREA)) %>% plot(col = grey(3:6 / 7)) nc %>% as.data.frame %>% summarise(mean(AREA)) } if (require(dplyr, quietly = TRUE)) { nc[c(1:100, 1:10), ] %>% distinct() %>% nrow() } if (require(tidyr, quietly = TRUE) && require(dplyr, quietly = TRUE) && "geometry" %in% names(nc)) { nc %>% select(SID74, SID79) %>% gather("VAR", "SID", -geometry) %>% summary() } if (require(tidyr, quietly = TRUE) && require(dplyr, quietly = TRUE) && "geometry" %in% names(nc)) { nc$row = 1:100 # needed for spread to work nc %>% select(SID74, SID79, geometry, row) %>% gather("VAR", "SID", -geometry, -row) %>% spread(VAR, SID) %>% head() } if (require(tidyr, quietly = TRUE) && require(dplyr, quietly = TRUE)) { storms.sf = st_as_sf(storms, coords = c("long", "lat"), crs = 4326) x <- storms.sf %>% group_by(name, year) %>% nest trs = lapply(x$data, function(tr) st_cast(st_combine(tr), "LINESTRING")[[1]]) %>% st_sfc(crs = 4326) trs.sf = st_sf(x[,1:2], trs) plot(trs.sf["year"], axes = TRUE) }
Can be used to create or modify attribute variables; for transforming geometries see
st_transform, and all other functions starting with st_
.
## S3 method for class 'sf' transform(`_data`, ...)
## S3 method for class 'sf' transform(`_data`, ...)
_data |
object of class |
... |
Further arguments of the form |
a = data.frame(x1 = 1:3, x2 = 5:7) st_geometry(a) = st_sfc(st_point(c(0,0)), st_point(c(1,1)), st_point(c(2,2))) transform(a, x1_sq = x1^2) transform(a, x1_x2 = x1*x2)
a = data.frame(x1 = 1:3, x2 = 5:7) st_geometry(a) = st_sfc(st_point(c(0,0)), st_point(c(1,1)), st_point(c(2,2))) transform(a, x1_sq = x1^2) transform(a, x1_x2 = x1*x2)
Checks whether a geometry is valid, or makes an invalid geometry valid
st_is_valid(x, ...) ## S3 method for class 'sfc' st_is_valid(x, ..., NA_on_exception = TRUE, reason = FALSE) ## S3 method for class 'sf' st_is_valid(x, ...) ## S3 method for class 'sfg' st_is_valid(x, ...) st_make_valid(x, ...) ## S3 method for class 'sfg' st_make_valid(x, ...) ## S3 method for class 'sfc' st_make_valid( x, ..., oriented = FALSE, s2_options = s2::s2_options(snap = s2::s2_snap_precision(1e+07), ...), geos_method = "valid_structure", geos_keep_collapsed = TRUE )
st_is_valid(x, ...) ## S3 method for class 'sfc' st_is_valid(x, ..., NA_on_exception = TRUE, reason = FALSE) ## S3 method for class 'sf' st_is_valid(x, ...) ## S3 method for class 'sfg' st_is_valid(x, ...) st_make_valid(x, ...) ## S3 method for class 'sfg' st_make_valid(x, ...) ## S3 method for class 'sfc' st_make_valid( x, ..., oriented = FALSE, s2_options = s2::s2_options(snap = s2::s2_snap_precision(1e+07), ...), geos_method = "valid_structure", geos_keep_collapsed = TRUE )
x |
object of class |
... |
passed on to s2_options |
NA_on_exception |
logical; if TRUE, for polygons that would otherwise raise a GEOS error (exception, e.g. for a POLYGON having more than zero but less than 4 points, or a LINESTRING having one point) return an |
reason |
logical; if |
oriented |
logical; only relevant if |
s2_options |
only relevant if |
geos_method |
character; either "valid_linework" (Original method, combines all rings into a set of noded lines and then extracts valid polygons from that linework) or "valid_structure" (Structured method, first makes all rings valid then merges shells and subtracts holes from shells to generate valid result. Assumes that holes and shells are correctly categorized.) (requires GEOS >= 3.10.1) |
geos_keep_collapsed |
logical; When this parameter is not set to |
For projected geometries, st_make_valid
uses the lwgeom_makevalid
method also used by the PostGIS command ST_makevalid
if the GEOS version linked to is smaller than 3.8.0, and otherwise the version shipped in GEOS; for geometries having ellipsoidal coordinates s2::s2_rebuild
is being used.
if s2_options
is not specified and x
has a non-zero precision set, then this precision value will be used as the value in s2_snap_precision
, passed on to s2_options
, rather than the 1e7 default.
st_is_valid
returns a logical vector indicating for each geometries of x
whether it is valid. st_make_valid
returns an object with a topologically valid geometry.
Object of the same class as x
p1 = st_as_sfc("POLYGON((0 0, 0 10, 10 0, 10 10, 0 0))") st_is_valid(p1) st_is_valid(st_sfc(st_point(0:1), p1[[1]]), reason = TRUE) library(sf) x = st_sfc(st_polygon(list(rbind(c(0,0),c(0.5,0),c(0.5,0.5),c(0.5,0),c(1,0),c(1,1),c(0,1),c(0,0))))) suppressWarnings(st_is_valid(x)) y = st_make_valid(x) st_is_valid(y) y %>% st_cast()
p1 = st_as_sfc("POLYGON((0 0, 0 10, 10 0, 10 10, 0 0))") st_is_valid(p1) st_is_valid(st_sfc(st_point(0:1), p1[[1]]), reason = TRUE) library(sf) x = st_sfc(st_polygon(list(rbind(c(0,0),c(0.5,0),c(0.5,0.5),c(0.5,0),c(1,0),c(1,1),c(0,1),c(0,0))))) suppressWarnings(st_is_valid(x)) y = st_make_valid(x) st_is_valid(y) y %>% st_cast()
vctrs methods for sf objects
vec_ptype2.sfc(x, y, ...) ## Default S3 method: vec_ptype2.sfc(x, y, ..., x_arg = "x", y_arg = "y") ## S3 method for class 'sfc' vec_ptype2.sfc(x, y, ...) vec_cast.sfc(x, to, ...) ## S3 method for class 'sfc' vec_cast.sfc(x, to, ...) ## Default S3 method: vec_cast.sfc(x, to, ...)
vec_ptype2.sfc(x, y, ...) ## Default S3 method: vec_ptype2.sfc(x, y, ..., x_arg = "x", y_arg = "y") ## S3 method for class 'sfc' vec_ptype2.sfc(x, y, ...) vec_cast.sfc(x, to, ...) ## S3 method for class 'sfc' vec_cast.sfc(x, to, ...) ## Default S3 method: vec_cast.sfc(x, to, ...)
x , y
|
Vector types. |
... |
These dots are for future extensions and must be empty. |
x_arg , y_arg
|
Argument names for |
to |
Type to cast to. If |