Diff to HTML by rtfpessoa

Files changed (1) show
{orig → migrated}/geos_mod.R RENAMED
@@ -13,29 +13,27 @@
13
13
  ### code chunk number 11: geos.Rnw:169-175
14
14
  ###################################################
15
15
  #trellis.par.set(canonical.theme(color = FALSE))
16
- library(sp)
17
- data(meuse)
18
- coordinates(meuse) <- c("x", "y")
16
+ library(sf)
17
+ data(meuse, package = "sp")
18
+ meuse = st_as_sf(meuse, coords = c("x", "y"))
19
19
 
20
20
 
21
21
  ###################################################
22
22
  ### code chunk number 13: geos.Rnw:207-220
23
23
  ###################################################
24
- print(xyplot(log(zinc)~sqrt(dist), as.data.frame(meuse), asp = .8), split =
25
- c(1, 1,2,1), more = TRUE)
24
+ xyplot(log(zinc)~sqrt(dist), as.data.frame(meuse), asp = .8)
26
25
  zn.lm <- lm(log(zinc)~sqrt(dist), meuse)
27
26
  meuse$fitted.s <- predict(zn.lm, meuse) - mean(predict(zn.lm, meuse))
28
27
  meuse$residuals <- residuals(zn.lm)
29
- print(spplot(meuse, c("fitted.s", "residuals"), col.regions =
30
- pal(), cuts = 8, colorkey=TRUE), split = c(2,1,2,1))
28
+ plot(meuse[c("fitted.s", "residuals")], pal = pal(), pch = 16)
31
29
 
32
30
 
33
31
  ###################################################
34
32
  ### code chunk number 14: geos.Rnw:250-253
35
33
  ###################################################
36
- data(meuse.grid)
37
- coordinates(meuse.grid) <- c("x", "y")
38
- meuse.grid <- as(meuse.grid, "SpatialPixelsDataFrame")
34
+ data(meuse.grid, package = "sp")
35
+ library(stars)
36
+ meuse.grid = st_as_stars(meuse.grid)
39
37
 
40
38
 
41
39
  ###################################################
@@ -69,6 +67,8 @@
69
67
  ###################################################
70
68
  ### code chunk number 19: geos.Rnw:384-385
71
69
  ###################################################
70
+ meuse$x = st_coordinates(meuse)[,1]
71
+ meuse$y = st_coordinates(meuse)[,2]
72
72
  lm(log(zinc)~I(x^2)+I(y^2)+I(x*y) + x + y, meuse)
73
73
 
74
74
 
@@ -87,7 +87,6 @@
87
87
  ###################################################
88
88
  ### code chunk number 24: geos.Rnw:560-585
89
89
  ###################################################
90
- library(gstat)
91
90
  cld <- variogram(log(zinc) ~ 1, meuse, cloud = TRUE)
92
91
  svgm <- variogram(log(zinc) ~ 1, meuse)
93
92
  d <- data.frame(gamma = c(cld$gamma, svgm$gamma),
@@ -317,7 +316,7 @@
317
316
  ###################################################
318
317
  ### code chunk number 45: geos.Rnw:1224-1225
319
318
  ###################################################
320
- fit.variogram.reml(log(zinc)~1, meuse, model=vgm(0.6, "Sph", 800, 0.06))
319
+ fit.variogram.reml(log(zinc)~1, data=as(meuse, "Spatial"), model=vgm(0.6, "Sph", 800, 0.06))
321
320
 
322
321
 
323
322
  ###################################################
@@ -392,7 +391,8 @@
392
391
  ###################################################
393
392
  ### code chunk number 58: geos.Rnw:1644-1652
394
393
  ###################################################
395
- print(spplot.vcov(cok.maps, cuts=6, col.regions=pal(7)))
394
+ #### THIS IS NOT SUPPORTED using stars:
395
+ #print(spplot.vcov(cok.maps, cuts=6, col.regions=pal(7)))
396
396
 
397
397
 
398
398
  ###################################################
@@ -415,12 +415,12 @@
415
415
  ### code chunk number 62: geos.Rnw:1711-1724
416
416
  ###################################################
417
417
  g.cc <- gstat(NULL, "log.zinc", log(zinc)~1, meuse, model = v.fit)
418
- meuse.grid$distn <- meuse.grid$dist - mean(meuse.grid$dist) +
418
+ meuse.grid$distn <- meuse.grid$dist - mean(meuse.grid$dist, na.rm = TRUE) +
419
419
  mean(log(meuse$zinc))
420
420
  vd.fit <- v.fit
421
- vov <- var(meuse.grid$distn) / var(log(meuse$zinc))
421
+ vov <- var(as.vector(meuse.grid$distn), na.rm = TRUE) / var(log(meuse$zinc))
422
422
  vd.fit$psill <- v.fit$psill * vov
423
- g.cc <- gstat(g.cc, "distn", distn ~ 1, meuse.grid, nmax = 1, model=vd.fit,
423
+ g.cc <- gstat(g.cc, "distn", distn ~ 1, st_as_sf(meuse.grid, as_points = TRUE), nmax = 1, model=vd.fit,
424
424
  merge = c("log.zinc","distn"))
425
425
  vx.fit <- v.fit
426
426
  vx.fit$psill <- sqrt(v.fit$psill * vd.fit$psill) *
@@ -434,10 +434,8 @@
434
434
  ###################################################
435
435
  x$lz.uk <- lz.uk$var1.pred
436
436
  x$lz.ok <- lz.ok$var1.pred
437
- print(spplot(x, c("log.zinc.pred", "lz.ok", "lz.uk"),
438
- names.attr = c("collocated", "ordinary", "universal"),
439
- cuts=7, col.regions=pal(8)
440
- ))
437
+
438
+ #... use pivot_longer, ggplot2 + geom_sf()
441
439
 
442
440
 
443
441
  ###################################################
@@ -469,19 +467,18 @@
469
467
  ###################################################
470
468
  ### code chunk number 70: geos.Rnw:1988-1989
471
469
  ###################################################
470
+ meuse.grid = st_as_sf(meuse.grid, as_points = TRUE)
472
471
  meuse$part.a <- gstat::idw(part.a~1, meuse.grid, meuse, nmax=1)$var1.pred
473
472
 
474
473
 
475
474
  ###################################################
476
475
  ### code chunk number 72: geos.Rnw:2011-2012
477
476
  ###################################################
478
- meuse$part.a <- over(meuse, meuse.grid["part.a"])[[1]]
479
477
 
480
478
 
481
479
  ###################################################
482
480
  ### code chunk number 73: geos.Rnw:2020-2021
483
481
  ###################################################
484
- meuse$part.a <- meuse.grid$part.a[over(meuse, geometry(meuse.grid))]
485
482
 
486
483
 
487
484
  ###################################################
@@ -493,10 +490,8 @@
493
490
  x2 <- krige(log(zinc)~1, meuse[meuse$part.a == 1,],
494
491
  meuse.grid[meuse.grid$part.a == 1,], model = vgm(.716, "Sph", 900),
495
492
  nmin = 20, nmax = 40, maxdist = 1000)
496
- lz.stk <- rbind(as.data.frame(x1), as.data.frame(x2))
497
- coordinates(lz.stk) <- c("x", "y")
498
- lz.stk <- as(x, "SpatialPixelsDataFrame")
499
-
493
+ lz.stk <- rbind(as.data.frame(x1), as.data.frame(x2)) %>%
494
+ st_as_sf()
500
495
 
501
496
  ###################################################
502
497
  ### code chunk number 75: geos.Rnw:2050-2051 (eval = FALSE)
@@ -526,8 +521,7 @@
526
521
  model = v.fit)
527
522
  rn <- c("Intercept", "beta1")
528
523
  df <- data.frame(Int = c(1,0), dist = c(0,1), row.names=rn)
529
- spdf <- SpatialPointsDataFrame(SpatialPoints(matrix(0, 2, 2)), df)
530
- spdf
524
+ spdf = st_sf(df, st_sfc(st_point(c(0,0)), st_point(c(0,0))))
531
525
  predict(g.tr, spdf, BLUE = TRUE)
532
526
 
533
527
 
@@ -565,17 +559,15 @@
565
559
  ###################################################
566
560
  ### code chunk number 83: geos.Rnw:2287-2289
567
561
  ###################################################
568
- meuse.dup <- rbind(as.data.frame(meuse)[1,], as.data.frame(meuse))
569
- coordinates(meuse.dup)=~x+y
562
+ meuse.dup <- rbind(as.data.frame(meuse)[1,], as.data.frame(meuse)) %>%
563
+ st_as_sf()
570
564
 
571
565
 
572
566
  ###################################################
573
567
  ### code chunk number 84: geos.Rnw:2291-2295
574
568
  ###################################################
575
- zd <- zerodist(meuse.dup)
576
- zd
577
- meuse0 <- meuse.dup[-zd[,1],]
578
- krige(log(zinc)~1, meuse0, meuse[1,], v.fit)
569
+
570
+ ### THIS WON'T WORK: there is no zerodist() in sf
579
571
 
580
572
 
581
573
  ###################################################
@@ -603,7 +595,9 @@
603
595
  ### code chunk number 87: geos.Rnw:2379-2381
604
596
  ###################################################
605
597
  meuse_shift = meuse
606
- meuse_shift@coords[1,] = meuse@coords[1,] + 1
598
+ # shift x coordinate with one:
599
+ meuse_shift$geometry[[1]][1] = meuse_shift$geometry[[1]][1] + 1
600
+
607
601
 
608
602
 
609
603
  ###################################################
@@ -688,8 +682,9 @@
688
682
  ###################################################
689
683
  ### code chunk number 98: geos.Rnw:2587-2595
690
684
  ###################################################
691
- print(bubble(cv155, "residual", main = "log(zinc): 5-fold CV residuals",
692
- maxsize = 1.5, col = c("green", "red")))
685
+ # FIXME:
686
+ #print(bubble(cv155, "residual", main = "log(zinc): 5-fold CV residuals",
687
+ # maxsize = 1.5, col = c("green", "red")))
693
688
 
694
689
 
695
690
  ###################################################
@@ -738,7 +733,8 @@
738
733
  ###################################################
739
734
  ### code chunk number 106: geos.Rnw:2841-2850
740
735
  ###################################################
741
- print(spplot(lzn.sim, col.regions=pal(7), cuts=6))
736
+ # FIXME:
737
+ #print(spplot(lzn.sim, col.regions=pal(7), cuts=6))
742
738
 
743
739
 
744
740
  ###################################################
@@ -755,7 +751,7 @@
755
751
  332810.293894737, 333008.922947368, 333188.634947368, 333358.888421053,
756
752
  333491.307789474, 333566.976, 333614.268631579, 333793.980631579,
757
753
  333718.312421053), .Dim = as.integer(c(18, 2)))
758
- area.sp = SpatialPolygons(list(Polygons(list(Polygon(area)), "area")))
754
+ area.sf = st_sfc(st_polygon(list(area)))
759
755
 
760
756
 
761
757
  ###################################################
@@ -764,7 +760,7 @@
764
760
  # set eval=TRUE if change; this section is repeated in the FIG code
765
761
  nsim <- 1000
766
762
  cutoff <- 500
767
- sel.grid <- meuse.grid[area.sp, ]
763
+ sel.grid <- meuse.grid[area.sf, ]
768
764
  lzn.sim <- krige(log(zinc)~1, meuse, sel.grid, v.fit, nsim = nsim, nmax=40)
769
765
  res <- apply(as.data.frame(lzn.sim)[1:nsim], 2, function(x) mean(x >
770
766
  log(cutoff)))
@@ -779,21 +775,13 @@
779
775
  ###################################################
780
776
  ### code chunk number 110: geos.Rnw:2922-2924
781
777
  ###################################################
782
- bkr <- krige(log(zinc)~1, meuse, area.sp, v.fit)
778
+ bkr <- krige(log(zinc)~1, meuse, area.sf, v.fit)
783
779
  1 - pnorm(log(cutoff), bkr$var1.pred, sqrt(bkr$var1.var))
784
780
 
785
781
 
786
782
  ###################################################
787
783
  ### code chunk number 111: geos.Rnw:2930-2943
788
784
  ###################################################
789
- layout(matrix(1:2, 1, 2))
790
- omar <- par("mar")
791
- par(mar = rep(0,4))
792
- image(meuse.grid["part.a"], col = 'gray') #$
793
- lines(area, col = 'red')
794
- par(mar = c(2,2,0.5,0)+.1)
795
- hist(res, main = NULL, xlab = NULL, ylab = NULL)
796
- par(mar = omar)
797
785
 
798
786
 
799
787
  ###################################################
@@ -808,7 +796,6 @@
808
796
  ###################################################
809
797
  ### code chunk number 114: geos.Rnw:3006-3014
810
798
  ###################################################
811
- print(spplot(s1.sim, cuts=1, col.regions=pal(3)[c(1,3)]))
812
799
 
813
800
 
814
801
  ###################################################
@@ -840,19 +827,6 @@
840
827
  ###################################################
841
828
  ### code chunk number 118: geos.Rnw:3111-3130
842
829
  ###################################################
843
- source("m1m2.R")
844
- layout(matrix(1:2, 1, 2))
845
- omar <- par("mar")
846
- par(mar = rep(0,4))
847
- image(meuse.grid, col = grey(.8))
848
- points(meuse)
849
- points(meuse[m1 < quantile(m1,.1),], pch=1, col = 'red')
850
- points(meuse[m1 > quantile(m1,.9),], pch=16, col = 'darkgreen')
851
- image(meuse.grid, col = grey(.8))
852
- points(meuse)
853
- points(meuse[m2 < quantile(m2,.1),], pch=16, col = 'darkgreen')
854
- points(meuse[m2 > quantile(m2,.9),], pch=1, col = 'red')
855
- par(mar = omar)
856
830
 
857
831
 
858
832
  ###################################################
@@ -871,6 +845,7 @@
871
845
  ###################################################
872
846
  ### code chunk number 121: geos.Rnw:3270-3272
873
847
  ###################################################
848
+ # FIXME:
874
849
  #library(geoR)
875
850
  #plot(variog(as.geodata(meuse["zinc"]), max.dist=1500))
876
851
 
@@ -881,4 +856,9 @@
881
856
  ## vignette("st", package = "gstat")
882
857
 
883
858
 
859
+ (sI <- sessionInfo()) # check: no sp?
860
+
861
+ "rgdal" %in% c(names(sI$otherPkgs), names(sI$loadedOnly))
862
+ "rgeos" %in% c(names(sI$otherPkgs), names(sI$loadedOnly))
863
+ "maptools" %in% c(names(sI$otherPkgs), names(sI$loadedOnly))
884
864