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
|
R version 3.6.2 (2019-12-12) -- "Dark and Stormy Night"
Copyright (C) 2019 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))
> round_trip = function(x, EWKB = FALSE, pureR = FALSE) {
+ if (inherits(x, "sfg"))
+ x = st_sfc(x)
+ wkb = st_as_binary(x, EWKB = EWKB, pureR = pureR)
+ class(wkb) = "WKB"
+ # print(wkb)
+ y = st_as_sfc(wkb, EWKB = EWKB, pureR = pureR)
+ a = all.equal(x, y)
+ if (length(a) == 1 && is.logical(a) && a)
+ TRUE
+ else {
+ print(x)
+ print(wkb)
+ print(y)
+ FALSE
+ }
+ }
>
> p3 = st_point(c(0,0,0))
> p3m = st_point(c(0,0,0), "XYM")
> p4 = st_point(c(0,0,0,0))
> p2 = st_point(c(0,0))
> ls = st_linestring(matrix(1:6,3))
> mp = st_multipoint(matrix(1:6,3))
>
> 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)
>
> pol1 = list(outer, hole1, hole2)
> pol2 = list(outer + 12, hole1 + 12)
> pol3 = list(outer + 24)
> mp1 = st_multipolygon(list(pol1,pol2,pol3))
>
> ml1 = st_multilinestring(list(outer, hole1, hole2))
> gc = st_geometrycollection(list(p2, ls, pl1, mp1))
>
> sapply(list(p3, p3m, p4, p2, ls, mp, pl1, mp1, ml1, gc), round_trip, EWKB = FALSE)
[1] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
> sapply(list(p3, p3m, p4, p2, ls, mp, pl1, mp1, ml1, gc), round_trip, EWKB = FALSE, pureR = TRUE)
[1] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
> sapply(list(p3, p3m, p4, p2, ls, mp, pl1, mp1, ml1, gc), round_trip, EWKB = TRUE)
[1] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
> sapply(list(p3, p3m, p4, p2, ls, mp, pl1, mp1, ml1, gc), round_trip, EWKB = TRUE, pureR = TRUE)
[1] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
>
> rawToHex(st_as_binary(st_multipoint(matrix(1:6,3))))
[1] "0104000000030000000101000000000000000000f03f0000000000001040010100000000000000000000400000000000001440010100000000000000000008400000000000001840"
> rawToHex(st_as_binary(st_sfc(st_point(c(0,1)), st_multipoint(matrix(1:6,3)))))
[1] "01010000000000000000000000000000000000f03f"
[2] "0104000000030000000101000000000000000000f03f0000000000001040010100000000000000000000400000000000001440010100000000000000000008400000000000001840"
> try(rawToHex("error"))
Error in rawToHex("error") :
not implemented for objects of class character
>
> # debug roundtrips sf -> GDAL -> sf;
> # the first WKT is what GDAL reports, and will lack M
> st_as_text(st_sfc(sf:::CPL_roundtrip(st_sfc(st_linestring(matrix(1:18,6,3),dim="XYZ")))))
LINESTRING (1 7 13,2 8 14,3 9 15,4 10 16,5 11 17,6 12 18)
[1] "LINESTRING Z (1 7 13, 2 8 14, 3 9 15, 4 10 16, 5 11 17, 6 12 18)"
> st_as_text(st_sfc(sf:::CPL_roundtrip(st_sfc(st_multipoint(matrix(1:18,6,3),dim="XYZ")))))
MULTIPOINT (1 7 13,2 8 14,3 9 15,4 10 16,5 11 17,6 12 18)
[1] "MULTIPOINT Z ((1 7 13), (2 8 14), (3 9 15), (4 10 16), (5 11 17), (6 12 18))"
> st_as_text(st_sfc(sf:::CPL_roundtrip(st_sfc(st_point(c(0,0,0), dim="XYZ")))))
POINT (0 0 0)
[1] "POINT Z (0 0 0)"
>
> if (sf:::CPL_gdal_version() >= "2.1.0") { # address GDAL/Fedora (gdal 2.0.2) error:
+ st_as_text(st_sfc(sf:::CPL_roundtrip(st_sfc(st_linestring(matrix(1:18,6,3),dim="XYM")))))
+ st_as_text(st_sfc(sf:::CPL_roundtrip(st_sfc(st_multipoint(matrix(1:18,6,3),dim="XYM")))))
+ st_as_text(st_sfc(sf:::CPL_roundtrip(st_sfc(st_point(c(0,0,0), dim="XYM")))))
+ } else {
+ "(output expected when gdal <= 2.1.0, e.g. CRAN/fedora)"
+ }
LINESTRING (1 7,2 8,3 9,4 10,5 11,6 12)
MULTIPOINT (1 7,2 8,3 9,4 10,5 11,6 12)
POINT (0 0)
[1] "POINT M (0 0 0)"
>
> proc.time()
user system elapsed
0.505 0.036 0.532
|
