Diff to HTML by rtfpessoa

Files changed (1) show
{orig → migrated_terra}/cm2_mod.R RENAMED
@@ -20,13 +20,13 @@
20
20
  ###################################################
21
21
  ### code chunk number 10: cm2.Rnw:269-270
22
22
  ###################################################
23
- library(rgeos)
23
+ library(sp)
24
+ library(terra)
24
25
 
25
26
 
26
27
  ###################################################
27
28
  ### code chunk number 13: cm2.Rnw:321-322
28
29
  ###################################################
29
- getScale()
30
30
 
31
31
 
32
32
  ###################################################
@@ -41,29 +41,27 @@
41
41
  ###################################################
42
42
  ### code chunk number 15: cm2.Rnw:400-403
43
43
  ###################################################
44
- library(rgdal)
45
- set_ReplCRS_warn(FALSE)
46
44
 
47
45
 
48
46
  ###################################################
49
47
  ### code chunk number 16: cm2.Rnw:406-409
50
48
  ###################################################
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"))
49
+ olinda <- vect("olinda1.shp")
50
+ crs(olinda) <- "EPSG:4674"
51
+ olinda_utm <- project(olinda, "EPSG:31985")
54
52
 
55
53
 
56
54
  ###################################################
57
55
  ### code chunk number 17: cm2.Rnw:412-413
58
56
  ###################################################
59
- set_ReplCRS_warn(TRUE)
60
57
 
61
58
 
62
59
  ###################################################
63
60
  ### code chunk number 18: cm2.Rnw:457-461
64
61
  ###################################################
65
- Area <- gArea(olinda_utm, byid=TRUE)
66
- olinda_utm$area <- sapply(slot(olinda_utm, "polygons"), slot, "area")
62
+ library(raster)
63
+ Area <- area(as(olinda_utm, "Spatial"))
64
+ olinda_utm$area <- sapply(slot(as(olinda_utm, "Spatial"), "polygons"), slot, "area")
67
65
  all.equal(unname(Area), olinda_utm$area)
68
66
  olinda_utm$dens <- olinda_utm$V014/(olinda_utm$area/1000000)
69
67
 
@@ -72,59 +70,55 @@
72
70
  ### code chunk number 19: cm2.Rnw:480-490
73
71
  ###################################################
74
72
  library(RColorBrewer)
75
- spplot(olinda_utm, "dens", at=c(0, 8000, 12000, 15000, 20000, 60000), col.regions=brewer.pal(6, "YlOrBr")[-1], col="grey30", lwd=0.5, colorkey=list(space="right", labels=list(cex=0.7), width=1))
73
+ spplot(as(olinda_utm, "Spatial"), "dens", at=c(0, 8000, 12000, 15000, 20000, 60000), col.regions=brewer.pal(6, "YlOrBr")[-1], col="grey30", lwd=0.5, colorkey=list(space="right", labels=list(cex=0.7), width=1))
76
74
 
77
75
 
78
76
  ###################################################
79
77
  ### code chunk number 20: cm2.Rnw:518-521
80
78
  ###################################################
81
- bounds <- gUnaryUnion(olinda_utm)
82
- gArea(olinda_utm)
83
- sapply(slot(slot(bounds, "polygons")[[1]], "Polygons"), slot, "area")
79
+ bounds <- aggregate(olinda_utm, dissolve=TRUE)
80
+ area(as(bounds, "Spatial"))
81
+ sapply(slot(slot(as(bounds, "Spatial"), "polygons")[[1]], "Polygons"), slot, "area")
84
82
 
85
83
 
86
84
  ###################################################
87
85
  ### code chunk number 21: cm2.Rnw:542-544
88
86
  ###################################################
89
- pols_overlap <- gOverlaps(olinda_utm, byid=TRUE)
90
- any(pols_overlap)
91
87
 
92
88
 
93
89
  ###################################################
94
90
  ### code chunk number 22: cm2.Rnw:566-573
95
91
  ###################################################
96
- oScale <- getScale()
97
- setScale(1e+4)
98
- pols_overlap <- gOverlaps(olinda_utm, byid=TRUE)
99
- any(pols_overlap)
100
- bounds <- gUnaryUnion(olinda_utm)
101
- setScale(oScale)
102
- sapply(slot(slot(bounds, "polygons")[[1]], "Polygons"), slot, "area")
103
92
 
104
93
 
105
94
  ###################################################
106
95
  ### code chunk number 23: cm2.Rnw:592-593
107
96
  ###################################################
108
- set_ReplCRS_warn(FALSE)
109
97
 
110
98
 
111
99
  ###################################################
112
100
  ### code chunk number 24: cm2.Rnw:596-604
113
101
  ###################################################
114
- pan <- readGDAL("L7_ETM8s.tif")
115
- proj4string(pan) <- CRS(proj4string(bounds))
116
- TM <- readGDAL("L7_ETMs.tif")
117
- proj4string(TM) <- CRS(proj4string(bounds))
102
+ pan <- rast("L7_ETM8s.tif")
103
+ crs(pan) <- crs(bounds)
104
+ TM <- rast("L7_ETMs.tif")
105
+ crs(TM) <- crs(bounds)
118
106
  names(TM) <- c("TM1", "TM2", "TM3", "TM4", "TM5", "TM7")
119
- dem <- readGDAL("olinda_dem_utm25s.tif")
120
- proj4string(dem) <- CRS(proj4string(bounds))
121
- is.na(dem$band1) <- dem$band1 <= 0
107
+ dem <- rast("olinda_dem_utm25s.tif")
108
+ crs(dem) <- crs(bounds)
109
+ #RH unnecessarily complicated and inefficient
110
+ # values(dem)[values(dem) <= 0] <- NA
111
+ # better
112
+ dem[dem <= 0] <- NA
113
+ # what I would do (in order of preference)
114
+ dem <- clamp(dem, 0, Inf, values=FALSE)
115
+ dem <- classify(dem, cbind(-Inf, 0, NA))
116
+ dem <- ifel(dem < 0, NA, dem)
122
117
 
123
118
 
124
119
  ###################################################
125
120
  ### code chunk number 25: cm2.Rnw:607-609
126
121
  ###################################################
127
- set_ReplCRS_warn(TRUE)
128
122
 
129
123
 
130
124
  ###################################################
@@ -160,7 +154,6 @@
160
154
  ###################################################
161
155
  ### code chunk number 29: cm2.Rnw:691-692
162
156
  ###################################################
163
- set_ReplCRS_warn(FALSE)
164
157
 
165
158
 
166
159
  ###################################################
@@ -172,33 +165,32 @@
172
165
  ###################################################
173
166
  ### code chunk number 31: cm2.Rnw:698-699
174
167
  ###################################################
175
- stream_utm <- readOGR("stream.shp", "stream")
168
+ stream_utm <- vect("stream.shp", "stream")
176
169
 
177
170
 
178
171
  ###################################################
179
172
  ### code chunk number 32: cm2.Rnw:701-702
180
173
  ###################################################
181
- proj4string(stream_utm) <- CRS("+init=epsg:31985")
174
+ crs(stream_utm) <- "EPSG:31985"
182
175
 
183
176
 
184
177
  ###################################################
185
178
  ### code chunk number 33: cm2.Rnw:704-705
186
179
  ###################################################
187
- set_ReplCRS_warn(TRUE)
188
180
 
189
181
 
190
182
  ###################################################
191
183
  ### code chunk number 34: cm2.Rnw:707-709
192
184
  ###################################################
193
185
  nrow(stream_utm)
194
- summary(gLength(stream_utm, byid=TRUE))
186
+ summary(SpatialLinesLengths(as(stream_utm, "Spatial")))
195
187
 
196
188
 
197
189
  ###################################################
198
190
  ### code chunk number 35: cm2.Rnw:726-728
199
191
  ###################################################
200
- t0 <- gTouches(stream_utm, byid=TRUE)
201
- any(t0)
192
+ t0 <- relate(stream_utm, relation="touches", pairs=FALSE)
193
+ any(t0 > 0)
202
194
 
203
195
 
204
196
  ###################################################
@@ -218,79 +210,63 @@
218
210
  ###################################################
219
211
  ### code chunk number 38: cm2.Rnw:768-772
220
212
  ###################################################
221
- lns <- gLineMerge(stream_utm, id=as.character(nComp$comp.id))
213
+ stream_utm$comp_id <- nComp$comp.id
214
+ lns <- aggregate(stream_utm, by="comp_id", fun=head, n=1)
222
215
  length(row.names(lns))
223
- summary(gLength(lns, byid=TRUE))
224
- all.equal(SpatialLinesLengths(lns), unname(gLength(lns, byid=TRUE)))
216
+ summary(lns)
217
+ summary(SpatialLinesLengths(as(lns, "Spatial")))
225
218
 
226
219
 
227
220
  ###################################################
228
221
  ### code chunk number 39: cm2.Rnw:791-794
229
222
  ###################################################
230
- GI <- gIntersection(lns, olinda_utm, byid=TRUE)
223
+ crs(lns) <- crs(olinda_utm)
224
+ GI <- intersect(lns, olinda_utm)
231
225
  class(GI)
232
- length(row.names(GI))
226
+ nrow(GI)
233
227
 
234
228
 
235
229
  ###################################################
236
230
  ### code chunk number 40: cm2.Rnw:817-827
237
231
  ###################################################
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
232
 
249
233
 
250
234
  ###################################################
251
235
  ### code chunk number 41: cm2.Rnw:846-848
252
236
  ###################################################
253
- tree <- gBinarySTRtreeQuery(lns, olinda_utm)
254
- table(sapply(tree, length))
255
237
 
256
238
 
257
239
  ###################################################
258
240
  ### code chunk number 42: cm2.Rnw:872-880
259
241
  ###################################################
260
- res1 <- numeric(length=length(tree))
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))
265
- }
266
- }
267
- all.equal(olinda_utm$stream_len, res1)
242
+ GIa <- aggregate(GI, by="ID", fun=head, n=1)
243
+ GIa$stream_len <- SpatialLinesLengths(as(GIa, "Spatial"))
244
+ olinda_utm_1 <- merge(olinda_utm, as.data.frame(GIa[, c("ID", "stream_len")]), by="ID", all=TRUE)
245
+ olinda_utm_1$stream_len[is.na(olinda_utm_1$stream_len)] <- 0
268
246
 
269
247
 
270
248
  ###################################################
271
249
  ### code chunk number 43: cm2.Rnw:899-901
272
250
  ###################################################
273
- buf50m <- gBuffer(lns, width=50)
274
- length(slot(buf50m, "polygons"))
251
+ buf50m <- buffer(lns, width=50)
252
+ nrow(buf50m)
275
253
 
276
254
 
277
255
  ###################################################
278
256
  ### code chunk number 44: cm2.Rnw:918-926
279
257
  ###################################################
280
- GI1 <- gIntersection(buf50m, olinda_utm, byid=TRUE)
281
- 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)
286
- olinda_utm$buf_area <- res
287
- olinda_utm$prop_50m <- olinda_utm$buf_area/olinda_utm$area
258
+ GI1 <- intersect(buf50m, olinda_utm_1)
259
+ GI1a <- aggregate(GI1, by="ID", fun=head, n=1)
260
+ GI1a$buf_area <- area(as(GI1a, "Spatial"))
261
+ olinda_utm_2 <- merge(olinda_utm_1, as.data.frame(GI1a[, c("ID", "buf_area")]), by="ID", all=TRUE)
262
+ olinda_utm_2$buf_area[is.na(olinda_utm_2$buf_area)] <- 0
263
+ olinda_utm_2$prop_50m <- olinda_utm_2$buf_area/olinda_utm_2$area
288
264
 
289
265
 
290
266
  ###################################################
291
267
  ### code chunk number 45: cm2.Rnw:942-943
292
268
  ###################################################
293
- stream_inside <- gIntersection(lns, bounds)
269
+ stream_inside <- intersect(lns, bounds)
294
270
 
295
271
 
296
272
  ###################################################
@@ -298,7 +274,7 @@
298
274
  ###################################################
299
275
  library(RColorBrewer)
300
276
  bl <- brewer.pal(5, "Blues")
301
- spplot(olinda_utm, "prop_50m", col.regions=colorRampPalette(bl)(20), col="transparent", sp.layout=list("sp.lines", stream_inside), colorkey=list(space="right", labels=list(cex=0.7), width=1))
277
+ spplot(as(olinda_utm_2, "Spatial"), "prop_50m", col.regions=colorRampPalette(bl)(20), col="transparent", sp.layout=list("sp.lines", as(stream_inside, "Spatial")), colorkey=list(space="right", labels=list(cex=0.7), width=1))
302
278
 
303
279
 
304
280
  ###################################################
@@ -349,8 +325,8 @@
349
325
  ###################################################
350
326
  ### code chunk number 55: cm2.Rnw:1189-1196
351
327
  ###################################################
352
- TM0 <- as(TM, "SpatialPixelsDataFrame")
353
- TM1 <- TM0[bounds, ]
328
+ TM0 <- as(brick(TM), "SpatialPixelsDataFrame")
329
+ TM1 <- TM0[as(bounds, "Spatial"), ]
354
330
  PC <- prcomp(as(TM1, "data.frame")[,1:6], center=TRUE, scale.=TRUE)
355
331
  PCout <- predict(PC)
356
332
  TM1$PC1 <- PCout[,1]
@@ -361,29 +337,28 @@
361
337
  ###################################################
362
338
  ### code chunk number 56: cm2.Rnw:1202-1210
363
339
  ###################################################
364
- spplot(TM1, c("PC1", "PC2"), at=seq(-17, 17, length.out=21), col.regions=rev(colorRampPalette(brewer.pal(10, "PiYG"))(20)), sp.layout=list("sp.lines", as(olinda_utm, "SpatialLines"), lwd=0.5), colorkey=list(space="right", labels=list(cex=0.8), width=1))
340
+ spplot(TM1, c("PC1", "PC2"), at=seq(-17, 17, length.out=21), col.regions=rev(colorRampPalette(brewer.pal(10, "PiYG"))(20)), sp.layout=list("sp.lines", as(as(olinda_utm, "Spatial"), "SpatialLines"), lwd=0.5), colorkey=list(space="right", labels=list(cex=0.8), width=1))
365
341
 
366
342
 
367
343
  ###################################################
368
344
  ### code chunk number 57: cm2.Rnw:1226-1230
369
345
  ###################################################
370
- o_mean <- over(olinda_utm, TM1[,c("PC1", "PC2", "PC3")])
346
+ o_mean <- over(as(olinda_utm, "Spatial"), TM1[,c("PC1", "PC2", "PC3")])
371
347
  str(o_mean)
372
348
  row.names(o_mean) <- row.names(olinda_utm)
373
- library(maptools)
374
- olinda_utmA <- spCbind(olinda_utm, o_mean)
349
+ olinda_utmA <- cbind(as(olinda_utm, "Spatial"), o_mean)
375
350
 
376
351
 
377
352
 
378
353
  ###################################################
379
354
  ### code chunk number 58: cm2.Rnw:1241-1248
380
355
  ###################################################
381
- o_median <- over(olinda_utm, TM1[,c("PC1", "PC2", "PC3")], fn=median)
356
+ o_median <- over(as(olinda_utm, "Spatial"), TM1[,c("PC1", "PC2", "PC3")], fn=median)
382
357
  row.names(o_median) <- row.names(olinda_utmA)
383
358
  names(o_median) <- paste(names(o_median), "med", sep="_")
384
- olinda_utmB <- spCbind(olinda_utmA, o_median)
359
+ olinda_utmB <- cbind(olinda_utmA, o_median)
385
360
  TM1$count <- 1
386
- o_count <- over(olinda_utm, TM1[,"count"], fn=sum)
361
+ o_count <- over(as(olinda_utm, "Spatial"), TM1[,"count"], fn=sum)
387
362
  olinda_utmB$count <- o_count$count
388
363
 
389
364
 
@@ -396,10 +371,10 @@
396
371
  ###################################################
397
372
  ### code chunk number 60: cm2.Rnw:1272-1277
398
373
  ###################################################
399
- o_dem_median <- over(olinda_utm, dem, fn=median)
400
- olinda_utmB$dem_median <- o_dem_median$band1
374
+ o_dem_median <- over(as(olinda_utm, "Spatial"), as(raster(dem), "SpatialGridDataFrame"), fn=median)
375
+ olinda_utmB$dem_median <- o_dem_median$olinda_dem_utm25s
401
376
  summary(olinda_utmB$dem_median)
402
- o_ndvi_median <- over(olinda_utm, TM1["ndvi"], fn=median)
377
+ o_ndvi_median <- over(as(olinda_utm, "Spatial"), TM1["ndvi"], fn=median)
403
378
  olinda_utmB$ndvi_median <- o_ndvi_median$ndvi
404
379
 
405
380
 
@@ -413,13 +388,13 @@
413
388
  ### code chunk number 64: cm2.Rnw:1314-1316 (eval = FALSE)
414
389
  ###################################################
415
390
  TMrs <- stack(TM1)
416
- e1 <- extract(TMrs, as(olinda_utm, "SpatialPolygons"), small=FALSE)
391
+ e1 <- extract(TMrs, as(as(olinda_utm, "Spatial"), "SpatialPolygons"), small=FALSE)
417
392
 
418
393
 
419
394
  ###################################################
420
395
  ### code chunk number 65: cm2.Rnw:1328-1329 (eval = FALSE)
421
396
  ###################################################
422
- e2 <- extract(raster(dem), as(olinda_utm, "SpatialPolygons"), small=FALSE)
397
+ e2 <- extract(raster(dem), as(as(olinda_utm, "Spatial"), "SpatialPolygons"), small=FALSE)
423
398
 
424
399
 
425
400
  ###################################################
@@ -443,9 +418,10 @@
443
418
  #o <- over(dem, bounds)
444
419
  #dem$band1 <- dem$band1*o
445
420
  set.seed(9876)
446
- p_r <- spsample(bounds, 1000, type="random")
421
+ p_r <- spsample(as(bounds, "Spatial"), 1000, type="random")
447
422
  length(p_r)
448
- dem <- dem[bounds,]
423
+ dem <- as(raster(dem), "SpatialGridDataFrame")
424
+ dem <- dem[as(bounds, "Spatial"),]
449
425
  dem_sp <- as(dem, "SpatialPixelsDataFrame")
450
426
  g_r <- spsample(dem_sp, 1000, type="random")
451
427
  length(g_r)
@@ -471,12 +447,12 @@
471
447
  title(main="grid_random")
472
448
  plot(g_rg, cex=0.15)
473
449
  title(main="grid_regular")
474
- plot(ecdf(p_r_dem$band1), verticals=TRUE, do.p=FALSE, ann=FALSE,
450
+ plot(ecdf(p_r_dem$olinda_dem_utm25s), verticals=TRUE, do.p=FALSE, ann=FALSE,
475
451
  col.hor="green", col.vert="green")
476
452
  title(main="ECDF")
477
- plot(ecdf(g_r_dem$band1), verticals=TRUE, do.p=FALSE,
453
+ plot(ecdf(g_r_dem$olinda_dem_utm25s), verticals=TRUE, do.p=FALSE,
478
454
  col.hor="blue", col.vert="blue", add=TRUE)
479
- plot(ecdf(g_rg_dem$band1), verticals=TRUE, do.p=FALSE,
455
+ plot(ecdf(g_rg_dem$olinda_dem_utm25s), verticals=TRUE, do.p=FALSE,
480
456
  col.hor="red", col.vert="red", add=TRUE)
481
457
  abline(h=c(0.25,0.5,0.75), lty=2, col="grey")
482
458
  legend("bottomright", c("polygon random", "grid random", "grid regular"),
@@ -486,11 +462,11 @@
486
462
  ###################################################
487
463
  ### code chunk number 72: cm2.Rnw:1468-1477
488
464
  ###################################################
489
- tab <- rbind(polygon_random=c(fivenum(p_r_dem$band1),
465
+ tab <- rbind(polygon_random=c(fivenum(p_r_dem$olinda_dem_utm25s),
490
466
  nrow(p_r_dem)),
491
- grid_random=c(fivenum(g_r_dem$band1),
467
+ grid_random=c(fivenum(g_r_dem$olinda_dem_utm25s),
492
468
  nrow(g_r_dem)),
493
- grid_regular=c(fivenum(g_rg_dem$band1),
469
+ grid_regular=c(fivenum(g_rg_dem$olinda_dem_utm25s),
494
470
  nrow(g_rg_dem)))
495
471
  colnames(tab) <- c("minimum", "lower-hinge", "median", "upper-hinge",
496
472
  "maximum", "n")
@@ -500,12 +476,14 @@
500
476
  ###################################################
501
477
  ### code chunk number 73: cm2.Rnw:1499-1505
502
478
  ###################################################
503
- o_sp <- as(olinda_utm, "SpatialPolygons")
479
+ o_sp <- as(as(olinda_utm, "Spatial"), "SpatialPolygons")
504
480
  whichPoly <- over(p_r, o_sp)
505
- whichPoly1 <- gContains(o_sp, p_r, byid=TRUE)
506
- whichPoly1a <- apply(unname(whichPoly1), 1, which)
481
+ p_r_v <- vect(p_r)
482
+ crs(p_r_v) <- crs(olinda_utm)
483
+ whichPoly1 <- relate(olinda_utm, p_r_v, relation="contains", pairs=FALSE)
484
+ whichPoly1a <- apply(unname(whichPoly1), 2, which)
507
485
  table(sapply(whichPoly1a, length))
508
- all.equal(whichPoly, whichPoly1a)
486
+ all.equal( unname(whichPoly),whichPoly1a)
509
487
 
510
488
 
511
489
  ###################################################
@@ -515,7 +493,13 @@
515
493
  p4s <- CRS("+proj=longlat +datum=WGS84")
516
494
  Hels <- SpatialPoints(hels, proj4string=p4s)
517
495
  d041224 <- as.POSIXct("2004-12-24", tz="EET")
518
- sunriset(Hels, d041224, direction="sunrise", POSIXct.out=TRUE)
496
+ solartime::computeSunriseHour(d041224, latDeg=hels[1,2], longDeg=hels[1,1])
497
+
498
+
499
+ (sI <- sessionInfo()) # check: no sp?
519
500
 
501
+ "rgdal" %in% c(names(sI$otherPkgs), names(sI$loadedOnly))
502
+ "rgeos" %in% c(names(sI$otherPkgs), names(sI$loadedOnly))
503
+ "maptools" %in% c(names(sI$otherPkgs), names(sI$loadedOnly))
520
504
 
521
505