R Under development (unstable) (2024-09-04 r87094 ucrt) -- "Unsuffered Consequences" Copyright (C) 2024 The R Foundation for Statistical Computing Platform: x86_64-w64-mingw32/x64 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. > Sys.setenv(TZ = "UTC") > > # 0. using sp: > > suppressPackageStartupMessages(library(sp)) > demo(meuse, ask = FALSE) demo(meuse) ---- ~~~~~ > require(sp) > crs = CRS("EPSG:28992") > data("meuse") > coordinates(meuse) <- ~x+y > proj4string(meuse) <- crs > data("meuse.grid") > coordinates(meuse.grid) <- ~x+y > gridded(meuse.grid) <- TRUE > proj4string(meuse.grid) <- crs > data("meuse.riv") > meuse.riv <- SpatialPolygons(list(Polygons(list(Polygon(meuse.riv)),"meuse.riv"))) > proj4string(meuse.riv) <- crs > data("meuse.area") > meuse.area = SpatialPolygons(list(Polygons(list(Polygon(meuse.area)), "area"))) > proj4string(meuse.area) <- crs > suppressPackageStartupMessages(library(gstat)) > v = variogram(log(zinc)~1, meuse) > (v.fit = fit.variogram(v, vgm(1, "Sph", 900, 1))) model psill range 1 Nug 0.05066243 0.0000 2 Sph 0.59060780 897.0209 > k_sp = krige(log(zinc)~1, meuse[-(1:5),], meuse[1:5,], v.fit) [using ordinary kriging] > k_sp_grd = krige(log(zinc)~1, meuse, meuse.grid, v.fit) [using ordinary kriging] > > # 1. using sf: > suppressPackageStartupMessages(library(sf)) > demo(meuse_sf, ask = FALSE, echo = FALSE) > # reloads meuse as data.frame, so > demo(meuse, ask = FALSE) demo(meuse) ---- ~~~~~ > require(sp) > crs = CRS("EPSG:28992") > data("meuse") > coordinates(meuse) <- ~x+y > proj4string(meuse) <- crs > data("meuse.grid") > coordinates(meuse.grid) <- ~x+y > gridded(meuse.grid) <- TRUE > proj4string(meuse.grid) <- crs > data("meuse.riv") > meuse.riv <- SpatialPolygons(list(Polygons(list(Polygon(meuse.riv)),"meuse.riv"))) > proj4string(meuse.riv) <- crs > data("meuse.area") > meuse.area = SpatialPolygons(list(Polygons(list(Polygon(meuse.area)), "area"))) > proj4string(meuse.area) <- crs > > v = variogram(log(zinc)~1, meuse_sf) > (v.fit = fit.variogram(v, vgm(1, "Sph", 900, 1))) model psill range 1 Nug 0.05066243 0.0000 2 Sph 0.59060780 897.0209 > k_sf = krige(log(zinc)~1, meuse_sf[-(1:5),], meuse_sf[1:5,], v.fit) [using ordinary kriging] > > all.equal(k_sp, as(k_sf, "Spatial"), check.attributes = FALSE) [1] TRUE > all.equal(k_sp, as(k_sf, "Spatial"), check.attributes = TRUE) [1] "Attributes: < Component \"bbox\": Attributes: < Component \"dimnames\": Component 1: 2 string mismatches > >" [2] "Attributes: < Component \"coords\": Attributes: < Component \"dimnames\": Component 2: 2 string mismatches > >" [3] "Attributes: < Component \"coords.nrs\": Numeric: lengths (2, 0) differ >" > > # 2. using stars for grid: > > suppressPackageStartupMessages(library(stars)) > st = st_as_stars(meuse.grid) > st_crs(st) Coordinate Reference System: User input: Amersfoort / RD New wkt: PROJCRS["Amersfoort / RD New", BASEGEOGCRS["Amersfoort", DATUM["Amersfoort", ELLIPSOID["Bessel 1841",6377397.155,299.1528128, LENGTHUNIT["metre",1]]], PRIMEM["Greenwich",0, ANGLEUNIT["degree",0.0174532925199433]], ID["EPSG",4289]], CONVERSION["RD New", METHOD["Oblique Stereographic", ID["EPSG",9809]], PARAMETER["Latitude of natural origin",52.1561605555556, ANGLEUNIT["degree",0.0174532925199433], ID["EPSG",8801]], PARAMETER["Longitude of natural origin",5.38763888888889, ANGLEUNIT["degree",0.0174532925199433], ID["EPSG",8802]], PARAMETER["Scale factor at natural origin",0.9999079, SCALEUNIT["unity",1], ID["EPSG",8805]], PARAMETER["False easting",155000, LENGTHUNIT["metre",1], ID["EPSG",8806]], PARAMETER["False northing",463000, LENGTHUNIT["metre",1], ID["EPSG",8807]]], CS[Cartesian,2], AXIS["easting (X)",east, ORDER[1], LENGTHUNIT["metre",1]], AXIS["northing (Y)",north, ORDER[2], LENGTHUNIT["metre",1]], USAGE[ SCOPE["Engineering survey, topographic mapping."], AREA["Netherlands - onshore, including Waddenzee, Dutch Wadden Islands and 12-mile offshore coastal zone."], BBOX[50.75,3.2,53.7,7.22]], ID["EPSG",28992]] > > # compare inputs: > sp = as(st, "Spatial") > fullgrid(meuse.grid) = TRUE > all.equal(sp, meuse.grid["dist"], check.attributes = FALSE) [1] "Names: Lengths (5, 1) differ (string compare on first 1)" [2] "Names: 1 string mismatch" > all.equal(sp, meuse.grid["dist"], check.attributes = TRUE, use.names = FALSE) [1] "Names: Lengths (5, 1) differ (string compare on first 1)" [2] "Names: 1 string mismatch" [3] "Attributes: < Component 3: Names: 1 string mismatch >" [4] "Attributes: < Component 3: Length mismatch: comparison on first 1 components >" [5] "Attributes: < Component 3: Component 1: Mean relative difference: 1.08298 >" [6] "Attributes: < Component 4: Attributes: < Component 2: names for current but not for target > >" [7] "Attributes: < Component 4: Attributes: < Component 3: names for current but not for target > >" > > # kriging: > st_crs(st) = st_crs(meuse_sf) = NA # GDAL roundtrip messes them up! > k_st = if (Sys.getenv("USER") == "travis") { + try(krige(log(zinc)~1, meuse_sf, st, v.fit)) + } else { + krige(log(zinc)~1, meuse_sf, st, v.fit) + } [using ordinary kriging] > k_st stars object with 2 dimensions and 2 attributes attribute(s): Min. 1st Qu. Median Mean 3rd Qu. Max. NA's var1.pred 4.7765547 5.2376293 5.5728839 5.7072287 6.1717619 7.4399911 5009 var1.var 0.0854949 0.1372864 0.1621838 0.1853319 0.2116152 0.5002756 5009 dimension(s): from to offset delta x/y x 1 78 178440 40 [x] y 1 104 333760 -40 [y] > > # handle factors, when going to stars? > k_sp_grd$cls = cut(k_sp_grd$var1.pred, c(0, 5, 6, 7, 8, 9)) > st_as_stars(k_sp_grd) stars object with 2 dimensions and 3 attributes attribute(s): var1.pred var1.var cls Min. :4.777 Min. :0.085 (0,5]: 316 1st Qu.:5.238 1st Qu.:0.137 (5,6]:1778 Median :5.573 Median :0.162 (6,7]: 962 Mean :5.707 Mean :0.185 (7,8]: 47 3rd Qu.:6.172 3rd Qu.:0.212 (8,9]: 0 Max. :7.440 Max. :0.500 NA's :5009 NA's :5009 NA's :5009 dimension(s): from to offset delta refsys x/y x 1 78 178440 40 Amersfoort / RD New [x] y 1 104 333760 -40 Amersfoort / RD New [y] > if (require(raster, quietly = TRUE)) { + print(st_as_stars(raster::stack(k_sp_grd))) # check + print(all.equal(st_redimension(st_as_stars(k_sp_grd)), st_as_stars(raster::stack(k_sp_grd)), check.attributes=FALSE)) + } stars object with 3 dimensions and 1 attribute attribute(s): Min. 1st Qu. Median Mean 3rd Qu. Max. NA's var1.pred 0.0854949 0.2116778 2 2.710347 5.237542 7.439991 15027 dimension(s): from to offset delta refsys values x 1 78 178440 40 Amersfoort / RD New NULL y 1 104 333760 -40 Amersfoort / RD New NULL band 1 3 NA NA NA var1.pred, var1.var , cls x/y x [x] y [y] band [1] TRUE > > suppressPackageStartupMessages(library(spacetime)) > > tm = as.POSIXct("2019-02-25 15:37:24 CET") > n = 4 > s = stars:::st_stars(list(foo = array(1:(n^3), rep(n,3))), + stars:::create_dimensions(list( + x = stars:::create_dimension(from = 1, to = n, offset = 10, delta = 0.5), + y = stars:::create_dimension(from = 1, to = n, offset = 0, delta = -0.7), + time = stars:::create_dimension(values = tm + 1:n)), + raster = stars:::get_raster(dimensions = c("x", "y"))) + ) > s stars object with 3 dimensions and 1 attribute attribute(s): Min. 1st Qu. Median Mean 3rd Qu. Max. foo 1 16.75 32.5 32.5 48.25 64 dimension(s): from to offset delta refsys x/y x 1 4 10 0.5 NA [x] y 1 4 0 -0.7 NA [y] time 1 4 2019-02-25 15:37:25 UTC 1 secs POSIXct > > as.data.frame(s) x y time foo 1 10.25 -0.35 2019-02-25 15:37:25 1 2 10.75 -0.35 2019-02-25 15:37:25 2 3 11.25 -0.35 2019-02-25 15:37:25 3 4 11.75 -0.35 2019-02-25 15:37:25 4 5 10.25 -1.05 2019-02-25 15:37:25 5 6 10.75 -1.05 2019-02-25 15:37:25 6 7 11.25 -1.05 2019-02-25 15:37:25 7 8 11.75 -1.05 2019-02-25 15:37:25 8 9 10.25 -1.75 2019-02-25 15:37:25 9 10 10.75 -1.75 2019-02-25 15:37:25 10 11 11.25 -1.75 2019-02-25 15:37:25 11 12 11.75 -1.75 2019-02-25 15:37:25 12 13 10.25 -2.45 2019-02-25 15:37:25 13 14 10.75 -2.45 2019-02-25 15:37:25 14 15 11.25 -2.45 2019-02-25 15:37:25 15 16 11.75 -2.45 2019-02-25 15:37:25 16 17 10.25 -0.35 2019-02-25 15:37:26 17 18 10.75 -0.35 2019-02-25 15:37:26 18 19 11.25 -0.35 2019-02-25 15:37:26 19 20 11.75 -0.35 2019-02-25 15:37:26 20 21 10.25 -1.05 2019-02-25 15:37:26 21 22 10.75 -1.05 2019-02-25 15:37:26 22 23 11.25 -1.05 2019-02-25 15:37:26 23 24 11.75 -1.05 2019-02-25 15:37:26 24 25 10.25 -1.75 2019-02-25 15:37:26 25 26 10.75 -1.75 2019-02-25 15:37:26 26 27 11.25 -1.75 2019-02-25 15:37:26 27 28 11.75 -1.75 2019-02-25 15:37:26 28 29 10.25 -2.45 2019-02-25 15:37:26 29 30 10.75 -2.45 2019-02-25 15:37:26 30 31 11.25 -2.45 2019-02-25 15:37:26 31 32 11.75 -2.45 2019-02-25 15:37:26 32 33 10.25 -0.35 2019-02-25 15:37:27 33 34 10.75 -0.35 2019-02-25 15:37:27 34 35 11.25 -0.35 2019-02-25 15:37:27 35 36 11.75 -0.35 2019-02-25 15:37:27 36 37 10.25 -1.05 2019-02-25 15:37:27 37 38 10.75 -1.05 2019-02-25 15:37:27 38 39 11.25 -1.05 2019-02-25 15:37:27 39 40 11.75 -1.05 2019-02-25 15:37:27 40 41 10.25 -1.75 2019-02-25 15:37:27 41 42 10.75 -1.75 2019-02-25 15:37:27 42 43 11.25 -1.75 2019-02-25 15:37:27 43 44 11.75 -1.75 2019-02-25 15:37:27 44 45 10.25 -2.45 2019-02-25 15:37:27 45 46 10.75 -2.45 2019-02-25 15:37:27 46 47 11.25 -2.45 2019-02-25 15:37:27 47 48 11.75 -2.45 2019-02-25 15:37:27 48 49 10.25 -0.35 2019-02-25 15:37:28 49 50 10.75 -0.35 2019-02-25 15:37:28 50 51 11.25 -0.35 2019-02-25 15:37:28 51 52 11.75 -0.35 2019-02-25 15:37:28 52 53 10.25 -1.05 2019-02-25 15:37:28 53 54 10.75 -1.05 2019-02-25 15:37:28 54 55 11.25 -1.05 2019-02-25 15:37:28 55 56 11.75 -1.05 2019-02-25 15:37:28 56 57 10.25 -1.75 2019-02-25 15:37:28 57 58 10.75 -1.75 2019-02-25 15:37:28 58 59 11.25 -1.75 2019-02-25 15:37:28 59 60 11.75 -1.75 2019-02-25 15:37:28 60 61 10.25 -2.45 2019-02-25 15:37:28 61 62 10.75 -2.45 2019-02-25 15:37:28 62 63 11.25 -2.45 2019-02-25 15:37:28 63 64 11.75 -2.45 2019-02-25 15:37:28 64 > plot(s, col = sf.colors(), axes = TRUE) > (s.stfdf = as(s, "STFDF")) An object of class "STFDF" Slot "data": foo 1 1 2 2 3 3 4 4 5 5 6 6 7 7 8 8 9 9 10 10 11 11 12 12 13 13 14 14 15 15 16 16 17 17 18 18 19 19 20 20 21 21 22 22 23 23 24 24 25 25 26 26 27 27 28 28 29 29 30 30 31 31 32 32 33 33 34 34 35 35 36 36 37 37 38 38 39 39 40 40 41 41 42 42 43 43 44 44 45 45 46 46 47 47 48 48 49 49 50 50 51 51 52 52 53 53 54 54 55 55 56 56 57 57 58 58 59 59 60 60 61 61 62 62 63 63 64 64 Slot "sp": Object of class SpatialPixels Grid topology: cellcentre.offset cellsize cells.dim x 10.25 0.5 4 y -2.45 0.7 4 SpatialPoints: x y [1,] 10.25 -0.35 [2,] 10.75 -0.35 [3,] 11.25 -0.35 [4,] 11.75 -0.35 [5,] 10.25 -1.05 [6,] 10.75 -1.05 [7,] 11.25 -1.05 [8,] 11.75 -1.05 [9,] 10.25 -1.75 [10,] 10.75 -1.75 [11,] 11.25 -1.75 [12,] 11.75 -1.75 [13,] 10.25 -2.45 [14,] 10.75 -2.45 [15,] 11.25 -2.45 [16,] 11.75 -2.45 Coordinate Reference System (CRS) arguments: NA Slot "time": timeIndex 2019-02-25 15:37:25 1 2019-02-25 15:37:26 2 2019-02-25 15:37:27 3 2019-02-25 15:37:28 4 Slot "endTime": [1] "2019-02-25 15:37:26 UTC" "2019-02-25 15:37:27 UTC" [3] "2019-02-25 15:37:28 UTC" "2019-02-25 15:37:29 UTC" > stplot(s.stfdf, scales = list(draw = TRUE)) > > (s2 = st_as_stars(s.stfdf)) stars object with 3 dimensions and 1 attribute attribute(s): Min. 1st Qu. Median Mean 3rd Qu. Max. foo 1 16.75 32.5 32.5 48.25 64 dimension(s): from to offset delta refsys x/y x 1 4 10 0.5 NA [x] y 1 4 -1.11e-16 -0.7 NA [y] time 1 4 2019-02-25 15:37:25 UTC 1 secs POSIXct > plot(s2, col = sf.colors(), axes = TRUE) > all.equal(s, s2, check.attributes = FALSE) [1] TRUE > > # multiple simulations: > data(meuse, package = "sp") > data(meuse.grid, package = "sp") > coordinates(meuse.grid) <- ~x+y > gridded(meuse.grid) <- TRUE > meuse.grid = st_as_stars(meuse.grid) > meuse_sf = st_as_sf(meuse, coords = c("x", "y")) > g = gstat(NULL, "zinc", zinc~1, meuse_sf, model = vgm(1, "Exp", 300), nmax = 10) > g = gstat(g, "lead", lead~1, meuse_sf, model = vgm(1, "Exp", 300), nmax = 10, fill.cross = TRUE) > set.seed(123) > ## IGNORE_RDIFF_BEGIN > (p = predict(g, meuse.grid, nsim = 5)) drawing 5 multivariate GLS realisations of beta... [using conditional Gaussian simulation] stars object with 3 dimensions and 2 attributes attribute(s): Min. 1st Qu. Median Mean 3rd Qu. Max. NA's zinc 101.7643 199.48523 300.2140 399.7954 544.0229 1817.2057 25045 lead 32.7799 73.01782 107.9339 133.1702 178.8162 638.9396 25045 dimension(s): from to offset delta values x/y x 1 78 178440 40 NULL [x] y 1 104 333760 -40 NULL [y] sample 1 5 NA NA sim1,...,sim5 > ## IGNORE_RDIFF_END > > proc.time() user system elapsed 10.68 0.70 11.37