@@ -105,7 +105,8 @@
|
|
105
105
|
### code chunk number 25: cm.Rnw:528-533
|
106
106
|
###################################################
|
107
107
|
m <- matrix(c(0,0,1,1), ncol=2, dimnames=list(NULL, c("min", "max")))
|
108
|
-
|
108
|
+
library(sf)
|
109
|
+
crs <- as(st_crs(as.character(NA)), "CRS")
|
109
110
|
crs
|
110
111
|
S <- Spatial(bbox=m, proj4string=crs)
|
111
112
|
S
|
@@ -115,7 +116,8 @@
|
|
115
116
|
### code chunk number 27: cm.Rnw:555-558
|
116
117
|
###################################################
|
117
118
|
bb <- matrix(c(350, 85, 370, 95), ncol=2, dimnames=list(NULL, c("min", "max")))
|
118
|
-
res <- try(Spatial(bb, proj4string=
|
119
|
+
res <- try(Spatial(bb, proj4string=as(sf::st_crs("+proj=longlat +ellps=WGS84"),
|
120
|
+
"CRS")))
|
119
121
|
|
120
122
|
|
121
123
|
###################################################
|
@@ -136,7 +138,7 @@
|
|
136
138
|
###################################################
|
137
139
|
### code chunk number 33: cm.Rnw:671-674
|
138
140
|
###################################################
|
139
|
-
llCRS <-
|
141
|
+
llCRS <- as(st_crs("+proj=longlat +ellps=WGS84"), "CRS")
|
140
142
|
CRAN_sp <- SpatialPoints(CRAN_mat, proj4string=llCRS)
|
141
143
|
summary(CRAN_sp)
|
142
144
|
|
@@ -151,7 +153,7 @@
|
|
151
153
|
### code chunk number 35: cm.Rnw:720-724
|
152
154
|
###################################################
|
153
155
|
proj4string(CRAN_sp)
|
154
|
-
proj4string(CRAN_sp) <-
|
156
|
+
proj4string(CRAN_sp) <- as(st_crs(as.character(NA)), "CRS")
|
155
157
|
proj4string(CRAN_sp)
|
156
158
|
proj4string(CRAN_sp) <- llCRS
|
157
159
|
|
@@ -269,21 +271,21 @@
|
|
269
271
|
turtle_df1$lon <- ifelse(turtle_df1$lon < 0, turtle_df1$lon+360, turtle_df1$lon)
|
270
272
|
turtle_sp <- turtle_df1[order(turtle_df1$timestamp),]
|
271
273
|
coordinates(turtle_sp) <- c("lon", "lat")
|
272
|
-
proj4string(turtle_sp) <-
|
274
|
+
proj4string(turtle_sp) <- as(st_crs("+proj=longlat +ellps=WGS84"), "CRS")
|
273
275
|
|
274
276
|
|
275
277
|
###################################################
|
276
278
|
### code chunk number 55: cm.Rnw:1031-1034
|
277
279
|
###################################################
|
278
|
-
library(
|
279
|
-
|
280
|
-
pac <-
|
280
|
+
library(maps)
|
281
|
+
pac <- map("world", fill=TRUE, plot=FALSE, wrap=c(0,360))
|
282
|
+
pac <- as(st_as_sf(pac), "Spatial")
|
281
283
|
|
282
284
|
|
283
285
|
###################################################
|
284
286
|
### code chunk number 56: cm.Rnw:1039-1051
|
285
287
|
###################################################
|
286
|
-
plot(pac
|
288
|
+
plot(pac, axes=TRUE, xlim=c(130,250), ylim=c(15,60), col="khaki2", xaxs="i", yaxs="i")
|
287
289
|
plot(turtle_sp, add=TRUE)
|
288
290
|
m_rle <- rle(months(turtle_sp$timestamp))
|
289
291
|
clen <- cumsum(m_rle$lengths[-length(m_rle$lengths)])-1
|
@@ -308,10 +310,9 @@
|
|
308
310
|
### code chunk number 59: cm.Rnw:1140-1146
|
309
311
|
###################################################
|
310
312
|
library(maps)
|
311
|
-
japan <- map("world", "japan", plot=FALSE)
|
312
|
-
p4s <-
|
313
|
-
|
314
|
-
SLjapan <- map2SpatialLines(japan, proj4string=p4s)
|
313
|
+
japan <- map("world", "japan", fill=FALSE, plot=FALSE)
|
314
|
+
p4s <- st_crs("+proj=longlat +ellps=WGS84")
|
315
|
+
SLjapan <- as(st_geometry(st_as_sf(japan, fill=FALSE, crs=p4s)), "Spatial")
|
315
316
|
str(SLjapan, max.level=2)
|
316
317
|
|
317
318
|
|
@@ -325,16 +326,22 @@
|
|
325
326
|
###################################################
|
326
327
|
### code chunk number 61: cm.Rnw:1189-1191
|
327
328
|
###################################################
|
328
|
-
|
329
|
+
data(volcano)
|
330
|
+
x = contourLines(volcano)
|
331
|
+
xx = lapply(x, function(x) cbind(x$x, x$y))
|
332
|
+
xxx = split(xx, sapply(x, "[[", "level"))
|
333
|
+
x2a = st_sfc(lapply(xxx, st_multilinestring))
|
334
|
+
x2 = st_as_sf(x2a, level=names(xxx))
|
335
|
+
volcano_sl <- as(x2, "Spatial")
|
336
|
+
#volcano_sl <- ContourLines2SLDF(contourLines(volcano))
|
329
337
|
t(slot(volcano_sl, "data"))
|
330
338
|
|
331
339
|
|
332
340
|
###################################################
|
333
341
|
### code chunk number 63: cm.Rnw:1217-1220
|
334
342
|
###################################################
|
335
|
-
|
336
|
-
auck_shore <-
|
337
|
-
summary(auck_shore)
|
343
|
+
# MapGen format defunct https://www.ngdc.noaa.gov/mgg/shorelines/shorelines.html
|
344
|
+
auck_shore <- as(st_geometry(st_read("auck_shore.gpkg")), "Spatial")
|
338
345
|
|
339
346
|
|
340
347
|
###################################################
|
@@ -351,7 +358,7 @@
|
|
351
358
|
Polygons(list(Polygon(slot(slot(x, "Lines")[[1]], "coords"))),
|
352
359
|
ID=slot(x, "ID"))
|
353
360
|
}),
|
354
|
-
proj4string=
|
361
|
+
proj4string=as(st_crs("+proj=longlat +ellps=WGS84"), "CRS"))
|
355
362
|
|
356
363
|
|
357
364
|
###################################################
|
@@ -403,7 +410,7 @@
|
|
403
410
|
Polygons(list(Polygon(slot(slot(x, "Lines")[[1]], "coords"))),
|
404
411
|
ID=slot(x, "ID"))
|
405
412
|
}),
|
406
|
-
proj4string=
|
413
|
+
proj4string=as(st_crs("+proj=longlat +ellps=WGS84"), "CRS"))
|
407
414
|
summary(islands_sp)
|
408
415
|
slot(islands_sp, "plotOrder")
|
409
416
|
order(sapply(slot(islands_sp, "polygons"),
|
@@ -415,10 +422,10 @@
|
|
415
422
|
###################################################
|
416
423
|
library(maps)
|
417
424
|
state.map <- map("state", plot=FALSE, fill=TRUE)
|
418
|
-
|
419
|
-
|
420
|
-
state.sp <-
|
421
|
-
|
425
|
+
state_sf <- st_as_sf(state.map, fill=TRUE, group=TRUE,
|
426
|
+
crs=st_crs("+proj=longlat +ellps=WGS84"))
|
427
|
+
state.sp <- as(st_geometry(state_sf), "Spatial")
|
428
|
+
row.names(state.sp) <- state_sf$ID
|
422
429
|
|
423
430
|
|
424
431
|
###################################################
|
@@ -477,23 +484,27 @@
|
|
477
484
|
###################################################
|
478
485
|
### code chunk number 82: cm.Rnw:1643-1644
|
479
486
|
###################################################
|
480
|
-
library(rgeos)
|
487
|
+
#library(rgeos)
|
481
488
|
|
482
489
|
|
483
490
|
###################################################
|
484
491
|
### code chunk number 83: cm.Rnw:1646-1648
|
485
492
|
###################################################
|
486
|
-
manitoulin_sp <-
|
493
|
+
comment(slot(manitoulin_sp, "polygons")[[1]]) <- rep("0", 19)
|
494
|
+
comment(manitoulin_sp) <- "TRUE"
|
495
|
+
manitoulin_sf <- st_make_valid(st_geometry(st_as_sf(manitoulin_sp)))
|
496
|
+
manitoulin_sp <- as(manitoulin_sf, "Spatial")
|
487
497
|
sapply(slot(manitoulin_sp, "polygons"), comment)
|
498
|
+
# eastern edge lake assigned here to exterior, minor differences
|
488
499
|
|
489
500
|
|
490
501
|
###################################################
|
491
502
|
### code chunk number 84: cm.Rnw:1674-1687
|
492
503
|
###################################################
|
493
|
-
plot(manitoulin_sp,
|
504
|
+
plot(manitoulin_sp, bg="lightsteelblue2", col="khaki2", usePolypath=TRUE, xaxs="i", yaxs="i")
|
494
505
|
text(t(sapply(slot(slot(manitoulin_sp, "polygons")[[1]], "Polygons"), function(x) slot(x, "labpt")))[-c(1,2),], label=high$polydata$level[-c(1,2)], col="black", font=2)
|
495
506
|
cmt <- unlist(strsplit(sapply(slot(manitoulin_sp, "polygons"), comment), " "))
|
496
|
-
plot(manitoulin_sp,
|
507
|
+
plot(manitoulin_sp, bg="lightsteelblue2", col="khaki2", usePolypath=TRUE, xaxs="i", yaxs="i")
|
497
508
|
text(t(sapply(slot(slot(manitoulin_sp, "polygons")[[1]], "Polygons"), function(x) slot(x, "labpt")))[-c(1,2),], label=cmt[-c(1,2)], col="black", font=2)
|
498
509
|
|
499
510
|
|
@@ -524,7 +535,7 @@
|
|
524
535
|
###################################################
|
525
536
|
### code chunk number 89: cm.Rnw:1772-1775
|
526
537
|
###################################################
|
527
|
-
p4s <-
|
538
|
+
p4s <- slot(manitoulin_sp, "proj4string")
|
528
539
|
manitoulin_SG <- SpatialGrid(manitoulin_grd, proj4string=p4s)
|
529
540
|
summary(manitoulin_SG)
|
530
541
|
|
@@ -609,57 +620,56 @@
|
|
609
620
|
###################################################
|
610
621
|
### code chunk number 100: cm.Rnw:1983-1984
|
611
622
|
###################################################
|
612
|
-
library(
|
623
|
+
library(stars)
|
613
624
|
|
614
625
|
|
615
626
|
###################################################
|
616
627
|
### code chunk number 101: cm.Rnw:1986-1987
|
617
628
|
###################################################
|
618
|
-
r <-
|
629
|
+
r <- read_stars("70042108.tif", proxy = TRUE)
|
619
630
|
|
620
631
|
|
621
632
|
###################################################
|
622
633
|
### code chunk number 103: cm.Rnw:1992-1998
|
623
634
|
###################################################
|
624
635
|
class(r)
|
625
|
-
inMemory(r)
|
626
636
|
object.size(r)
|
627
|
-
|
628
|
-
|
629
|
-
|
637
|
+
max(as.vector(st_as_stars(r)[[1]]), na.rm = TRUE)
|
638
|
+
min(as.vector(st_as_stars(r)[[1]]), na.rm = TRUE)
|
639
|
+
inherits(r, "stars_proxy") # out-of-memory?
|
640
|
+
object.size(r)
|
641
|
+
object.size(st_as_stars(r))
|
630
642
|
|
631
643
|
|
632
644
|
###################################################
|
633
645
|
### code chunk number 104: cm.Rnw:2012-2024
|
634
646
|
###################################################
|
635
|
-
out
|
636
|
-
|
637
|
-
|
638
|
-
|
639
|
-
|
640
|
-
v[v <= 0] <- NA
|
641
|
-
writeValues(out, v, bs$row[i])
|
642
|
-
}
|
643
|
-
out <- writeStop(out)
|
644
|
-
cellStats(out, min)
|
645
|
-
cellStats(out, max)
|
646
|
-
inMemory(out)
|
647
|
+
out = r
|
648
|
+
out[ r <= 0 ] = NA
|
649
|
+
dout <- st_as_stars(out)
|
650
|
+
max(as.vector(dout[[1]]), na.rm = TRUE)
|
651
|
+
min(as.vector(dout[[1]]), na.rm = TRUE)
|
647
652
|
|
648
653
|
|
649
654
|
###################################################
|
650
655
|
### code chunk number 106: cm.Rnw:2044-2052
|
651
656
|
###################################################
|
652
|
-
|
653
|
-
|
657
|
+
r1 <- as(dout, "Spatial")
|
658
|
+
image(r1, col = terrain.colors(100))
|
659
|
+
plot(auck_gshhs, add = TRUE)
|
654
660
|
|
655
661
|
|
656
662
|
###################################################
|
657
663
|
### code chunk number 107: cm.Rnw:2069-2073
|
658
664
|
###################################################
|
659
|
-
r1 <- as(out, "SpatialGridDataFrame")
|
660
665
|
summary(r1)
|
661
|
-
r2 <- as(
|
666
|
+
r2 <- as(dout, "SpatRaster")
|
662
667
|
summary(r2)
|
663
668
|
|
669
|
+
(sI <- sessionInfo()) # check: no sp?
|
670
|
+
|
671
|
+
"rgdal" %in% c(names(sI$otherPkgs), names(sI$loadedOnly))
|
672
|
+
"rgeos" %in% c(names(sI$otherPkgs), names(sI$loadedOnly))
|
673
|
+
"maptools" %in% c(names(sI$otherPkgs), names(sI$loadedOnly))
|
664
674
|
|
665
675
|
|