Diff to HTML by rtfpessoa

Files changed (1) show
{orig → migrated}/cm2_mod.R RENAMED
@@ -20,13 +20,12 @@
20
20
  ###################################################
21
21
  ### code chunk number 10: cm2.Rnw:269-270
22
22
  ###################################################
23
- library(rgeos)
24
-
23
+ library(sp)
24
+ library(sf)
25
25
 
26
26
  ###################################################
27
27
  ### code chunk number 13: cm2.Rnw:321-322
28
28
  ###################################################
29
- getScale()
30
29
 
31
30
 
32
31
  ###################################################
@@ -41,28 +40,31 @@
41
40
  ###################################################
42
41
  ### code chunk number 15: cm2.Rnw:400-403
43
42
  ###################################################
44
- library(rgdal)
45
- set_ReplCRS_warn(FALSE)
46
43
 
47
44
 
48
45
  ###################################################
49
46
  ### code chunk number 16: cm2.Rnw:406-409
50
47
  ###################################################
51
- olinda <- readOGR("olinda1.shp", "olinda1", integer64="allow.loss")
52
- proj4string(olinda) <- CRS("+init=epsg:4674")
53
- olinda_utm <- spTransform(olinda, CRS("+init=epsg:31985"))
48
+ olinda_sf <- st_read("olinda1.shp")
49
+ st_precision(olinda_sf)
50
+ olinda <- as(olinda_sf, "Spatial")
51
+ proj4string(olinda) <- CRS("EPSG:4674")
52
+ olinda_utm <- as(st_transform(st_as_sf(olinda), st_crs(CRS("EPSG:31985"))), "Spatial")
53
+ olinda_utm_sf <- st_as_sf(olinda_utm)
54
+ st_precision(olinda_utm_sf)
55
+ st_precision(olinda_utm_sf) <- 1e8
56
+ st_precision(olinda_utm_sf)
54
57
 
55
58
 
56
59
  ###################################################
57
60
  ### code chunk number 17: cm2.Rnw:412-413
58
61
  ###################################################
59
- set_ReplCRS_warn(TRUE)
60
62
 
61
63
 
62
64
  ###################################################
63
65
  ### code chunk number 18: cm2.Rnw:457-461
64
66
  ###################################################
65
- Area <- gArea(olinda_utm, byid=TRUE)
67
+ Area <- st_area(olinda_utm_sf)
66
68
  olinda_utm$area <- sapply(slot(olinda_utm, "polygons"), slot, "area")
67
69
  all.equal(unname(Area), olinda_utm$area)
68
70
  olinda_utm$dens <- olinda_utm$V014/(olinda_utm$area/1000000)
@@ -78,53 +80,60 @@
78
80
  ###################################################
79
81
  ### code chunk number 20: cm2.Rnw:518-521
80
82
  ###################################################
81
- bounds <- gUnaryUnion(olinda_utm)
82
- gArea(olinda_utm)
83
+ bounds <- as(st_union(olinda_utm_sf), "Spatial")
84
+ st_area(st_union(olinda_utm_sf))
83
85
  sapply(slot(slot(bounds, "polygons")[[1]], "Polygons"), slot, "area")
84
86
 
85
87
 
86
88
  ###################################################
87
89
  ### code chunk number 21: cm2.Rnw:542-544
88
90
  ###################################################
89
- pols_overlap <- gOverlaps(olinda_utm, byid=TRUE)
91
+ pols_overlap <- st_overlaps(olinda_utm_sf, sparse=FALSE)
90
92
  any(pols_overlap)
91
93
 
92
94
 
93
95
  ###################################################
94
96
  ### code chunk number 22: cm2.Rnw:566-573
95
97
  ###################################################
96
- oScale <- getScale()
97
- setScale(1e+4)
98
- pols_overlap <- gOverlaps(olinda_utm, byid=TRUE)
98
+ oScale <- st_precision(olinda_utm_sf)
99
+ st_precision(olinda_utm_sf) <- 1e4
100
+ pols_overlap <- st_overlaps(olinda_utm_sf, sparse=FALSE)
99
101
  any(pols_overlap)
100
- bounds <- gUnaryUnion(olinda_utm)
101
- setScale(oScale)
102
+ bounds_sf <- st_union(olinda_utm_sf)
103
+ st_precision(bounds_sf) <- oScale
104
+ bounds <- as(bounds_sf, "Spatial")
105
+ st_precision(olinda_utm_sf) <- oScale
102
106
  sapply(slot(slot(bounds, "polygons")[[1]], "Polygons"), slot, "area")
103
107
 
104
108
 
105
109
  ###################################################
106
110
  ### code chunk number 23: cm2.Rnw:592-593
107
111
  ###################################################
108
- set_ReplCRS_warn(FALSE)
109
112
 
110
113
 
111
114
  ###################################################
112
115
  ### code chunk number 24: cm2.Rnw:596-604
113
116
  ###################################################
114
- pan <- readGDAL("L7_ETM8s.tif")
115
- proj4string(pan) <- CRS(proj4string(bounds))
116
- TM <- readGDAL("L7_ETMs.tif")
117
- proj4string(TM) <- CRS(proj4string(bounds))
117
+ library(stars)
118
+ pan <- as(read_stars("L7_ETM8s.tif"), "Spatial")
119
+ slot(pan, "proj4string") <- slot(bounds, "proj4string")
120
+ TM0 <- read_stars("L7_ETMs.tif")
121
+ for (i in 1:dim(TM0)[3]) {
122
+ TMi <- as(TM0[,,,i, drop=TRUE], "Spatial")
123
+ if (i == 1) TM <- TMi
124
+ else TM <- cbind(TM, TMi)
125
+ }
126
+ slot(TM, "proj4string") <- slot(bounds, "proj4string")
118
127
  names(TM) <- c("TM1", "TM2", "TM3", "TM4", "TM5", "TM7")
119
- dem <- readGDAL("olinda_dem_utm25s.tif")
120
- proj4string(dem) <- CRS(proj4string(bounds))
128
+ dem <- as(read_stars("olinda_dem_utm25s.tif"), "Spatial")
129
+ slot(dem, "proj4string") <- slot(bounds, "proj4string")
130
+ names(dem) <- "band1"
121
131
  is.na(dem$band1) <- dem$band1 <= 0
122
132
 
123
133
 
124
134
  ###################################################
125
135
  ### code chunk number 25: cm2.Rnw:607-609
126
136
  ###################################################
127
- set_ReplCRS_warn(TRUE)
128
137
 
129
138
 
130
139
  ###################################################
@@ -160,7 +169,7 @@
160
169
  ###################################################
161
170
  ### code chunk number 29: cm2.Rnw:691-692
162
171
  ###################################################
163
- set_ReplCRS_warn(FALSE)
172
+ #set_ReplCRS_warn(FALSE)
164
173
 
165
174
 
166
175
  ###################################################
@@ -172,32 +181,31 @@
172
181
  ###################################################
173
182
  ### code chunk number 31: cm2.Rnw:698-699
174
183
  ###################################################
175
- stream_utm <- readOGR("stream.shp", "stream")
176
-
184
+ stream_utm <- as(st_read("stream.shp"), "Spatial")
177
185
 
178
186
  ###################################################
179
187
  ### code chunk number 32: cm2.Rnw:701-702
180
188
  ###################################################
181
- proj4string(stream_utm) <- CRS("+init=epsg:31985")
189
+ slot(stream_utm, "proj4string") <- CRS("EPSG:31985")
182
190
 
183
191
 
184
192
  ###################################################
185
193
  ### code chunk number 33: cm2.Rnw:704-705
186
194
  ###################################################
187
- set_ReplCRS_warn(TRUE)
188
195
 
189
196
 
190
197
  ###################################################
191
198
  ### code chunk number 34: cm2.Rnw:707-709
192
199
  ###################################################
193
200
  nrow(stream_utm)
194
- summary(gLength(stream_utm, byid=TRUE))
201
+ stream_utm_sf <- st_as_sf(stream_utm)
202
+ summary(st_length(stream_utm_sf))
195
203
 
196
204
 
197
205
  ###################################################
198
206
  ### code chunk number 35: cm2.Rnw:726-728
199
207
  ###################################################
200
- t0 <- gTouches(stream_utm, byid=TRUE)
208
+ t0 <- st_touches(stream_utm_sf, sparse=FALSE)
201
209
  any(t0)
202
210
 
203
211
 
@@ -218,16 +226,18 @@
218
226
  ###################################################
219
227
  ### code chunk number 38: cm2.Rnw:768-772
220
228
  ###################################################
221
- lns <- gLineMerge(stream_utm, id=as.character(nComp$comp.id))
229
+ lns <- as(aggregate(stream_utm_sf, list(as.character(nComp$comp.id)), head, n=1), "Spatial")
222
230
  length(row.names(lns))
223
- summary(gLength(lns, byid=TRUE))
224
- all.equal(SpatialLinesLengths(lns), unname(gLength(lns, byid=TRUE)))
225
-
231
+ lns_sf <- st_as_sf(lns)
232
+ st_precision(lns_sf) <- 1e8
233
+ lens <- units::set_units(st_length(lns_sf), NULL)
234
+ summary(lens)
235
+ all.equal(SpatialLinesLengths(lns), lens)
226
236
 
227
237
  ###################################################
228
238
  ### code chunk number 39: cm2.Rnw:791-794
229
239
  ###################################################
230
- GI <- gIntersection(lns, olinda_utm, byid=TRUE)
240
+ GI <- as(st_intersection(lns_sf, olinda_utm_sf), "Spatial")
231
241
  class(GI)
232
242
  length(row.names(GI))
233
243
 
@@ -235,22 +245,12 @@
235
245
  ###################################################
236
246
  ### code chunk number 40: cm2.Rnw:817-827
237
247
  ###################################################
238
- res <- numeric(nrow(olinda_utm))
239
- head(row.names(GI))
240
- range(as.integer(row.names(olinda_utm)))
241
- rnGI <- as.integer(sapply(strsplit(row.names(GI), " "), "[", 2))+1
242
- length(rnGI) == length(unique(rnGI))
243
- lens <- gLength(GI, byid=TRUE)
244
- tlens <- tapply(lens, rnGI, sum)
245
- res[as.integer(names(tlens))] <- unname(tlens)
246
- olinda_utm$stream_len <- res
247
- summary(olinda_utm$stream_len)
248
248
 
249
249
 
250
250
  ###################################################
251
251
  ### code chunk number 41: cm2.Rnw:846-848
252
252
  ###################################################
253
- tree <- gBinarySTRtreeQuery(lns, olinda_utm)
253
+ tree <- st_intersects(olinda_utm_sf, lns_sf)
254
254
  table(sapply(tree, length))
255
255
 
256
256
 
@@ -259,30 +259,30 @@
259
259
  ###################################################
260
260
  res1 <- numeric(length=length(tree))
261
261
  for (i in seq(along=res1)) {
262
- if (!is.null(tree[[i]])) {
263
- gi <- gIntersection(lns[tree[[i]]], olinda_utm[i,])
264
- res1[i] <- ifelse(is.null(gi), 0, gLength(gi))
262
+ if (!(length(tree[[i]]) == 0L)) {
263
+ gi <- st_intersection(olinda_utm_sf[i,], lns_sf[tree[[i]],])
264
+ res1[i] <- units::set_units(st_length(gi), NULL)
265
265
  }
266
266
  }
267
- all.equal(olinda_utm$stream_len, res1)
267
+ olinda_utm$stream_len <- res1
268
268
 
269
269
 
270
270
  ###################################################
271
271
  ### code chunk number 43: cm2.Rnw:899-901
272
272
  ###################################################
273
- buf50m <- gBuffer(lns, width=50)
273
+ buf50m_sf <- st_union(st_buffer(lns_sf, dist=50))
274
+ st_precision(buf50m_sf) <- 1e8
275
+ buf50m <- as(buf50m_sf, "Spatial")
274
276
  length(slot(buf50m, "polygons"))
275
277
 
276
278
 
277
279
  ###################################################
278
280
  ### code chunk number 44: cm2.Rnw:918-926
279
281
  ###################################################
280
- GI1 <- gIntersection(buf50m, olinda_utm, byid=TRUE)
282
+ GI1_sf <- st_intersection(olinda_utm_sf, buf50m_sf)
283
+ GI1_sf$area <- units::set_units(st_area(GI1_sf), NULL)
281
284
  res <- numeric(length(slot(olinda_utm, "polygons")))
282
- head(row.names(GI1))
283
- rnGI <- as.integer(sapply(strsplit(row.names(GI1), " "), "[", 2))+1
284
- length(rnGI) == length(unique(rnGI))
285
- res[rnGI] <- gArea(GI1, byid=TRUE)
285
+ res[match(GI1_sf$ID, olinda_utm$ID)] <- GI1_sf$area
286
286
  olinda_utm$buf_area <- res
287
287
  olinda_utm$prop_50m <- olinda_utm$buf_area/olinda_utm$area
288
288
 
@@ -290,8 +290,7 @@
290
290
  ###################################################
291
291
  ### code chunk number 45: cm2.Rnw:942-943
292
292
  ###################################################
293
- stream_inside <- gIntersection(lns, bounds)
294
-
293
+ stream_inside <- as(st_intersection(lns_sf, bounds_sf), "Spatial")
295
294
 
296
295
  ###################################################
297
296
  ### code chunk number 46: cm2.Rnw:960-970
@@ -367,11 +366,10 @@
367
366
  ###################################################
368
367
  ### code chunk number 57: cm2.Rnw:1226-1230
369
368
  ###################################################
369
+ slot(TM1, "proj4string") <- slot(olinda_utm, "proj4string")
370
370
  o_mean <- over(olinda_utm, TM1[,c("PC1", "PC2", "PC3")])
371
371
  str(o_mean)
372
- row.names(o_mean) <- row.names(olinda_utm)
373
- library(maptools)
374
- olinda_utmA <- spCbind(olinda_utm, o_mean)
372
+ olinda_utmA <- cbind(olinda_utm, o_mean)
375
373
 
376
374
 
377
375
 
@@ -379,9 +377,8 @@
379
377
  ### code chunk number 58: cm2.Rnw:1241-1248
380
378
  ###################################################
381
379
  o_median <- over(olinda_utm, TM1[,c("PC1", "PC2", "PC3")], fn=median)
382
- row.names(o_median) <- row.names(olinda_utmA)
383
380
  names(o_median) <- paste(names(o_median), "med", sep="_")
384
- olinda_utmB <- spCbind(olinda_utmA, o_median)
381
+ olinda_utmB <- cbind(olinda_utmA, o_median)
385
382
  TM1$count <- 1
386
383
  o_count <- over(olinda_utm, TM1[,"count"], fn=sum)
387
384
  olinda_utmB$count <- o_count$count
@@ -396,6 +393,7 @@
396
393
  ###################################################
397
394
  ### code chunk number 60: cm2.Rnw:1272-1277
398
395
  ###################################################
396
+ slot(dem, "proj4string") <- slot(olinda_utm, "proj4string")
399
397
  o_dem_median <- over(olinda_utm, dem, fn=median)
400
398
  olinda_utmB$dem_median <- o_dem_median$band1
401
399
  summary(olinda_utmB$dem_median)
@@ -445,6 +443,7 @@
445
443
  set.seed(9876)
446
444
  p_r <- spsample(bounds, 1000, type="random")
447
445
  length(p_r)
446
+ slot(bounds, "proj4string") <- slot(dem, "proj4string")
448
447
  dem <- dem[bounds,]
449
448
  dem_sp <- as(dem, "SpatialPixelsDataFrame")
450
449
  g_r <- spsample(dem_sp, 1000, type="random")
@@ -502,20 +501,24 @@
502
501
  ###################################################
503
502
  o_sp <- as(olinda_utm, "SpatialPolygons")
504
503
  whichPoly <- over(p_r, o_sp)
505
- whichPoly1 <- gContains(o_sp, p_r, byid=TRUE)
506
- whichPoly1a <- apply(unname(whichPoly1), 1, which)
504
+ whichPoly1a <- st_within(st_as_sf(p_r), st_as_sf(o_sp))
507
505
  table(sapply(whichPoly1a, length))
508
- all.equal(whichPoly, whichPoly1a)
506
+ all.equal(unname(whichPoly), unlist(whichPoly1a))
509
507
 
510
508
 
511
509
  ###################################################
512
510
  ### code chunk number 74: cm2.Rnw:1542-1547
513
511
  ###################################################
514
512
  hels <- matrix(c(24.97, 60.17), nrow=1)
515
- p4s <- CRS("+proj=longlat +datum=WGS84")
513
+ p4s <- as(st_crs("EPSG:4326"), "CRS")
516
514
  Hels <- SpatialPoints(hels, proj4string=p4s)
517
515
  d041224 <- as.POSIXct("2004-12-24", tz="EET")
518
- sunriset(Hels, d041224, direction="sunrise", POSIXct.out=TRUE)
516
+ solartime::computeSunriseHour(d041224, latDeg=hels[1,2], longDeg=hels[1,1])
517
+
519
518
 
519
+ (sI <- sessionInfo()) # check: no sp?
520
520
 
521
+ "rgdal" %in% c(names(sI$otherPkgs), names(sI$loadedOnly))
522
+ "rgeos" %in% c(names(sI$otherPkgs), names(sI$loadedOnly))
523
+ "maptools" %in% c(names(sI$otherPkgs), names(sI$loadedOnly))
521
524