forked from r-spatial/stars
-
Notifications
You must be signed in to change notification settings - Fork 0
/
nc.Rout.save
277 lines (272 loc) · 10.8 KB
/
nc.Rout.save
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
R version 4.0.3 (2020-10-10) -- "Bunny-Wunnies Freak Out"
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.
> # examples with the sf::nc polygon dataset NOT ABOUT NETCDF's, see e.g. stars.R netcdf.R
> suppressPackageStartupMessages(library(sf))
> nc = read_sf(system.file("gpkg/nc.gpkg", package="sf"))
> m = st_set_geometry(nc, NULL)
> n = as.matrix(m[c("BIR74", "SID74", "NWBIR74", "BIR79", "SID79", "NWBIR79")])
> dim(n) = c(county = 100, var = 3, year = 2)
> dimnames(n) = list(county = nc$NAME, var = c("BIR", "SID", "NWBIR"), year = c(1974, 1979))
> suppressPackageStartupMessages(library(stars))
> (st = st_as_stars(pop = n))
stars object with 3 dimensions and 1 attribute
attribute(s):
pop
Min. : 0
1st Qu.: 8
Median : 538
Mean : 1657
3rd Qu.: 1784
Max. :30757
dimension(s):
from to offset delta refsys point values
county 1 100 NA NA NA NA Ashe,...,Brunswick
var 1 3 NA NA NA NA BIR , SID , NWBIR
year 1 2 NA NA NA NA 1974, 1979
> foo <- st %>% st_set_dimensions(1, st_geometry(nc)) # %>% st_set_dimensions(3, c(1974, 1979))
> st %>% st_set_dimensions(1, st_geometry(nc)) %>% st_set_dimensions(names = c("geometries", "var", "year"))
stars object with 3 dimensions and 1 attribute
attribute(s):
pop
Min. : 0
1st Qu.: 8
Median : 538
Mean : 1657
3rd Qu.: 1784
Max. :30757
dimension(s):
from to offset delta refsys point
geometries 1 100 NA NA NAD27 FALSE
var 1 3 NA NA NA NA
year 1 2 NA NA NA NA
values
geometries MULTIPOLYGON (((-81.47276 3...,...,MULTIPOLYGON (((-78.65572 3...
var BIR , SID , NWBIR
year 1974, 1979
> foo
stars object with 3 dimensions and 1 attribute
attribute(s):
pop
Min. : 0
1st Qu.: 8
Median : 538
Mean : 1657
3rd Qu.: 1784
Max. :30757
dimension(s):
from to offset delta refsys point
county 1 100 NA NA NAD27 FALSE
var 1 3 NA NA NA NA
year 1 2 NA NA NA NA
values
county MULTIPOLYGON (((-81.47276 3...,...,MULTIPOLYGON (((-78.65572 3...
var BIR , SID , NWBIR
year 1974, 1979
> st_bbox(foo)
xmin ymin xmax ymax
-84.32385 33.88199 -75.45698 36.58965
> (x = st_as_sf(foo))
Simple feature collection with 100 features and 6 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:
pop.BIR.1974 pop.SID.1974 pop.NWBIR.1974 pop.BIR.1979 pop.SID.1979
Ashe 1091 1 10 1364 0
Alleghany 487 0 10 542 3
Surry 3188 5 208 3616 6
Currituck 508 1 123 830 2
Northampton 1421 9 1066 1606 3
Hertford 1452 7 954 1838 5
Camden 286 0 115 350 2
Gates 420 0 254 594 2
Warren 968 4 748 1190 2
Stokes 1612 1 160 2038 5
pop.NWBIR.1979 county
Ashe 19 MULTIPOLYGON (((-81.47276 3...
Alleghany 12 MULTIPOLYGON (((-81.23989 3...
Surry 260 MULTIPOLYGON (((-80.45634 3...
Currituck 145 MULTIPOLYGON (((-76.00897 3...
Northampton 1197 MULTIPOLYGON (((-77.21767 3...
Hertford 1237 MULTIPOLYGON (((-76.74506 3...
Camden 139 MULTIPOLYGON (((-76.00897 3...
Gates 371 MULTIPOLYGON (((-76.56251 3...
Warren 844 MULTIPOLYGON (((-78.30876 3...
Stokes 176 MULTIPOLYGON (((-80.02567 3...
> frac = function(x) x[2] / x[1]
> frac2 = function(x) c(sidsr = x[2] / x[1], nwbr = x[3] / x[1])
> frac2an = function(x) c(x[2] / x[1], x[3] / x[1])
> st_apply(foo, c(1,3), frac)
stars object with 2 dimensions and 1 attribute
attribute(s):
frac
Min. :0.000000
1st Qu.:0.001198
Median :0.001931
Mean :0.002042
3rd Qu.:0.002599
Max. :0.009554
dimension(s):
from to offset delta refsys point
county 1 100 NA NA NAD27 FALSE
year 1 2 NA NA NA NA
values
county MULTIPOLYGON (((-81.47276 3...,...,MULTIPOLYGON (((-78.65572 3...
year 1974, 1979
> st_apply(foo, c(1,3), frac2)
stars object with 3 dimensions and 1 attribute
attribute(s):
pop
Min. :0.000000
1st Qu.:0.001923
Median :0.005120
Mean :0.156758
3rd Qu.:0.306530
Max. :0.772727
dimension(s):
from to offset delta refsys point
frac2 1 2 NA NA NA NA
county 1 100 NA NA NAD27 FALSE
year 1 2 NA NA NA NA
values
frac2 sidsr.SID , nwbr.NWBIR
county MULTIPOLYGON (((-81.47276 3...,...,MULTIPOLYGON (((-78.65572 3...
year 1974, 1979
> st_apply(foo, c(1,3), frac2an)
stars object with 3 dimensions and 1 attribute
attribute(s):
pop
Min. :0.000000
1st Qu.:0.001923
Median :0.005120
Mean :0.156758
3rd Qu.:0.306530
Max. :0.772727
dimension(s):
from to offset delta refsys point
frac2an 1 2 NA NA NA NA
county 1 100 NA NA NAD27 FALSE
year 1 2 NA NA NA NA
values
frac2an SID , NWBIR
county MULTIPOLYGON (((-81.47276 3...,...,MULTIPOLYGON (((-78.65572 3...
year 1974, 1979
> library(abind)
> (x = aperm(st_apply(foo, c(1,3), frac2), c(2,3,1)))
stars object with 3 dimensions and 1 attribute
attribute(s):
pop
Min. :0.000000
1st Qu.:0.001923
Median :0.005120
Mean :0.156758
3rd Qu.:0.306530
Max. :0.772727
dimension(s):
from to offset delta refsys point
county 1 100 NA NA NAD27 FALSE
year 1 2 NA NA NA NA
frac2 1 2 NA NA NA NA
values
county MULTIPOLYGON (((-81.47276 3...,...,MULTIPOLYGON (((-78.65572 3...
year 1974, 1979
frac2 sidsr.SID , nwbr.NWBIR
> y = aperm(st_apply(foo, c(1,3), frac2), c("county","year","frac2"))
> all.equal(st_dimensions(x), st_dimensions(y))
[1] TRUE
>
> split(foo, 2)
stars object with 2 dimensions and 3 attributes
attribute(s):
BIR SID NWBIR
Min. : 248 Min. : 0.000 Min. : 1
1st Qu.: 1177 1st Qu.: 2.000 1st Qu.: 206
Median : 2265 Median : 5.000 Median : 742
Mean : 3762 Mean : 7.515 Mean : 1202
3rd Qu.: 4451 3rd Qu.: 9.000 3rd Qu.: 1316
Max. :30757 Max. :57.000 Max. :11631
dimension(s):
from to offset delta refsys point
county 1 100 NA NA NAD27 FALSE
year 1 2 NA NA NA NA
values
county MULTIPOLYGON (((-81.47276 3...,...,MULTIPOLYGON (((-78.65572 3...
year 1974, 1979
> split(foo, 3)
stars object with 2 dimensions and 2 attributes
attribute(s):
1974 1979
Min. : 0.0 Min. : 0
1st Qu.: 8.0 1st Qu.: 9
Median : 501.5 Median : 614
Mean : 1452.4 Mean : 1862
3rd Qu.: 1580.5 3rd Qu.: 2040
Max. :21588.0 Max. :30757
dimension(s):
from to offset delta refsys point
county 1 100 NA NA NAD27 FALSE
var 1 3 NA NA NA NA
values
county MULTIPOLYGON (((-81.47276 3...,...,MULTIPOLYGON (((-78.65572 3...
var BIR , SID , NWBIR
>
> st_crs(foo)
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["geodetic latitude (Lat)",north,
ORDER[1],
ANGLEUNIT["degree",0.0174532925199433]],
AXIS["geodetic longitude (Lon)",east,
ORDER[2],
ANGLEUNIT["degree",0.0174532925199433]],
USAGE[
SCOPE["unknown"],
AREA["North America - NAD27"],
BBOX[7.15,167.65,83.17,-47.74]],
ID["EPSG",4267]]
> plot(foo)
>
> # subset vector cube:
> foo[nc[1]]
although coordinates are longitude/latitude, st_intersects assumes that they are planar
stars object with 3 dimensions and 1 attribute
attribute(s):
pop
Min. : 0
1st Qu.: 8
Median : 538
Mean : 1657
3rd Qu.: 1784
Max. :30757
dimension(s):
from to offset delta refsys point
county 1 100 NA NA NAD27 FALSE
var 1 3 NA NA NA NA
year 1 2 NA NA NA NA
values
county MULTIPOLYGON (((-81.47276 3...,...,MULTIPOLYGON (((-78.65572 3...
var BIR , SID , NWBIR
year 1974, 1979
>
> proc.time()
user system elapsed
1.128 0.208 1.033