R version 4.1.3 (2022-03-10) -- "One Push-Up"
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.
> Sys.setenv(TZ="UTC")
> suppressPackageStartupMessages(library(sf))
> if ("GPKG" %in% st_drivers()$name) {
+ tst = st_read(system.file("gpkg/nc.gpkg", package="sf"), "nc.gpkg", crs = 4267, quiet = TRUE)
+ tst = st_read(system.file("gpkg/nc.gpkg", package="sf"), quiet = TRUE) # default layer name
+ }
>
> tst = st_read(system.file("shape/nc.shp", package="sf"), quiet = TRUE) # no layer specified
>
> # data from https://github.com/edzer/sfr/issues/6
> tst = st_read(system.file("shape/storms_xyz.shp", package="sf"), quiet = TRUE)
> class(st_geometry(tst))
[1] "sfc_LINESTRING" "sfc"
> class(st_geometry(tst)[[1]])
[1] "XYZ" "LINESTRING" "sfg"
> tst = st_read(system.file("shape/storms_xyzm.shp", package="sf"), quiet = TRUE)
> class(st_geometry(tst))
[1] "sfc_LINESTRING" "sfc"
> class(st_geometry(tst)[[1]])
[1] "XYM" "LINESTRING" "sfg"
> tst = st_read(system.file("shape/storms_xyz_feature.shp", package="sf"), quiet = TRUE)
> class(st_geometry(tst))
[1] "sfc_LINESTRING" "sfc"
> class(st_geometry(tst)[[1]])
[1] "XYZ" "LINESTRING" "sfg"
> tst = st_read(system.file("shape/storms_xyzm_feature.shp", package="sf"), quiet = TRUE)
> class(st_geometry(tst))
[1] "sfc_LINESTRING" "sfc"
> class(st_geometry(tst)[[1]])
[1] "XYM" "LINESTRING" "sfg"
>
> if ("GPKG" %in% st_drivers()$name) { # read Int64
+ print(st_read(system.file("gpkg/tl.gpkg", package="sf"), quiet = TRUE)$AWATER)
+ print(st_read(system.file("gpkg/tl.gpkg", package="sf"), quiet = TRUE, int64_as_string = TRUE)$AWATER)
+ }
[1] 1028678842
[1] "1028678842"
>
> # see https://github.com/edzer/sfr/issues/45 :
> if ("OSM" %in% st_drivers()$name && Sys.info()['sysname'] != "Darwin") {
+ osm = system.file("osm/overpass.osm", package="sf")
+ print(st_layers(osm))
+ suppressWarnings(print(st_layers(osm, do_count = TRUE)))
+ suppressWarnings(print(st_read(osm, "multipolygons", quiet = TRUE)))
+ }
Driver: OSM
Available layers:
layer_name geometry_type features fields crs_name
1 points Point NA 10 WGS 84
2 lines Line String NA 9 WGS 84
3 multilinestrings Multi Line String NA 4 WGS 84
4 multipolygons Multi Polygon NA 25 WGS 84
5 other_relations Geometry Collection NA 4 WGS 84
Driver: OSM
Available layers:
layer_name geometry_type features fields crs_name
1 points Point 0 10 WGS 84
2 lines Line String 0 9 WGS 84
3 multilinestrings Multi Line String 0 4 WGS 84
4 multipolygons Multi Polygon 3 25 WGS 84
5 other_relations Geometry Collection 0 4 WGS 84
Simple feature collection with 3 features and 25 fields
Geometry type: MULTIPOLYGON
Dimension: XY
Bounding box: xmin: -1.556254 ymin: 53.79519 xmax: -1.550716 ymax: 53.80341
Geodetic CRS: WGS 84
osm_id osm_way_id name type aeroway amenity admin_level barrier boundary
1 <NA> 263958111 <NA> <NA> <NA> <NA> <NA> <NA> <NA>
2 <NA> 301461256 <NA> <NA> <NA> <NA> <NA> <NA> <NA>
3 <NA> 462227487 <NA> <NA> <NA> <NA> <NA> <NA> <NA>
building craft geological historic land_area landuse leisure man_made
1 <NA> <NA> <NA> <NA> <NA> <NA> park <NA>
2 <NA> <NA> <NA> <NA> <NA> <NA> park <NA>
3 <NA> <NA> <NA> <NA> <NA> <NA> park <NA>
military natural office place shop sport tourism other_tags
1 <NA> <NA> <NA> <NA> <NA> <NA> <NA> <NA>
2 <NA> <NA> <NA> <NA> <NA> <NA> <NA> <NA>
3 <NA> <NA> <NA> <NA> <NA> <NA> <NA> <NA>
geometry
1 MULTIPOLYGON (((-1.550719 5...
2 MULTIPOLYGON (((-1.555763 5...
3 MULTIPOLYGON (((-1.551502 5...
>
> # layer opening option:
> st_read(system.file("shape/nc.shp", package="sf"), quiet = TRUE,
+ options = c("ADJUST_TYPE=YES", "ENCODING=CPL_ENC_UTF8"))
Simple feature collection with 100 features and 14 fields
Geometry type: MULTIPOLYGON
Dimension: XY
Bounding box: xmin: -84.32385 ymin: 33.88199 xmax: -75.45698 ymax: 36.58965
Geodetic CRS: NAD27
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 geometry
1 10 1364 0 19 MULTIPOLYGON (((-81.47276 3...
2 10 542 3 12 MULTIPOLYGON (((-81.23989 3...
3 208 3616 6 260 MULTIPOLYGON (((-80.45634 3...
4 123 830 2 145 MULTIPOLYGON (((-76.00897 3...
5 1066 1606 3 1197 MULTIPOLYGON (((-77.21767 3...
6 954 1838 5 1237 MULTIPOLYGON (((-76.74506 3...
7 115 350 2 139 MULTIPOLYGON (((-76.00897 3...
8 254 594 2 371 MULTIPOLYGON (((-76.56251 3...
9 748 1190 2 844 MULTIPOLYGON (((-78.30876 3...
10 160 2038 5 176 MULTIPOLYGON (((-80.02567 3...
>
> x <- st_sf(a = 1:2, geom = st_sfc(st_point(0:1), st_multipoint(matrix(1:4,2,2))))
>
> try(st_layers("foo")) # cannot open datasource
Cannot open data source foo
Error in CPL_get_layers(dsn, options, do_count) : Open failed.
> try(st_read("foo")) # cannot open datasource
Error : Cannot open "foo"; The file doesn't seem to exist.
> try(st_read("x.gpkg", "xyz")) # cannot open layer
Error : Cannot open "x.gpkg"; The file doesn't seem to exist.
> try(st_write(c("foo", "bar")))
Error in UseMethod("st_write") :
no applicable method for 'st_write' applied to an object of class "character"
> try(st_write(x, c("foo", "bar")))
Error in guess_driver(dns) : length(dsn) == 1 is not TRUE
> try(st_write(x, "foo", driver = "foo"))
driver `foo' not available.
Error in CPL_write_ogr(obj, dsn, layer, driver, as.character(dataset_options), :
Driver not available.
> if (Sys.getenv("USER") == "travis") {
+ try(st_write(x, "/x", driver = "ESRI Shapefile"))
+ st_write(x, "xxx.gpkg")
+ st_write(x, "xxx.gpkg", append = TRUE, quiet = FALSE) # appends to layer
+ y <- st_sf(b = 1:2, geom = st_sfc(st_point(0:1), st_multipoint(matrix(1:4,2,2))))
+ try(st_write(y, "xxx.gpkg"))
+ }
>
> geom = st_sfc(st_point(0:1), st_multipoint(matrix(1:4,2,2)))
> st_write(geom, "geom.gpkg")
writing GPKG: substituting LOCAL_CS["Undefined Cartesian SRS"] for missing CRS
Writing layer `geom' to data source `geom.gpkg' using driver `GPKG'
Writing 2 features with 0 fields and geometry type Unknown (any).
> st_write(geom, "geom1.gpkg", layer = "foo")
writing GPKG: substituting LOCAL_CS["Undefined Cartesian SRS"] for missing CRS
Writing layer `foo' to data source `geom1.gpkg' using driver `GPKG'
Writing 2 features with 0 fields and geometry type Unknown (any).
>
> df <- data.frame(
+ a = c(0, 1, NA, -Inf, Inf),
+ b = c("a", "b", NA, "c", ""),
+ c = c(as.Date("2001-01-01"),
+ as.Date(c(NA, -99, 0, 1), origin = "1970-01-01")),
+ d = c(as.POSIXct("2001-01-01"),
+ as.POSIXct(c(NA, -99, 0, 1), origin = "1970-01-01")),
+ x = 1:5,
+ y = 1:5)
>
> x = st_as_sf(df, coords = c("x", "y"))
>
> if ("GPKG" %in% st_drivers()$name) {
+ st_write(x, "x2.gpkg", quiet = TRUE)
+ y = st_read("x2.gpkg", quiet = TRUE)
+ print(y)
+ }
writing GPKG: substituting LOCAL_CS["Undefined Cartesian SRS"] for missing CRS
Simple feature collection with 5 features and 4 fields
Geometry type: POINT
Dimension: XY
Bounding box: xmin: 1 ymin: 1 xmax: 5 ymax: 5
Projected CRS: Undefined Cartesian SRS
a b c d geom
1 0 a 2001-01-01 2001-01-01 00:00:00 POINT (1 1)
2 1 b <NA> <NA> POINT (2 2)
3 NA <NA> 1969-09-24 1969-12-31 23:58:21 POINT (3 3)
4 -Inf c 1970-01-01 1970-01-01 00:00:00 POINT (4 4)
5 Inf 1970-01-02 1970-01-01 00:00:01 POINT (5 5)
>
> if ("SQLite" %in% st_drivers()$name && require(RSQLite)) {
+ db = system.file("sqlite/meuse.sqlite", package = "sf")
+ dbcon <- dbConnect(dbDriver("SQLite"), db)
+ m = dbReadTable(dbcon, "meuse.sqlite")
+ m$GEOMETRY = st_as_sfc(m$GEOMETRY, spatialite = FALSE) # ISO wkb
+ print(st_sf(m), n = 3)
+ # or:
+ (s = st_read(dbcon, "meuse.sqlite"))[1:3,]
+ dbDisconnect(dbcon)
+
+ db = system.file("sqlite/nc.sqlite", package = "sf")
+ dbcon <- dbConnect(dbDriver("SQLite"), db)
+ m = dbReadTable(dbcon, "nc.sqlite")
+ m$GEOMETRY = st_as_sfc(m$GEOMETRY, spatialite = FALSE) # ISO wkb
+ print(st_sf(m), n = 3)
+ dbDisconnect(dbcon)
+
+ db = system.file("sqlite/b.sqlite", package = "sf") # has an INT8 field
+ b = st_read(db, quiet = TRUE)
+ print(b)
+ b = st_read(db, int64_as_string = TRUE, quiet = TRUE)
+ print(b)
+
+ if (.Platform$endian == "little") {
+ db = system.file("sqlite/test3.sqlite", package = "sf")
+ dbcon <- dbConnect(dbDriver("SQLite"), db)
+ m = dbReadTable(dbcon, "HighWays")
+ m$Geometry = st_as_sfc(m$Geometry, spatialite = TRUE) # spatialite wkb
+ print(st_sf(m), n = 3)
+ m = dbReadTable(dbcon, "Towns")
+ m$Geometry = st_as_sfc(m$Geometry, spatialite = TRUE) # spatialite wkb
+ print(st_sf(m), n = 3)
+ m = dbReadTable(dbcon, "Regions")
+ m$Geometry = st_as_sfc(m$Geometry, spatialite = TRUE) # spatialite wkb
+ print(st_sf(m), n = 3)
+ }
+ }
Loading required package: RSQLite
Simple feature collection with 155 features and 13 fields
Geometry type: POINT
Dimension: XY
Bounding box: xmin: 178605 ymin: 329714 xmax: 181390 ymax: 333611
CRS: NA
First 3 features:
ogc_fid cadmium copper lead zinc elev dist om ffreq soil lime
1 1 11.7 85 299 1022 7.909 0.00135803 13.6 1 1 1
2 2 8.6 81 277 1141 6.983 0.01222430 14.0 1 1 1
3 3 6.5 68 199 640 7.800 0.10302900 13.0 1 1 1
landuse dist.m GEOMETRY
1 Ah 50 POINT (181072 333611)
2 Ah 30 POINT (181025 333558)
3 Ah 150 POINT (181165 333537)
Simple feature collection with 100 features and 15 fields
Geometry type: MULTIPOLYGON
Dimension: XY
Bounding box: xmin: -84.32385 ymin: 33.88199 xmax: -75.45698 ymax: 36.58965
CRS: NA
First 3 features:
ogc_fid area perimeter cnty_ cnty_id name fips fipsno cress_id bir74
1 1 0.114 1.442 1825 1825 Ashe 37009 37009 5 1091
2 2 0.061 1.231 1827 1827 Alleghany 37005 37005 3 487
3 3 0.143 1.630 1828 1828 Surry 37171 37171 86 3188
sid74 nwbir74 bir79 sid79 nwbir79 GEOMETRY
1 1 10 1364 0 19 MULTIPOLYGON (((-81.47276 3...
2 0 10 542 3 12 MULTIPOLYGON (((-81.23989 3...
3 5 208 3616 6 260 MULTIPOLYGON (((-80.45634 3...
Integer64 values larger than 9.0072e+15 lost significance after conversion to double;
use argument int64_as_string = TRUE to import them lossless, as character
Simple feature collection with 1 feature and 2 fields
Geometry type: POINT
Dimension: XY
Bounding box: xmin: 0 ymin: 1 xmax: 0 ymax: 1
CRS: NA
a bigint GEOMETRY
1 1 4.611686e+18 POINT (0 1)
Simple feature collection with 1 feature and 2 fields
Geometry type: POINT
Dimension: XY
Bounding box: xmin: 0 ymin: 1 xmax: 0 ymax: 1
CRS: NA
a bigint GEOMETRY
1 1 4611686018427387904 POINT (0 1)
Simple feature collection with 10 features and 2 fields
Geometry type: LINESTRING
Dimension: XY
Bounding box: xmin: 645003 ymin: 4822568 xmax: 748805.6 ymax: 4992001
Projected CRS: WGS 84 / UTM zone 32N
First 3 features:
PK_UID Name Geometry
1 1 Unknown LINESTRING (668540.7 485826...
2 2 Unknown LINESTRING (671365.9 485417...
3 3 Unknown LINESTRING (671417.6 485417...
Simple feature collection with 8101 features and 6 fields
Geometry type: POINT
Dimension: XY
Bounding box: xmin: 319224 ymin: 3934674 xmax: 1308585 ymax: 5214373
Projected CRS: WGS 84 / UTM zone 32N
First 3 features:
PK_UID Name Peoples LocalCounc County Region
1 1 Brozolo 435 1 0 0
2 2 Campiglione-Fenile 1284 1 0 0
3 3 Canischio 274 1 0 0
Geometry
1 POINT (427002.8 4996361)
2 POINT (367470.5 4962414)
3 POINT (390084.1 5025552)
Simple feature collection with 10 features and 2 fields
Geometry type: MULTIPOLYGON
Dimension: XY
Bounding box: xmin: 378359.8 ymin: 4198234 xmax: 1211486 ymax: 5036803
Projected CRS: WGS 84 / UTM zone 32N
First 3 features:
PK_UID Name Geometry
1 1 VENETO MULTIPOLYGON (((752912.3 50...
2 2 VENETO MULTIPOLYGON (((751705.5 50...
3 3 VENETO MULTIPOLYGON (((766247.8 50...
>
> csv = system.file("csv/pt.csv", package = "sf")
> identical(st_read(quiet = TRUE, csv, options = "AUTODETECT_TYPE=Yes")$Int64[3], NA_real_)
[1] TRUE
> identical(st_read(quiet = TRUE, csv, int64_as_string = TRUE, stringsAsFactors = FALSE,
+ options = "AUTODETECT_TYPE=Yes")$Int64[3], NA_character_)
[1] TRUE
> identical(st_read(quiet = TRUE, csv, options = "AUTODETECT_TYPE=Yes")$Int32[3], NA_integer_)
[1] TRUE
>
> if ("GML" %in% st_drivers()$name) {
+ gml = system.file("gml/fmi_test.gml", package = "sf")
+ print(dim(st_read(gml, quiet = TRUE)))
+ gml = system.file("gml/20170930_OB_530964_UKSH.xml.gz", package = "sf")
+ print(dim(st_read(gml, layer = "Parcely", quiet = TRUE)))
+ print(dim(st_read(gml, layer = "Parcely", int64_as_string=TRUE, quiet = TRUE)))
+ }
[1] 22 11
[1] 644 25
[1] 644 25
>
> # logical:
> if ("GPKG" %in% st_drivers()$name) {
+ tst = st_read(system.file("gpkg/nc.gpkg", package="sf"), quiet = TRUE) # default layer name
+ tst$bool = tst$NWBIR79 > 800 # logical
+ tst$bool[1:3] = NA
+ st_write(tst, "tst__.gpkg")
+ tst2 = st_read("tst__.gpkg")
+ print(identical(tst$bool, tst2$bool))
+ }
Writing layer `tst__' to data source `tst__.gpkg' using driver `GPKG'
Writing 100 features with 15 fields and geometry type Multi Polygon.
Reading layer `tst__' from data source
`/home/edzer/git/sf.Rcheck/tests/tst__.gpkg' using driver `GPKG'
Simple feature collection with 100 features and 15 fields
Geometry type: MULTIPOLYGON
Dimension: XY
Bounding box: xmin: -84.32385 ymin: 33.88199 xmax: -75.45698 ymax: 36.58965
Geodetic CRS: NAD27
[1] TRUE
>
> # spatial filter:
> nc = read_sf(system.file("gpkg/nc.gpkg", package="sf"))
> wkt = st_as_text(st_geometry(nc[1,]))
> wkt
[1] "MULTIPOLYGON (((-81.47276 36.23436, -81.54084 36.27251, -81.56198 36.27359, -81.63306 36.34069, -81.74107 36.39178, -81.69828 36.47178, -81.7028 36.51934, -81.67 36.58965, -81.3453 36.57286, -81.34754 36.53791, -81.32478 36.51368, -81.31332 36.4807, -81.26624 36.43721, -81.26284 36.40504, -81.24069 36.37942, -81.23989 36.36536, -81.26424 36.35241, -81.32899 36.3635, -81.36137 36.35316, -81.36569 36.33905, -81.35413 36.29972, -81.36745 36.2787, -81.40639 36.28505, -81.41233 36.26729, -81.43104 36.26072, -81.45289 36.23959, -81.47276 36.23436)))"
> nc_filtered = read_sf(system.file("gpkg/nc.gpkg", package="sf"), wkt_filter = wkt)
> try(read_sf(system.file("gpkg/nc.gpkg", package="sf"), wkt_filter = "wrong"))
Cannot create geometry from: wrong
Error in CPL_read_ogr(dsn, layer, query, as.character(options), quiet, :
wkt parse error.
> nc_filtered$NAME
[1] "Watauga" "Ashe" "Wilkes" "Alleghany"
>
> proc.time()
user system elapsed
1.409 0.087 1.500
Generated by dwww version 1.15 on Sun Jun 30 09:52:51 CEST 2024.