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
|
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.
> suppressPackageStartupMessages(library(sf))
>
> tif = system.file("tif/geomatrix.tif", package = "sf")
>
> gdal_metadata(tif)
[1] "AREA_OR_POINT=Point"
> gdal_metadata(tif, NA_character_)
[1] "IMAGE_STRUCTURE" "DERIVED_SUBDATASETS" ""
> try(gdal_metadata(tif, "wrongDomain"))
Error in gdal_metadata(tif, "wrongDomain") :
domain_item[1] not found in available metadata domains
> gdal_metadata(tif, c("IMAGE_STRUCTURE"))
[[1]]
[1] "BAND"
attr(,"class")
[1] "gdal_metadata"
> try(length(gdal_metadata(tif, c("DERIVED_SUBDATASETS")))) # fails on Fedora 26
[1] 2
>
> if (require(stars)) {
+ tif = system.file("tif/geomatrix.tif", package = "sf")
+ r = read_stars(tif)
+ d = (st_dimensions(r))
+ gt = c(1841001.75, 1.5, -5, 1144003.25, -5, -1.5)
+ x1 = st_as_sfc(d, as_points = TRUE, use_cpp = TRUE, geotransform = gt)
+ x2 = st_as_sfc(d, as_points = TRUE, use_cpp = FALSE, geotransform = gt)
+ print(identical(x1, x2))
+ y1 = st_as_sfc(d, as_points = FALSE, use_cpp = TRUE, geotransform = gt)
+ y2 = st_as_sfc(d, as_points = FALSE, use_cpp = FALSE, geotransform = gt)
+ print(identical(y1, y2))
+
+ # rectilinear grid:
+ m = matrix(1:20, nrow = 5, ncol = 4)
+ x = c(0,0.5,1,2,4,5)
+ y = c(0.3,0.5,1,2,2.2)
+ r = st_as_stars(list(m = m), dimensions = st_dimensions(x = x, y = y, .raster = c("x", "y")))
+ print(st_as_sfc(st_dimensions(r), as_points = TRUE))
+ print(st_as_sfc(st_dimensions(r), as_points = FALSE))
+
+ # curvilinear grid:
+ lon = st_as_stars(matrix(1:5, 4, 5, byrow = TRUE))
+ lat = st_as_stars(matrix(1:4, 4, 5))
+ ll = c(X1 = lon, X2 = lat)
+ curv = st_as_stars(st_as_stars(t(m)), curvilinear = setNames(ll, c("X1", "X2")))
+ print(st_as_sfc(st_dimensions(curv), as_points = TRUE))
+ print(st_as_sfc(st_dimensions(curv), as_points = FALSE))
+
+ demo(nc, echo = FALSE, ask = FALSE)
+ print(x <- st_rasterize(nc)) # default grid:
+ print(p <- st_as_sf(x, as_points = FALSE)) # polygonize: follow raster boundaries
+ print(p <- st_as_sf(x, as_points = FALSE, use_integer = TRUE)) # polygonize integers: follow raster boundaries
+ print(try(p <- st_as_sf(x, as_points = TRUE))) # polygonize: contour, requies GDAL >= 2.4.0
+ if (utils::packageVersion("stars") >= "0.2-1") {
+ write_stars(read_stars(tif), tempfile(fileext = ".tif"))
+ write_stars(read_stars(tif, proxy = TRUE), tempfile(fileext = ".tif"))
+ write_stars(read_stars(tif, proxy = TRUE), tempfile(fileext = ".tif"), chunk_size = c(200,200))
+ na.tif = read_stars(system.file("tif/na.tif", package = "stars"))
+ write_stars(na.tif, "na.tif")
+ write_stars(na.tif, "na.tif", NA_value = -999)
+ na.tif = read_stars(system.file("tif/na.tif", package = "stars"), NA_value = -999)
+ write_stars(na.tif, "na.tif")
+ write_stars(na.tif, "na.tif", NA_value = -999)
+ na.tif = read_stars(system.file("tif/na.tif", package = "stars"), NA_value = -999, proxy = TRUE)
+ write_stars(na.tif, "na.tif")
+ write_stars(na.tif, "na.tif", NA_value = -999)
+ }
+ # https://github.com/mtennekes/tmap/issues/368
+ if (utils::packageVersion("stars") > "0.4-0") {
+ lc = system.file('tif/lc.tif', package = 'stars')
+ if (lc != "") {
+ r = read_stars(lc, RAT = "Land Cover Class")
+ r <- droplevels(r)
+ }
+ }
+ }
Loading required package: stars
Loading required package: abind
[1] TRUE
[1] TRUE
Geometry set for 20 features
geometry type: POINT
dimension: XY
bbox: xmin: 0.25 ymin: 0.4 xmax: 4.5 ymax: 2.1
CRS: NA
First 5 geometries:
POINT (0.25 0.4)
POINT (0.75 0.4)
POINT (1.5 0.4)
POINT (3 0.4)
POINT (4.5 0.4)
Geometry set for 20 features
geometry type: POLYGON
dimension: XY
bbox: xmin: 0 ymin: 0.3 xmax: 5 ymax: 2.2
CRS: NA
First 5 geometries:
POLYGON ((0 0.3, 0.5 0.3, 0.5 0.5, 0 0.5, 0 0.3))
POLYGON ((0.5 0.3, 1 0.3, 1 0.5, 0.5 0.5, 0.5 0...
POLYGON ((1 0.3, 2 0.3, 2 0.5, 1 0.5, 1 0.3))
POLYGON ((2 0.3, 4 0.3, 4 0.5, 2 0.5, 2 0.3))
POLYGON ((4 0.3, 5 0.3, 5 0.5, 4 0.5, 4 0.3))
Geometry set for 20 features
geometry type: POINT
dimension: XY
bbox: xmin: 1 ymin: 1 xmax: 5 ymax: 4
geographic CRS: WGS 84
First 5 geometries:
POINT (1 1)
POINT (1 2)
POINT (1 3)
POINT (1 4)
POINT (2 1)
Geometry set for 20 features
geometry type: POLYGON
dimension: XY
bbox: xmin: 0.5 ymin: 0.5 xmax: 5.5 ymax: 4.5
geographic CRS: WGS 84
First 5 geometries:
POLYGON ((0.5 0.5, 0.5 1.5, 1.5 1.5, 1.5 0.5, 0...
POLYGON ((0.5 1.5, 0.5 2.5, 1.5 2.5, 1.5 1.5, 0...
POLYGON ((0.5 2.5, 0.5 3.5, 1.5 3.5, 1.5 2.5, 0...
POLYGON ((0.5 3.5, 0.5 4.5, 1.5 4.5, 1.5 3.5, 0...
POLYGON ((1.5 0.5, 1.5 1.5, 2.5 1.5, 2.5 0.5, 1...
Reading layer `nc.gpkg' from data source `/home/edzer/git/sf.Rcheck/sf/gpkg/nc.gpkg' using driver `GPKG'
Simple feature collection with 100 features and 14 fields
Attribute-geometry relationship: 0 constant, 8 aggregate, 6 identity
geometry type: MULTIPOLYGON
dimension: XY
bbox: xmin: -84.32385 ymin: 33.88199 xmax: -75.45698 ymax: 36.58965
geographic CRS: NAD27
stars object with 2 dimensions and 1 attribute
attribute(s):
AREA
Min. :0.042
1st Qu.:0.108
Median :0.142
Mean :0.145
3rd Qu.:0.181
Max. :0.241
NA's :30904
dimension(s):
from to offset delta refsys point values
x 1 461 -84.3239 0.0192484 NAD27 FALSE NULL [x]
y 1 141 36.5896 -0.0192484 NAD27 FALSE NULL [y]
Simple feature collection with 34097 features and 1 field
geometry type: POLYGON
dimension: XY
bbox: xmin: -84.32385 ymin: 33.87563 xmax: -75.45034 ymax: 36.58965
geographic CRS: NAD27
First 10 features:
AREA geometry
1 0.114 POLYGON ((-81.66757 36.5896...
2 0.114 POLYGON ((-81.64833 36.5896...
3 0.114 POLYGON ((-81.62908 36.5896...
4 0.114 POLYGON ((-81.60983 36.5896...
5 0.114 POLYGON ((-81.59058 36.5896...
6 0.114 POLYGON ((-81.57133 36.5896...
7 0.114 POLYGON ((-81.55208 36.5896...
8 0.114 POLYGON ((-81.53283 36.5896...
9 0.114 POLYGON ((-81.51359 36.5896...
10 0.114 POLYGON ((-81.49434 36.5896...
Simple feature collection with 34097 features and 1 field
geometry type: POLYGON
dimension: XY
bbox: xmin: -84.32385 ymin: 33.87563 xmax: -75.45034 ymax: 36.58965
geographic CRS: NAD27
First 10 features:
AREA geometry
1 0.114 POLYGON ((-81.66757 36.5896...
2 0.114 POLYGON ((-81.64833 36.5896...
3 0.114 POLYGON ((-81.62908 36.5896...
4 0.114 POLYGON ((-81.60983 36.5896...
5 0.114 POLYGON ((-81.59058 36.5896...
6 0.114 POLYGON ((-81.57133 36.5896...
7 0.114 POLYGON ((-81.55208 36.5896...
8 0.114 POLYGON ((-81.53283 36.5896...
9 0.114 POLYGON ((-81.51359 36.5896...
10 0.114 POLYGON ((-81.49434 36.5896...
Simple feature collection with 34097 features and 1 field
geometry type: POINT
dimension: XY
bbox: xmin: -84.31423 ymin: 33.88525 xmax: -75.45997 ymax: 36.58003
geographic CRS: NAD27
First 10 features:
AREA geometry
1 0.114 POINT (-81.65795 36.58003)
2 0.114 POINT (-81.6387 36.58003)
3 0.114 POINT (-81.61945 36.58003)
4 0.114 POINT (-81.6002 36.58003)
5 0.114 POINT (-81.58096 36.58003)
6 0.114 POINT (-81.56171 36.58003)
7 0.114 POINT (-81.54246 36.58003)
8 0.114 POINT (-81.52321 36.58003)
9 0.114 POINT (-81.50396 36.58003)
10 0.114 POINT (-81.48471 36.58003)
>
> r = gdal_read(tif)
> gt = c(0,1,0,0,0,1)
> gdal_inv_geotransform(gt)
[1] 0 1 0 0 0 1
> rc = expand.grid(x=1:3, y = 1:3)
> #(xy = xy_from_colrow(rc, gt))
> #xy_from_colrow(xy, gt, inverse = TRUE)
> crs <- gdal_crs(tif)
>
> try(gdal_metadata("foo"))
[1] NA
> gdal_metadata(tif)
[1] "AREA_OR_POINT=Point"
>
> proc.time()
user system elapsed
1.344 0.103 1.438
|
