1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
|
R version 4.0.0 (2020-04-24) -- "Arbor Day"
Copyright (C) 2020 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
1 points Point NA 10
2 lines Line String NA 9
3 multilinestrings Multi Line String NA 4
4 multipolygons Multi Polygon NA 25
5 other_relations Geometry Collection NA 4
Driver: OSM
Available layers:
layer_name geometry_type features fields
1 points Point 0 10
2 lines Line String 0 9
3 multilinestrings Multi Line String 0 4
4 multipolygons Multi Polygon 3 25
5 other_relations Geometry Collection 0 4
Simple feature collection with 3 features and 25 fields
geometry type: MULTIPOLYGON
dimension: XY
bbox: xmin: -1.556254 ymin: 53.79519 xmax: -1.550716 ymax: 53.80341
geographic 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
bbox: xmin: -84.32385 ymin: 33.88199 xmax: -75.45698 ymax: 36.58965
geographic 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 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 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)
+ }
Simple feature collection with 5 features and 4 fields
geometry type: POINT
dimension: XY
bbox: xmin: 1 ymin: 1 xmax: 5 ymax: 5
CRS: NA
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
bbox: 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
bbox: 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
bbox: 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
bbox: 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
bbox: 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
bbox: 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
bbox: 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
bbox: xmin: -84.32385 ymin: 33.88199 xmax: -75.45698 ymax: 36.58965
geographic 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.320 0.047 1.357
|
