R Under development (unstable) (2024-11-03 r87286 ucrt) -- "Unsuffered Consequences" Copyright (C) 2024 The R Foundation for Statistical Computing Platform: x86_64-w64-mingw32/x64 R is free software and comes with ABSOLUTELY NO WARRANTY. You are welcome to redistribute it under certain conditions. Type 'license()' or 'licence()' for distribution details. R is a collaborative project with many contributors. Type 'contributors()' for more information and 'citation()' on how to cite R or R packages in publications. Type 'demo()' for some demos, 'help()' for on-line help, or 'help.start()' for an HTML browser interface to help. Type 'q()' to quit R. > suppressPackageStartupMessages(library(sf)) > > p = st_point(c(1/3,1/6)) > st_sfc(p, precision = 1000) Geometry set for 1 feature Geometry type: POINT Dimension: XY Bounding box: xmin: 0.3333333 ymin: 0.1666667 xmax: 0.3333333 ymax: 0.1666667 CRS: NA Precision: 1000 POINT (0.3333333 0.1666667) > st_as_sfc(st_as_binary(st_sfc(p, precision = 0L))) Geometry set for 1 feature Geometry type: POINT Dimension: XY Bounding box: xmin: 0.3333333 ymin: 0.1666667 xmax: 0.3333333 ymax: 0.1666667 CRS: NA POINT (0.3333333 0.1666667) > st_as_sfc(st_as_binary(st_sfc(p, precision = 1000))) Geometry set for 1 feature Geometry type: POINT Dimension: XY Bounding box: xmin: 0.333 ymin: 0.167 xmax: 0.333 ymax: 0.167 CRS: NA POINT (0.333 0.167) > st_as_sfc(st_as_binary(st_sfc(p, precision = 1000000))) Geometry set for 1 feature Geometry type: POINT Dimension: XY Bounding box: xmin: 0.333333 ymin: 0.166667 xmax: 0.333333 ymax: 0.166667 CRS: NA POINT (0.333333 0.166667) > st_as_sfc(st_as_binary(st_sfc(p, precision = 10L))) Geometry set for 1 feature Geometry type: POINT Dimension: XY Bounding box: xmin: 0.3 ymin: 0.2 xmax: 0.3 ymax: 0.2 CRS: NA POINT (0.3 0.2) > st_as_sfc(st_as_binary(st_sfc(p, precision = -1))) Geometry set for 1 feature Geometry type: POINT Dimension: XY Bounding box: xmin: 0.3333333 ymin: 0.1666667 xmax: 0.3333333 ymax: 0.1666667 CRS: NA POINT (0.3333333 0.1666667) > > d = data.frame(a = 1:2) > d$geom = c("POINT(0 0)", "POINT(1 1)") > > st_as_sf(d, wkt = "geom") Simple feature collection with 2 features and 1 field Geometry type: POINT Dimension: XY Bounding box: xmin: 0 ymin: 0 xmax: 1 ymax: 1 CRS: NA a geom 1 1 POINT (0 0) 2 2 POINT (1 1) > st_as_sf(d, wkt = 2) Simple feature collection with 2 features and 1 field Geometry type: POINT Dimension: XY Bounding box: xmin: 0 ymin: 0 xmax: 1 ymax: 1 CRS: NA a geom 1 1 POINT (0 0) 2 2 POINT (1 1) > st_as_sf(d, wkt = "geom", remove = FALSE) Simple feature collection with 2 features and 2 fields Geometry type: POINT Dimension: XY Bounding box: xmin: 0 ymin: 0 xmax: 1 ymax: 1 CRS: NA a geom geometry 1 1 POINT(0 0) POINT (0 0) 2 2 POINT(1 1) POINT (1 1) > > st_as_sfc(c("POINT(0 0)", "POINT(1 1)")) Geometry set for 2 features Geometry type: POINT Dimension: XY Bounding box: xmin: 0 ymin: 0 xmax: 1 ymax: 1 CRS: NA POINT (0 0) POINT (1 1) > st_as_sfc(c("POINT(0 0)", "POINT(1 1)", "POLYGON((0 0,1 1,0 1,0 0))")) Geometry set for 3 features Geometry type: GEOMETRY Dimension: XY Bounding box: xmin: 0 ymin: 0 xmax: 1 ymax: 1 CRS: NA POINT (0 0) POINT (1 1) POLYGON ((0 0, 1 1, 0 1, 0 0)) > st_as_sfc(character(0)) Geometry set for 0 features Bounding box: xmin: NA ymin: NA xmax: NA ymax: NA CRS: NA > x = st_as_sfc(character(0), 4326) > y = st_as_sfc(character(0), crs = 4326) > all.equal(x, y) [1] TRUE > st_as_sfc(c("POINT(0 0)", "POINT(1 1)", "POLYGON((0 0,1 1,0 1,0 0))"), + "+proj=longlat +datum=WGS84") Geometry set for 3 features Geometry type: GEOMETRY Dimension: XY Bounding box: xmin: 0 ymin: 0 xmax: 1 ymax: 1 Geodetic CRS: +proj=longlat +datum=WGS84 POINT (0 0) POINT (1 1) POLYGON ((0 0, 1 1, 0 1, 0 0)) > dg = st_as_sf(d, wkt = "geom") > print(dg, n = 1) Simple feature collection with 2 features and 1 field Geometry type: POINT Dimension: XY Bounding box: xmin: 0 ymin: 0 xmax: 1 ymax: 1 CRS: NA First 1 features: a geom 1 1 POINT (0 0) > head(st_as_sf(d, wkt = "geom"), 1) Simple feature collection with 1 feature and 1 field Geometry type: POINT Dimension: XY Bounding box: xmin: 0 ymin: 0 xmax: 0 ymax: 0 CRS: NA a geom 1 1 POINT (0 0) > > d$geom = st_as_sfc(d$geom) > d1 = d > attr(d1, "sf_col") = "geom" > st_geometry(d1) = d$geom > > d$geometry = d$geom # second geometry list-column > if (require(testthat, quietly = TRUE)) { + expect_warning(st_geometry(d) <- d$geom) + } > d Simple feature collection with 2 features and 1 field Active geometry column: geom Geometry type: POINT Dimension: XY Bounding box: xmin: 0 ymin: 0 xmax: 1 ymax: 1 CRS: NA a geom geometry 1 1 POINT (0 0) POINT (0 0) 2 2 POINT (1 1) POINT (1 1) > > x = st_sfc(list(st_point(0:1), st_point(0:1)), crs = 4326) > # don't warn when replacing crs with identical value: > st_sfc(x, crs = 4326) Geometry set for 2 features Geometry type: POINT Dimension: XY Bounding box: xmin: 0 ymin: 1 xmax: 0 ymax: 1 Geodetic CRS: WGS 84 POINT (0 1) POINT (0 1) > y = st_sfc(x, crs = "+proj=longlat +datum=WGS84 +no_defs") > # but do when it changes: > y = st_sfc(x, crs = 3857) Warning message: st_crs<- : replacing crs does not reproject data; use st_transform for that > > p = st_point(0:1) > st_cast(p, "MULTIPOINT") MULTIPOINT ((0 1)) > mp = st_multipoint(rbind(c(0,1), c(2,2))) > st_cast(mp, "POINT") POINT (0 1) Warning message: In st_cast.MULTIPOINT(mp, "POINT") : point from first coordinate only > st_cast(mp, "MULTIPOINT") MULTIPOINT ((0 1), (2 2)) > > # geometry collection to its elements: > st_cast(st_geometrycollection(list(mp)), "POINT") POINT (0 1) Warning message: In st_cast.MULTIPOINT(x[[1]], to, ...) : point from first coordinate only > st_cast(st_geometrycollection(list(mp)), "MULTIPOINT") MULTIPOINT ((0 1), (2 2)) > st_cast(st_geometrycollection(list(p,mp)), "MULTIPOINT") MULTIPOINT ((0 1)) Warning message: In st_cast.GEOMETRYCOLLECTION(st_geometrycollection(list(p, mp)), : only first part of geometrycollection is retained > > mp = st_multipoint(rbind(c(0,1))) > x = st_sfc(p, mp) > st_cast(x, "POINT") Geometry set for 2 features Geometry type: POINT Dimension: XY Bounding box: xmin: 0 ymin: 1 xmax: 0 ymax: 1 CRS: NA POINT (0 1) POINT (0 1) Warning message: In st_cast.MULTIPOINT(X[[i]], ...) : point from first coordinate only > > sf = st_sf(a = 3:2, geom = x) > st_cast(sf, "POINT") Simple feature collection with 2 features and 1 field Geometry type: POINT Dimension: XY Bounding box: xmin: 0 ymin: 1 xmax: 0 ymax: 1 CRS: NA a geom 1 3 POINT (0 1) 2 2 POINT (0 1) Warning message: In st_cast.MULTIPOINT(X[[i]], ...) : point from first coordinate only > > > x %>% st_cast("POINT") Geometry set for 2 features Geometry type: POINT Dimension: XY Bounding box: xmin: 0 ymin: 1 xmax: 0 ymax: 1 CRS: NA POINT (0 1) POINT (0 1) Warning message: In st_cast.MULTIPOINT(X[[i]], ...) : point from first coordinate only > > # points: > mp = st_multipoint(rbind(c(0,1))) # single-point multipoint > st_sfc(p,mp) %>% st_cast("POINT") Geometry set for 2 features Geometry type: POINT Dimension: XY Bounding box: xmin: 0 ymin: 1 xmax: 0 ymax: 1 CRS: NA POINT (0 1) POINT (0 1) Warning message: In st_cast.MULTIPOINT(X[[i]], ...) : point from first coordinate only > st_sfc(p,mp) %>% st_cast("MULTIPOINT") Geometry set for 2 features Geometry type: MULTIPOINT Dimension: XY Bounding box: xmin: 0 ymin: 1 xmax: 0 ymax: 1 CRS: NA MULTIPOINT ((0 1)) MULTIPOINT ((0 1)) > > # lines: > pts = rbind(c(0,0), c(1,1), c(2,1)) > st_sfc(st_linestring(pts), st_multilinestring(list(pts))) %>% st_cast("LINESTRING") Geometry set for 2 features Geometry type: LINESTRING Dimension: XY Bounding box: xmin: 0 ymin: 0 xmax: 2 ymax: 1 CRS: NA LINESTRING (0 0, 1 1, 2 1) LINESTRING (0 0, 1 1, 2 1) > st_sfc(st_linestring(pts), st_multilinestring(list(pts))) %>% st_cast("MULTILINESTRING") Geometry set for 2 features Geometry type: MULTILINESTRING Dimension: XY Bounding box: xmin: 0 ymin: 0 xmax: 2 ymax: 1 CRS: NA MULTILINESTRING ((0 0, 1 1, 2 1)) MULTILINESTRING ((0 0, 1 1, 2 1)) > > # polygons: > pts = rbind(c(0,0), c(1,1), c(0,1), c(0,0)) > st_sfc(st_polygon(list(pts)), st_multipolygon(list(list(pts)))) %>% st_cast("POLYGON") Geometry set for 2 features Geometry type: POLYGON Dimension: XY Bounding box: xmin: 0 ymin: 0 xmax: 1 ymax: 1 CRS: NA POLYGON ((0 0, 1 1, 0 1, 0 0)) POLYGON ((0 0, 1 1, 0 1, 0 0)) > st_sfc(st_polygon(list(pts)), st_multipolygon(list(list(pts)))) %>% st_cast("MULTIPOLYGON") Geometry set for 2 features Geometry type: MULTIPOLYGON Dimension: XY Bounding box: xmin: 0 ymin: 0 xmax: 1 ymax: 1 CRS: NA MULTIPOLYGON (((0 0, 1 1, 0 1, 0 0))) MULTIPOLYGON (((0 0, 1 1, 0 1, 0 0))) > > > st_sfc(st_geometrycollection(list(p)), st_geometrycollection(list(mp))) %>% st_cast() Geometry set for 2 features Geometry type: GEOMETRY Dimension: XY Bounding box: xmin: 0 ymin: 1 xmax: 0 ymax: 1 CRS: NA POINT (0 1) MULTIPOINT ((0 1)) > st_sfc(st_geometrycollection(list(p)), st_geometrycollection(list(mp))) %>% + st_cast() %>% + st_cast("POINT") Geometry set for 2 features Geometry type: POINT Dimension: XY Bounding box: xmin: 0 ymin: 1 xmax: 0 ymax: 1 CRS: NA POINT (0 1) POINT (0 1) Warning message: In st_cast.MULTIPOINT(X[[i]], ...) : point from first coordinate only > > p = rbind(c(0,0),c(1,0),c(1,1),c(0,1),c(0,0)) > pol = st_polygon(list(p)) > # plot(pol) > try(plot(st_polygonize(pol))) # --> breaks Error in st_polygonize.sfc(st_sfc(x)) : inherits(x, "sfc_LINESTRING") || inherits(x, "sfc_MULTILINESTRING") is not TRUE > st_length(st_sfc(st_point(c(0,0)))) [1] 0 > > try(as(st_sfc(st_linestring(matrix(1:9,3))), "Spatial")) Error in StopZ(zm) : sp supports Z dimension only for POINT and MULTIPOINT. use `st_zm(...)` to coerce to XY dimensions > > # check conus is present: > x = st_sfc(st_point(c(-90,35)), st_point(c(-80,36)), + crs = "+proj=longlat +datum=NAD27") > y = st_transform(x, 3857) > > ## IGNORE_RDIFF_BEGIN > sf_extSoftVersion()[1:3] GEOS GDAL proj.4 "3.12.1" "3.8.4" "9.3.1" > ## IGNORE_RDIFF_END > > # Ops.sfc: > ls = st_sfc(st_linestring(rbind(c(0,0),c(0,1)))) > ls * 2 Geometry set for 1 feature Geometry type: LINESTRING Dimension: XY Bounding box: xmin: 0 ymin: 0 xmax: 0 ymax: 2 CRS: NA LINESTRING (0 0, 0 2) > ls - 2 Geometry set for 1 feature Geometry type: LINESTRING Dimension: XY Bounding box: xmin: -2 ymin: -2 xmax: -2 ymax: -1 CRS: NA LINESTRING (-2 -2, -2 -1) > (ls + 2) %% 3 Geometry set for 1 feature Geometry type: LINESTRING Dimension: XY Bounding box: xmin: 2 ymin: 0 xmax: 2 ymax: 2 CRS: NA LINESTRING (2 2, 2 0) > ls / ls Geometry set for 1 feature (with 1 geometry empty) Geometry type: GEOMETRYCOLLECTION Dimension: XY Bounding box: xmin: NA ymin: NA xmax: NA ymax: NA CRS: NA GEOMETRYCOLLECTION EMPTY > p_ = st_point(0:1) > ll = st_sfc(ls[[1]], p_) > ll & st_sfc(p_) Geometry set for 2 features Geometry type: POINT Dimension: XY Bounding box: xmin: 0 ymin: 1 xmax: 0 ymax: 1 CRS: NA POINT (0 1) POINT (0 1) > ll | st_sfc(p_) Geometry set for 2 features Geometry type: GEOMETRY Dimension: XY Bounding box: xmin: 0 ymin: 0 xmax: 0 ymax: 1 CRS: NA LINESTRING (0 0, 0 1) POINT (0 1) > ll %/% st_sfc(p_) Geometry set for 2 features (with 1 geometry empty) Geometry type: GEOMETRY Dimension: XY Bounding box: xmin: 0 ymin: 0 xmax: 0 ymax: 1 CRS: NA LINESTRING (0 0, 0 1) GEOMETRYCOLLECTION EMPTY > ll == st_sfc(p_) [1] FALSE TRUE > ll != st_sfc(p_) [1] TRUE FALSE > > > str(x) sfc_POINT of length 2; first list element: 'XY' num [1:2] -90 35 > nc = st_read(system.file("shape/nc.shp", package="sf"), quiet = TRUE) > str(nc) Classes 'sf' and 'data.frame': 100 obs. of 15 variables: $ AREA : num 0.114 0.061 0.143 0.07 0.153 0.097 0.062 0.091 0.118 0.124 ... $ PERIMETER: num 1.44 1.23 1.63 2.97 2.21 ... $ CNTY_ : num 1825 1827 1828 1831 1832 ... $ CNTY_ID : num 1825 1827 1828 1831 1832 ... $ NAME : chr "Ashe" "Alleghany" "Surry" "Currituck" ... $ FIPS : chr "37009" "37005" "37171" "37053" ... $ FIPSNO : num 37009 37005 37171 37053 37131 ... $ CRESS_ID : int 5 3 86 27 66 46 15 37 93 85 ... $ BIR74 : num 1091 487 3188 508 1421 ... $ SID74 : num 1 0 5 1 9 7 0 0 4 1 ... $ NWBIR74 : num 10 10 208 123 1066 ... $ BIR79 : num 1364 542 3616 830 1606 ... $ SID79 : num 0 3 6 2 3 5 2 2 2 5 ... $ NWBIR79 : num 19 12 260 145 1197 ... $ geometry :sfc_MULTIPOLYGON of length 100; first list element: List of 1 ..$ :List of 1 .. ..$ : num [1:27, 1:2] -81.5 -81.5 -81.6 -81.6 -81.7 ... ..- attr(*, "class")= chr [1:3] "XY" "MULTIPOLYGON" "sfg" - attr(*, "sf_column")= chr "geometry" - attr(*, "agr")= Factor w/ 3 levels "constant","aggregate",..: NA NA NA NA NA NA NA NA NA NA ... ..- attr(*, "names")= chr [1:14] "AREA" "PERIMETER" "CNTY_" "CNTY_ID" ... > bb = st_as_sfc(st_bbox(nc)) > format(st_bbox(nc)) [1] "((-84.32385,33.88199),(-75.45698,36.58965))" > > st_agr("constant") [1] constant Levels: constant aggregate identity > st_agr() [1] Levels: constant aggregate identity > x <- st_sf(a = 1:2, b = 3:4, geom = x, agr = c("constant", "aggregate")) > y <- x %>% st_set_agr("constant") > y Simple feature collection with 2 features and 2 fields Attribute-geometry relationships: constant (2) Geometry type: POINT Dimension: XY Bounding box: xmin: -90 ymin: 35 xmax: -80 ymax: 36 Geodetic CRS: +proj=longlat +datum=NAD27 a b geom 1 1 3 POINT (-90 35) 2 2 4 POINT (-80 36) > > sf1 <- st_sf(a = c("x", "y"), geom = st_sfc(st_point(3:4), st_point(3:4))) > sf1[names(sf1)] Simple feature collection with 2 features and 1 field Geometry type: POINT Dimension: XY Bounding box: xmin: 3 ymin: 4 xmax: 3 ymax: 4 CRS: NA a geom 1 x POINT (3 4) 2 y POINT (3 4) > > st_bbox(sf1) xmin ymin xmax ymax 3 4 3 4 > bb = st_bbox(nc) > bb xmin ymin xmax ymax -84.32385 33.88199 -75.45698 36.58965 > st_crs(bb) Coordinate Reference System: User input: NAD27 wkt: GEOGCRS["NAD27", DATUM["North American Datum 1927", ELLIPSOID["Clarke 1866",6378206.4,294.978698213898, LENGTHUNIT["metre",1]]], PRIMEM["Greenwich",0, ANGLEUNIT["degree",0.0174532925199433]], CS[ellipsoidal,2], AXIS["latitude",north, ORDER[1], ANGLEUNIT["degree",0.0174532925199433]], AXIS["longitude",east, ORDER[2], ANGLEUNIT["degree",0.0174532925199433]], ID["EPSG",4267]] > st_bbox(c(xmin = 16.1, xmax = 16.6, ymin = 48.6, ymax = 47.9), crs = st_crs(4326)) xmin ymin xmax ymax 16.1 48.6 16.6 47.9 > st_bbox(c(xmin = 16.1, xmax = 16.6, ymin = 48.6, ymax = 47.9), crs = 4326) xmin ymin xmax ymax 16.1 48.6 16.6 47.9 > > bb$xrange xmin xmax -84.32385 -75.45698 > bb$yrange ymin ymax 33.88199 36.58965 > bb$xmin xmin -84.32385 > bb$ymin ymin 33.88199 > bb$xmax xmax -75.45698 > bb$ymax ymax 36.58965 > try(bb$foo) Error in `$.bbox`(bb, foo) : unsupported name > > # merge: > 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) Simple feature collection with 2 features and 3 fields Geometry type: POINT Dimension: XY Bounding box: xmin: 0 ymin: 0 xmax: 1 ymax: 1 CRS: NA b a x geometry 1 5 1 b POINT (0 0) 2 6 2 c POINT (1 1) > merge(a, b, all = TRUE) Simple feature collection with 4 features and 3 fields (with 1 geometry empty) Geometry type: GEOMETRY Dimension: XY Bounding box: xmin: 0 ymin: 0 xmax: 2 ymax: 2 CRS: NA b a x geometry 1 2 NA a GEOMETRYCOLLECTION EMPTY 2 5 1 b POINT (0 0) 3 6 2 c POINT (1 1) 4 7 3 POINT (2 2) > > if (require(dplyr, quietly = TRUE)) { + # joins: + inner_join(a, b) + left_join(a, b) + right_join(a, b) + full_join(a, b) + semi_join(a, b) + anti_join(a, b) + left_join(a, data.frame(b, geometry = 1), by = "b") + } Attaching package: 'dplyr' The following object is masked from 'package:testthat': matches The following objects are masked from 'package:stats': filter, lag The following objects are masked from 'package:base': intersect, setdiff, setequal, union Joining with `by = join_by(b)` Joining with `by = join_by(b)` Joining with `by = join_by(b)` Joining with `by = join_by(b)` Joining with `by = join_by(b)` Joining with `by = join_by(b)` Simple feature collection with 3 features and 4 fields Geometry type: POINT Dimension: XY Bounding box: xmin: 0 ymin: 0 xmax: 2 ymax: 2 CRS: NA a b x geometry.y geometry.x 1 1 5 b 1 POINT (0 0) 2 2 6 c 1 POINT (1 1) 3 3 7 NA POINT (2 2) > > # st_joins: > 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) Simple feature collection with 4 features and 2 fields Geometry type: POINT Dimension: XY Bounding box: xmin: 1 ymin: 1 xmax: 3 ymax: 3 CRS: NA a.x a.y geom 1 1 NA POINT (1 1) 2 2 12 POINT (2 2) 2.1 2 13 POINT (2 2) 3 3 14 POINT (3 3) > st_join(a, b, left = FALSE) Simple feature collection with 3 features and 2 fields Geometry type: POINT Dimension: XY Bounding box: xmin: 2 ymin: 2 xmax: 3 ymax: 3 CRS: NA a.x a.y geom 2 2 12 POINT (2 2) 2.1 2 13 POINT (2 2) 3 3 14 POINT (3 3) > # st_join, largest = TRUE: > nc <- st_transform(st_read(system.file("shape/nc.shp", package="sf"), quiet = TRUE), 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),] > st_join(nc, gr, largest = TRUE) Simple feature collection with 100 features and 16 fields Geometry type: MULTIPOLYGON Dimension: XY Bounding box: xmin: 406265 ymin: 48359.7 xmax: 3052877 ymax: 1044143 Projected CRS: NAD83 / North Carolina (ftUS) First 10 features: AREA PERIMETER CNTY_ CNTY_ID NAME FIPS FIPSNO CRESS_ID BIR74 SID74 1 0.114 1.442 1825 1825 Ashe 37009 37009 5 1091 1 2 0.061 1.231 1827 1827 Alleghany 37005 37005 3 487 0 3 0.143 1.630 1828 1828 Surry 37171 37171 86 3188 5 4 0.070 2.968 1831 1831 Currituck 37053 37053 27 508 1 5 0.153 2.206 1832 1832 Northampton 37131 37131 66 1421 9 6 0.097 1.670 1833 1833 Hertford 37091 37091 46 1452 7 7 0.062 1.547 1834 1834 Camden 37029 37029 15 286 0 8 0.091 1.284 1835 1835 Gates 37073 37073 37 420 0 9 0.118 1.421 1836 1836 Warren 37185 37185 93 968 4 10 0.124 1.428 1837 1837 Stokes 37169 37169 85 1612 1 NWBIR74 BIR79 SID79 NWBIR79 label col geometry 1 10 1364 0 19 A 4 #fb80724d MULTIPOLYGON (((1270813 913... 2 10 542 3 12 A 4 #fb80724d MULTIPOLYGON (((1340553 959... 3 208 3616 6 260 A 5 #80b1d34d MULTIPOLYGON (((1570586 910... 4 123 830 2 145 A 10 #bc80bd4d MULTIPOLYGON (((2881206 948... 5 1066 1606 3 1197 A 8 #fccde54d MULTIPOLYGON (((2525700 911... 6 954 1838 5 1237 A 9 #d9d9d94d MULTIPOLYGON (((2665112 911... 7 115 350 2 139 A 10 #bc80bd4d MULTIPOLYGON (((2881206 948... 8 254 594 2 371 A 9 #d9d9d94d MULTIPOLYGON (((2717988 951... 9 748 1190 2 844 A 8 #fccde54d MULTIPOLYGON (((2203888 914... 10 160 2038 5 176 A 5 #80b1d34d MULTIPOLYGON (((1697618 911... Warning message: attribute variables are assumed to be spatially constant throughout all geometries > > # rbind: > x = st_sf(a = 1:2, geom = st_sfc(list(st_point(0:1), st_point(0:1)), crs = 4326)) > rbind(x, x, x) Simple feature collection with 6 features and 1 field Geometry type: POINT Dimension: XY Bounding box: xmin: 0 ymin: 1 xmax: 0 ymax: 1 Geodetic CRS: WGS 84 a geom 1 1 POINT (0 1) 2 2 POINT (0 1) 3 1 POINT (0 1) 4 2 POINT (0 1) 5 1 POINT (0 1) 6 2 POINT (0 1) > nc2 = rbind(nc[1:50, ], nc[51:100, ]) > all.equal(nc, nc2) [1] TRUE > > # st_sample: > suppressWarnings(RNGversion("3.5.3")) > set.seed(131) > options(digits=6) > x = st_sfc(st_polygon(list(rbind(c(0,1),c(90,1),c(90,90),c(0,90),c(0,1)))), crs = st_crs(4326)) > (p <- st_sample(x, 10)) Geometry set for 10 features Geometry type: POINT Dimension: XY Bounding box: xmin: 11.2448 ymin: 3.16385 xmax: 82.3451 ymax: 60.1703 Geodetic CRS: WGS 84 First 5 geometries: POINT (18.5793 25.2416) POINT (11.2448 20.8596) POINT (26.3946 60.1703) POINT (33.8202 19.146) POINT (76.1712 32.1029) > p <- st_sample(x[[1]], 10) # sfg method > x = st_sfc(st_polygon(list(rbind(c(0,0),c(90,0),c(90,90),c(0,90),c(0,0))))) # NOT long/lat: > p <- st_sample(x, 10) > 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") # lwgeom breaks on this > # (p <- st_sample(x, 10)) > pt = st_multipoint(matrix(1:20,,2)) > st_sample(p, 3) Geometry set for 1 feature Geometry type: MULTIPOINT Dimension: XY Bounding box: xmin: 41.0557 ymin: 27.8024 xmax: 80.9558 ymax: 65.5424 CRS: NA MULTIPOINT ((57.2118 52.329), (80.9558 65.5424)... > try(st_sample(p, 3.3)) Geometry set for 1 feature Geometry type: MULTIPOINT Dimension: XY Bounding box: xmin: 19.3415 ymin: 9.17624 xmax: 80.9558 ymax: 65.5424 CRS: NA MULTIPOINT ((19.3415 18.6622), (57.1389 9.17624... Warning message: In st_sample.sfc(p, 3.3) : size is not an integer > 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) Geometry set for 4 features (with 1 geometry empty) Geometry type: MULTIPOINT Dimension: XY Bounding box: xmin: 0 ymin: 0 xmax: 0.0914728 ymax: 1 CRS: NA MULTIPOINT ((0 0.137868), (0 0.777635), (0 0.40... MULTIPOINT ((0.0914728 0), (0.073852 0), (0.003... MULTIPOINT ((0.0488222 1), (0.0716508 1), (0.08... MULTIPOINT EMPTY > st_sample(nc[1:2,], size = c(10,20)) Geometry set for 30 features Geometry type: POINT Dimension: XY Bounding box: xmin: 1216510 ymin: 935852 xmax: 1428880 ymax: 1031490 Projected CRS: NAD83 / North Carolina (ftUS) First 5 geometries: POINT (1264558 935852) POINT (1256623 949369) POINT (1222665 972624) POINT (1333528 968263) POINT (1310837 958630) > # try with LINES, LongLat, should generate a warning: > nc[1:2,] %>% st_transform(4326) %>% st_cast("MULTILINESTRING") %>% st_sample(size = c(10,20)) although coordinates are longitude/latitude, st_sample assumes that they are planar although coordinates are longitude/latitude, st_sample assumes that they are planar Geometry set for 2 features Geometry type: MULTIPOINT Dimension: XY Bounding box: xmin: -81.6893 ymin: 36.2491 xmax: -80.914 ymax: 36.5726 Geodetic CRS: WGS 84 MULTIPOINT ((-81.2401 36.3717), (-81.6893 36.36... MULTIPOINT ((-81.2629 36.4076), (-81.1375 36.56... > st_sample(ls, 80, type = "regular") Geometry set for 4 features (with 1 geometry empty) Geometry type: MULTIPOINT Dimension: XY Bounding box: xmin: 0 ymin: 0 xmax: 0.095454 ymax: 1 CRS: NA MULTIPOINT ((0 0.0104441), (0 0.0254442), (0 0.... MULTIPOINT ((0.000452334 0), (0.0154525 0), (0.... MULTIPOINT ((0.00545321 1), (0.0204533 1), (0.0... MULTIPOINT EMPTY > p_sample = lapply(1:10, function(i) st_sample(nc[i, ], 100, exact = FALSE)) > lengths(p_sample) [1] 104 106 110 96 98 120 87 105 104 99 > p_sample_exact = lapply(1:10, function(i) st_sample(nc[i, ], 100, exact = TRUE)) > lengths(p_sample_exact) [1] 100 100 100 100 100 100 100 100 100 100 > #plot(nc$geometry[1]) > #plot(p_sample[[1]], add = TRUE) > #plot(p_sample_exact[[1]], add = TRUE) > > if (require(dplyr, quietly = TRUE)) { + #class(st_bind_cols(nc, as.data.frame(nc)[1:3])) + print(class(dplyr::bind_cols(nc, as.data.frame(nc)[1:3]))) + } New names: • `AREA` -> `AREA...1` • `PERIMETER` -> `PERIMETER...2` • `CNTY_` -> `CNTY_...3` • `AREA` -> `AREA...16` • `PERIMETER` -> `PERIMETER...17` • `CNTY_` -> `CNTY_...18` [1] "sf" "data.frame" > class(rbind(nc, nc)) [1] "sf" "data.frame" > class(cbind(nc, nc)) [1] "sf" "data.frame" > > x = st_sfc(st_point(0:1), st_point(2:3)) > x[c(NA,1,NA,2,NA)] Geometry set for 5 features (with 3 geometries empty) Geometry type: POINT Dimension: XY Bounding box: xmin: 0 ymin: 1 xmax: 2 ymax: 3 CRS: NA POINT EMPTY POINT (0 1) POINT EMPTY POINT (2 3) POINT EMPTY > > # jitter > 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') > st_jitter(st_sfc(st_point(0:1)), amount = .1) Geometry set for 1 feature Geometry type: POINT Dimension: XY Bounding box: xmin: -0.0500922 ymin: 0.992953 xmax: -0.0500922 ymax: 0.992953 CRS: NA POINT (-0.0500922 0.992953) > > # st_bbox: > if (suppressPackageStartupMessages(require(sp, quietly = TRUE)) && require(raster, quietly = TRUE)) { + demo(meuse, ask = FALSE, echo = FALSE) + suppressWarnings(st_bbox(meuse)) + crs = suppressWarnings(st_crs(meuse)) + suppressWarnings(st_bbox(raster(meuse.grid))) + st_bbox(extent(raster())) + } Attaching package: 'raster' The following object is masked from 'package:dplyr': select xmin ymin xmax ymax -180 -90 180 90 > > # st_to_s2 > if (FALSE) { # stops working with GDAL 2.3.0 / PROJ 5.0.1: + x = sf:::st_to_s2(nc) + x1 = st_geometry(x) + cc = st_coordinates(x1) + summary(sqrt(cc[,1]^2+cc[,2]^2+cc[,3]^2)) + } > > # check_ring_dir > m = rbind(c(0,0), c(0,1), c(1,1), c(1,0), c(0,0)) > mi = m[nrow(m):1,] > pol = st_polygon(list(m * 10, m + .5, mi + 1.5, mi + 3.5, m + 5, mi + 6.5)) > st_sfc(pol) Geometry set for 1 feature Geometry type: POLYGON Dimension: XY Bounding box: xmin: 0 ymin: 0 xmax: 10 ymax: 10 CRS: NA POLYGON ((0 0, 0 10, 10 10, 10 0, 0 0), (0.5 0.... > x = st_sfc(pol, check_ring_dir=TRUE) > y = st_sf(a = 1, geom = st_sfc(pol), check_ring_dir=TRUE) > str(x) sfc_POLYGON of length 1; first list element: List of 6 $ : num [1:5, 1:2] 0 10 10 0 0 0 0 10 10 0 $ : num [1:5, 1:2] 0.5 0.5 1.5 1.5 0.5 0.5 1.5 1.5 0.5 0.5 $ : num [1:5, 1:2] 1.5 1.5 2.5 2.5 1.5 1.5 2.5 2.5 1.5 1.5 $ : num [1:5, 1:2] 3.5 3.5 4.5 4.5 3.5 3.5 4.5 4.5 3.5 3.5 $ : num [1:5, 1:2] 5 5 6 6 5 5 6 6 5 5 $ : num [1:5, 1:2] 6.5 6.5 7.5 7.5 6.5 6.5 7.5 7.5 6.5 6.5 - attr(*, "class")= chr [1:3] "XY" "POLYGON" "sfg" > x = st_sfc(st_polygon(), st_polygon(), check_ring_dir=TRUE) > str(x) sfc_POLYGON of length 2; first list element: list() - attr(*, "class")= chr [1:3] "XY" "POLYGON" "sfg" > # empty ring/zero area: > x = st_sfc(st_polygon(list(m[c(1,3,1),])), check_ring_dir=TRUE) > > mp = st_multipolygon(list(pol, pol)) > try(x <- st_sfc(mp, st_polygon(), check_ring_dir=TRUE)) Error in check_ring_dir(lst) : check_ring_dir: not supported for class sfc_GEOMETRY > x <- st_sfc(mp, pol) %>% st_cast("MULTIPOLYGON") %>% st_sfc(check_ring_dir=TRUE) > x Geometry set for 2 features Geometry type: MULTIPOLYGON Dimension: XY Bounding box: xmin: 0 ymin: 0 xmax: 10 ymax: 10 CRS: NA MULTIPOLYGON (((0 0, 10 0, 10 10, 0 10, 0 0), (... MULTIPOLYGON (((0 0, 10 0, 10 10, 0 10, 0 0), (... > str(x) sfc_MULTIPOLYGON of length 2; first list element: List of 2 $ :List of 6 ..$ : num [1:5, 1:2] 0 10 10 0 0 0 0 10 10 0 ..$ : num [1:5, 1:2] 0.5 0.5 1.5 1.5 0.5 0.5 1.5 1.5 0.5 0.5 ..$ : num [1:5, 1:2] 1.5 1.5 2.5 2.5 1.5 1.5 2.5 2.5 1.5 1.5 ..$ : num [1:5, 1:2] 3.5 3.5 4.5 4.5 3.5 3.5 4.5 4.5 3.5 3.5 ..$ : num [1:5, 1:2] 5 5 6 6 5 5 6 6 5 5 ..$ : num [1:5, 1:2] 6.5 6.5 7.5 7.5 6.5 6.5 7.5 7.5 6.5 6.5 ..- attr(*, "class")= chr [1:3] "XY" "POLYGON" "sfg" $ :List of 6 ..$ : num [1:5, 1:2] 0 10 10 0 0 0 0 10 10 0 ..$ : num [1:5, 1:2] 0.5 0.5 1.5 1.5 0.5 0.5 1.5 1.5 0.5 0.5 ..$ : num [1:5, 1:2] 1.5 1.5 2.5 2.5 1.5 1.5 2.5 2.5 1.5 1.5 ..$ : num [1:5, 1:2] 3.5 3.5 4.5 4.5 3.5 3.5 4.5 4.5 3.5 3.5 ..$ : num [1:5, 1:2] 5 5 6 6 5 5 6 6 5 5 ..$ : num [1:5, 1:2] 6.5 6.5 7.5 7.5 6.5 6.5 7.5 7.5 6.5 6.5 ..- attr(*, "class")= chr [1:3] "XY" "POLYGON" "sfg" - attr(*, "class")= chr [1:3] "XY" "MULTIPOLYGON" "sfg" > > x = st_sfc(st_linestring(rbind(c(-179,0),c(179,0))), crs = 4326) > st_wrap_dateline(st_sf(a = 1, geometry = x)) Simple feature collection with 1 feature and 1 field Geometry type: MULTILINESTRING Dimension: XY Bounding box: xmin: -180 ymin: 0 xmax: 180 ymax: 0 Geodetic CRS: WGS 84 a geometry 1 1 MULTILINESTRING ((-179 0, -... > st_wrap_dateline(x) Geometry set for 1 feature Geometry type: MULTILINESTRING Dimension: XY Bounding box: xmin: -180 ymin: 0 xmax: 180 ymax: 0 Geodetic CRS: WGS 84 MULTILINESTRING ((-179 0, -180 0), (180 0, 179 0)) > st_wrap_dateline(x[[1]]) MULTILINESTRING ((-179 0, -180 0), (180 0, 179 0)) > > geo <- c("{\"geodesic\":true,\"type\":\"Point\",\"coordinates\":[-118.68152563269095,36.43764870908927]}", + "{\"geodesic\":true,\"type\":\"Point\",\"coordinates\":[-118.67408758213843,36.43366018922779]}", + "{\"geodesic\":true,\"type\":\"Point\",\"coordinates\":[-118.67708346361097,36.44208638659282]}", + "{\"geodesic\":true,\"type\":\"Point\",\"coordinates\":[-118.67886661944996,36.44110273135671]}", + "{\"geodesic\":true,\"type\":\"Point\",\"coordinates\":[-118.68089232041565,36.44173155205561]}") > st_as_sfc(geo, GeoJSON = TRUE) Geometry set for 5 features Geometry type: POINT Dimension: XY Bounding box: xmin: -118.682 ymin: 36.4337 xmax: -118.674 ymax: 36.4421 Geodetic CRS: WGS 84 POINT (-118.682 36.4376) POINT (-118.674 36.4337) POINT (-118.677 36.4421) POINT (-118.679 36.4411) POINT (-118.681 36.4417) > st_as_sfc(geo, GeoJSON = TRUE, crs = 4326) Geometry set for 5 features Geometry type: POINT Dimension: XY Bounding box: xmin: -118.682 ymin: 36.4337 xmax: -118.674 ymax: 36.4421 Geodetic CRS: WGS 84 POINT (-118.682 36.4376) POINT (-118.674 36.4337) POINT (-118.677 36.4421) POINT (-118.679 36.4411) POINT (-118.681 36.4417) > > st_as_sfc(st_as_binary(st_sfc(st_point(0:1)))[[1]], crs = 4326) Geometry set for 1 feature Geometry type: POINT Dimension: XY Bounding box: xmin: 0 ymin: 1 xmax: 0 ymax: 1 Geodetic CRS: WGS 84 POINT (0 1) > > x = nc > x$geom = NULL > class(x) [1] "sf" "data.frame" > > st_as_sfc(list(st_point(0:1)), crs = 4326) Geometry set for 1 feature Geometry type: POINT Dimension: XY Bounding box: xmin: 0 ymin: 1 xmax: 0 ymax: 1 Geodetic CRS: WGS 84 POINT (0 1) > > # crop: > box = c(xmin = 0, ymin = 0, xmax = 1, ymax = 1) > > pol = st_sfc(st_buffer(st_point(c(.5, .5)), .65)) > pol_sf = st_sf(a=1, geom=pol) > > st_crop(pol, box) Geometry set for 1 feature Geometry type: POLYGON Dimension: XY Bounding box: xmin: 0 ymin: 0 xmax: 1 ymax: 1 CRS: NA POLYGON ((0.983044 0.0650651, 0.959619 0.040380... > st_crop(pol, st_bbox(box)) Geometry set for 1 feature Geometry type: POLYGON Dimension: XY Bounding box: xmin: 0 ymin: 0 xmax: 1 ymax: 1 CRS: NA POLYGON ((0.983044 0.0650651, 0.959619 0.040380... > st_crop(pol_sf, box) Simple feature collection with 1 feature and 1 field Geometry type: POLYGON Dimension: XY Bounding box: xmin: 0 ymin: 0 xmax: 1 ymax: 1 CRS: NA a geom 1 1 POLYGON ((0.983044 0.065065... Warning message: attribute variables are assumed to be spatially constant throughout all geometries > st_crop(pol_sf, st_bbox(box)) Simple feature collection with 1 feature and 1 field Geometry type: POLYGON Dimension: XY Bounding box: xmin: 0 ymin: 0 xmax: 1 ymax: 1 CRS: NA a geom 1 1 POLYGON ((0.983044 0.065065... Warning message: attribute variables are assumed to be spatially constant throughout all geometries > > # new sample methods: > x = st_sfc(st_polygon(list(rbind(c(0,0),c(90,0),c(90,90),c(0,90),c(0,0))))) # NOT long/lat: > p <- st_sample(x, 10, type = "regular") > p <- st_sample(x, 10, type = "hexagonal") > > all.equal(st_drop_geometry(pol_sf), st_set_geometry(pol_sf, NULL)) [1] TRUE > > # https://github.com/r-spatial/sf/issues/1024 > shape1 <-st_sfc(st_polygon(list(rbind(c(0,0), c(1,0), c(3,2), c(2,4), c(1,4), c(0,0))))) > shape2 <- st_sfc(st_polygon()) > shape3 <- st_sfc(st_polygon()) > > shape4 = st_intersection(shape2, shape3) # has zero features > > st_difference(shape1, shape4) Geometry set for 1 feature Geometry type: POLYGON Dimension: XY Bounding box: xmin: 0 ymin: 0 xmax: 3 ymax: 4 CRS: NA POLYGON ((0 0, 1 0, 3 2, 2 4, 1 4, 0 0)) > st_difference(shape4, shape1) Geometry set for 0 features Bounding box: xmin: NA ymin: NA xmax: NA ymax: NA CRS: NA > st_sym_difference(shape1, shape4) Geometry set for 1 feature Geometry type: POLYGON Dimension: XY Bounding box: xmin: 0 ymin: 0 xmax: 3 ymax: 4 CRS: NA POLYGON ((0 0, 1 0, 3 2, 2 4, 1 4, 0 0)) > st_union(shape1, shape4) Geometry set for 1 feature Geometry type: POLYGON Dimension: XY Bounding box: xmin: 0 ymin: 0 xmax: 3 ymax: 4 CRS: NA POLYGON ((0 0, 1 0, 3 2, 2 4, 1 4, 0 0)) > st_union(shape4, shape1) Geometry set for 1 feature Geometry type: POLYGON Dimension: XY Bounding box: xmin: 0 ymin: 0 xmax: 3 ymax: 4 CRS: NA POLYGON ((0 0, 1 0, 3 2, 2 4, 1 4, 0 0)) > > # transform empty: > tr = st_sf(geom=st_sfc()) %>% st_set_crs(3587) %>% st_transform(4326) > > # NA values are converted to empty; #1114: > x <- data.frame(name=LETTERS) > y <- data.frame(name=LETTERS[1:13], letters[14:26]) > y$geometry <- st_sfc(st_point(c(0,0))) > y <- st_sf(y) > out = merge(x, y, all.x=TRUE) > class(out) [1] "data.frame" > > st_as_sf(st_sfc(st_point(0:1))) Simple feature collection with 1 feature and 0 fields Geometry type: POINT Dimension: XY Bounding box: xmin: 0 ymin: 1 xmax: 0 ymax: 1 CRS: NA x 1 POINT (0 1) > > # st_exterior_ring(): > 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) > mpl1 = st_multipolygon(list(pl1,pl1+20)) > > spl1 = st_as_sfc(list(pl1),crs=4326) > smpl1 = st_as_sfc(list(mpl1),crs=4326) > > st_exterior_ring(spl1[[1]]) POLYGON ((0 0, 10 0, 10 10, 0 10, 0 0)) > st_exterior_ring(spl1) Geometry set for 1 feature Geometry type: POLYGON Dimension: XY Bounding box: xmin: 0 ymin: 0 xmax: 10 ymax: 10 Geodetic CRS: WGS 84 POLYGON ((0 0, 10 0, 10 10, 0 10, 0 0)) > st_exterior_ring(st_sf(a = 1, geom = spl1)) Simple feature collection with 1 feature and 1 field Geometry type: POLYGON Dimension: XY Bounding box: xmin: 0 ymin: 0 xmax: 10 ymax: 10 Geodetic CRS: WGS 84 a geom 1 1 POLYGON ((0 0, 10 0, 10 10,... > st_exterior_ring(smpl1[[1]]) MULTIPOLYGON (((0 0, 10 0, 10 10, 0 10, 0 0)), ((20 20, 30 20, 30 30, 20 30, 20 20))) > st_exterior_ring(st_sfc(smpl1)) Geometry set for 1 feature Geometry type: MULTIPOLYGON Dimension: XY Bounding box: xmin: 0 ymin: 0 xmax: 30 ymax: 30 Geodetic CRS: WGS 84 MULTIPOLYGON (((0 0, 10 0, 10 10, 0 10, 0 0)), ... > st_exterior_ring(st_sf(a = 1, geom = st_sfc(smpl1))) Simple feature collection with 1 feature and 1 field Geometry type: MULTIPOLYGON Dimension: XY Bounding box: xmin: 0 ymin: 0 xmax: 30 ymax: 30 Geodetic CRS: WGS 84 a geom 1 1 MULTIPOLYGON (((0 0, 10 0, ... > > '{"type":"Polygon","coordinates":[[]]}' |> read_sf() |> st_is_empty() [1] TRUE > # '{"type":"Polygon","coordinates":[]}' |> read_sf() |> st_is_empty() # breaks on GDAL < 3.9 or so > '{"type":"MultiPolygon","coordinates":[[[]]]}' |> read_sf() |> st_is_empty() [1] TRUE > '{"type":"MultiPolygon","coordinates":[[]]}' |> read_sf() |> st_is_empty() [1] TRUE > > proc.time() user system elapsed 7.23 0.23 7.46