R version 4.2.1 (2022-06-23) -- "Funny-Looking Kid"
Copyright (C) 2022 The R Foundation for Statistical Computing
Platform: x86_64-pc-linux-gnu (64-bit)
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
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)
>
> sf_extSoftVersion()[1:3]
GEOS GDAL proj.4
"3.10.2" "3.4.3" "8.2.1"
>
> # 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] <NA>
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 relationship: 2 constant, 0 aggregate, 0 identity
Geometry type: POINT
Dimension: XY
Bounding box: xmin: -90 ymin: 35 xmax: -80 ymax: 36
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 <NA> 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, by = "b"
Joining, by = "b"
Joining, by = "b"
Joining, by = "b"
Joining, by = "b"
Joining, 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> 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: 406262.2 ymin: 48374.87 xmax: 3052887 ymax: 1044158
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 (((1270814 913...
2 10 542 3 12 A 4 #fb80724d MULTIPOLYGON (((1340555 959...
3 208 3616 6 260 A 5 #80b1d34d MULTIPOLYGON (((1570591 910...
4 123 830 2 145 A 10 #bc80bd4d MULTIPOLYGON (((2881208 948...
5 1066 1606 3 1197 A 8 #fccde54d MULTIPOLYGON (((2525701 911...
6 954 1838 5 1237 A 9 #d9d9d94d MULTIPOLYGON (((2665113 911...
7 115 350 2 139 A 10 #bc80bd4d MULTIPOLYGON (((2881208 948...
8 254 594 2 371 A 9 #d9d9d94d MULTIPOLYGON (((2717989 951...
9 748 1190 2 844 A 8 #fccde54d MULTIPOLYGON (((2203890 914...
10 160 2038 5 176 A 5 #80b1d34d MULTIPOLYGON (((1697624 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: 935867 xmax: 1428880 ymax: 1031510
Projected CRS: NAD83 / North Carolina (ftUS)
First 5 geometries:
POINT (1264558 935867)
POINT (1256623 949384)
POINT (1222665 972639)
POINT (1333530 968278)
POINT (1310838 958645)
> # 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.2492 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.1374 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 (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
Warning message:
In showSRID(uprojargs, format = "PROJ", multiline = "NO", prefer_proj = prefer_proj) :
Discarded datum Amersfoort in Proj4 definition
>
> # 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)
>
> proc.time()
user system elapsed
8.907 0.852 8.874
Generated by dwww version 1.15 on Sun Jun 30 09:21:38 CEST 2024.