R Under development (unstable) (2023-10-29 r85433 ucrt) -- "Unsuffered Consequences" Copyright (C) 2023 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. > set.seed(1501) > #----------------------------- > ## IGNORE_RDIFF_BEGIN > library(rtop) > if (interactive()) options(error = recover) > # Read directly from shape-files in data directory > rpath = system.file("extdata",package="rtop") > library(sf) Linking to GEOS 3.11.2, GDAL 3.7.2, PROJ 9.3.0; sf_use_s2() is TRUE > observations = st_read(rpath, "observations") Reading layer `observations' from data source `D:\RCompile\CRANincoming\R-devel\lib\rtop\extdata' using driver `ESRI Shapefile' Simple feature collection with 57 features and 7 fields Geometry type: POLYGON Dimension: XY Bounding box: xmin: 392041.1 ymin: 454983.1 xmax: 512819.5 ymax: 543984.4 Projected CRS: Lambert_Conformal_Conic > predictionLocations = st_read(rpath, "predictionLocations") Reading layer `predictionLocations' from data source `D:\RCompile\CRANincoming\R-devel\lib\rtop\extdata' using driver `ESRI Shapefile' Simple feature collection with 235 features and 5 fields Geometry type: POLYGON Dimension: XY Bounding box: xmin: 393416.8 ymin: 454919.1 xmax: 519903.3 ymax: 543984.4 Projected CRS: Lambert_Conformal_Conic > ## IGNORE_RDIFF_END > observations = as(observations, "Spatial") > predictionLocations = as(predictionLocations, "Spatial") > #Finding a few prediction locations of them > > observations = observations[1:30,] > predictionLocations = predictionLocations[1:2,] > > observations$obs = observations$QSUMMER_OB/observations$AREASQKM > > # Setting some parameters > params = list(gDist = TRUE, cloud = FALSE, rresol = 25, hresol = 3, debug.level = -1) > # Build an object > rtopObj = createRtopObject(observations,predictionLocations, params = params, formulaString = "obs ~ 1" ) > # Fit a variogram (function also creates it) > rtopObj = rtopFitVariogram(rtopObj, iprint = -1) > print(rtopObj$variogramModel, 3) $model [1] "Ex1" $params [1] 2.51e-04 4.26e+05 0.00e+00 1.01e-05 9.23e-01 attr(,"class") [1] "rtopVariogramModel" attr(,"SSErr") [1] 0.557 attr(,"criterion") [1] 0.0981 > #rtopObj = checkVario(rtopObj) > rtopObj2 = rtopKrige(rtopObj, cv = TRUE) [1] "Sampling points from 30 areas" > > > print(attr(rtopObj2$varMatObs,"variogramModel"), 3) $model [1] "Ex1" $params [1] 2.51e-04 4.26e+05 0.00e+00 1.01e-05 9.23e-01 attr(,"class") [1] "rtopVariogramModel" attr(,"SSErr") [1] 0.557 attr(,"criterion") [1] 0.0981 > > rtopObj3 = rtopKrige(rtopObj) [1] "Sampling points from 30 areas" [1] "Sampling points from 2 areas" [1] "Creating prediction semivariance matrix. This can take some time." > > > varmat = varMat(observations, predictionLocations, variogramModel = rtopObj$variogramModel, + gDistEst = TRUE, gDistPred = TRUE, rresol = 25, hresol = 3) [1] "Sampling points from 30 areas" [1] "Sampled on average 54.2 points from 30 areas" [1] "Sampling points from 2 areas" [1] "Sampled on average 77 points from 2 areas" > > all.equal(varmat$varMatObs, rtopObj2$varMatObs) [1] TRUE > rtopObj4 = rtopKrige(rtopObj2) [1] "Sampling points from 2 areas" [1] "Creating prediction semivariance matrix. This can take some time." > > #debug(rtop:::rtopDisc.SpatialPolygons) > # rtopObj5 = rtopKrige(rtopObj, params = list(cnAreas = 5, cDlim = 10, nclus = 2)) > > print(summary(rtopObj2$predictions)) Length Class Mode 30 SpatialPolygonsDataFrame S4 > print(summary(rtopObj3$predictions)) Length Class Mode 2 SpatialPolygonsDataFrame S4 > print(summary(rtopObj4$predictions)) Length Class Mode 2 SpatialPolygonsDataFrame S4 > print(all.equal(rtopObj4$predictions, rtopObj3$predictions)) [1] TRUE > #spplot(rtopObj$predictions,col.regions = bpy.colors(), c("var1.pred","var1.var")) > > # Cross-validation > #spplot(rtopObj2$predictions,col.regions = bpy.colors(), c("observed","var1.pred")) > print(cor(rtopObj2$predictions$observed,rtopObj2$predictions$var1.pred)) [1] 0.1678744 > > > > > set.seed(1501) > library(intamap) Loading required package: sp > useRtopWithIntamap() Loading optional package: intamap > ## IGNORE_RDIFF_BEGIN > output = interpolate(observations,predictionLocations, + optList = list(formulaString = obs~1, gDist = TRUE, cloud = FALSE, nmax = 10, rresol = 25, hresol = 3), + methodName = "rtop", iprint = -1) R 2023-10-30 11:47:03.248372 interpolating 30 observations, 2 prediction locations [1] "standardization of projections is not possible when usergdal = FALSE" [1] "rgdal is about to be retired. After this, some of the checks on projections in the intamap package will disappear" Checking object ... OK [1] "Sampling points from 30 areas" [1] "Sampled on average 54.2 points from 30 areas" [1] "Sampling points from 2 areas" [1] "Sampled on average 77 points from 2 areas" [1] "Creating prediction semivariance matrix. This can take some time." Warning messages: 1: In predictTime(nObs = dim(observations)[1], nPred = nPred, formulaString = formulaString, : using standard model for estimating time. For better platform spesific predictions, please run timeModels <- generateTimeModels() and save the workspace 2: In predictTime(nObs = dim(observations)[1], nPred = nPred, formulaString = formulaString, : Could not find time model for method spatialPredict_rtop 3: In interpolate(observations, predictionLocations, optList = list(formulaString = obs ~ : was not able to estimate prediction time for methodName rtop > ## IGNORE_RDIFF_END > > > print(all.equal(rtopObj4$predictions@data$var1.pred, output$predictions@data$var1.pred)) [1] TRUE > print(all.equal(rtopObj4$predictions@data$var1.var, output$predictions@data$var1.var)) [1] TRUE > > > # Updating variogramModel > > rtopObj5 = varMat(rtopObj4) > rtopObj6 = updateRtopVariogram(rtopObj5, exp = 1.5, action = "mult") > rtopObj7 = varMat(rtopObj6) [1] "Creating prediction semivariance matrix. This can take some time." > > > > > # observations$obs = log(observations$obs) > > # Setting some parameters > # Build an object > rtopObj = createRtopObject(observations,predictionLocations, params = params, formulaString = "obs~1") > # Fit a variogram (function also creates it) > rtopObj = rtopFitVariogram(rtopObj, iprint = -1) > #rtopObj = checkVario(rtopObj) > > rtopObj10 = rtopSim(rtopObj, nsim = 5, logdist = TRUE, debug.level = -1) [1] "Sampling points from 30 areas" [1] "Sampling points from 2 areas" [1] "Creating prediction semivariance matrix. This can take some time." [1] "1. simulation of 2 areas" [1] "2. simulation of 2 areas" [1] "3. simulation of 2 areas" [1] "4. simulation of 2 areas" [1] "5. simulation of 2 areas" > rtopObj11 = rtopObj > rtopObj11$predictionLocations = rtopObj11$observations > #rtopObj11$observations = NULL > rtopObj11$observations$unc = var(rtopObj10$observations$obs)*min(rtopObj10$observations$area)/rtopObj10$observations$area > rtopObj11$predictionLocations$replaceNumber = 1:dim(rtopObj11$predictionLocations)[1] > rtopObj12 = rtopSim(rtopObj11, nsim = 10, replace = TRUE, debug.level = -1) [1] "Sampling points from 30 areas" [1] "Sampling points from 30 areas" [1] "Creating prediction semivariance matrix. This can take some time." [1] "1. simulation of 30 areas" [1] "2. simulation of 30 areas" [1] "3. simulation of 30 areas" [1] "4. simulation of 30 areas" [1] "5. simulation of 30 areas" [1] "6. simulation of 30 areas" [1] "7. simulation of 30 areas" [1] "8. simulation of 30 areas" [1] "9. simulation of 30 areas" [1] "10. simulation of 30 areas" > > print(rtopObj10$simulations@data, digits = 3) ID EZGID AREASQKM XSTATION YSTATION area sim1 sim2 sim3 sim4 1 76 76 36.0 490602 523496 35982138 0.0123 0.0106 0.0112 0.0117 2 77 77 38.5 490602 523496 38487054 0.0124 0.0120 0.0121 0.0133 sim5 1 0.0112 2 0.0112 > print(rtopObj12$simulations@data, digits = 3) ID EZGID AREASQKM XSTATION YSTATION QSUMMER_OB obs area 1 60 60 44.0 444255 519555 0.367 0.00836 4.40e+07 2 113 113 13.6 410740 471559 0.168 0.01228 1.36e+07 3 227 227 37.4 507123 505187 0.424 0.01134 3.74e+07 4 550 550 22.8 456944 464894 0.304 0.01337 2.28e+07 5 688 688 19.8 425108 474095 0.240 0.01212 1.98e+07 6 752 752 27.6 470975 496570 0.294 0.01066 2.76e+07 7 765 765 56.2 434609 491830 0.473 0.00841 5.62e+07 8 849 849 26.1 430218 486103 0.244 0.00935 2.61e+07 9 862 862 14.2 420165 481802 0.164 0.01153 1.42e+07 10 863 863 19.9 426344 482108 0.141 0.00710 1.99e+07 11 864 864 54.2 410554 481291 0.681 0.01257 5.42e+07 12 992 992 29.9 400882 475017 0.481 0.01610 2.99e+07 13 1221 1221 64.3 486654 512168 0.581 0.00904 6.43e+07 14 1222 1222 16.7 441490 535299 0.505 0.03029 1.67e+07 15 1632 1632 53.6 429735 480468 0.541 0.01009 5.36e+07 16 1698 1698 59.5 422167 489280 0.567 0.00953 5.95e+07 17 1939 1939 94.0 486398 524930 1.116 0.01187 9.40e+07 18 1953 1953 61.6 469348 507840 0.549 0.00891 6.16e+07 19 1960 1960 50.1 445002 470803 0.404 0.00805 5.01e+07 20 2062 2062 28.3 469204 507994 0.318 0.01123 2.83e+07 21 2067 2067 69.5 454679 494878 0.401 0.00577 6.95e+07 22 2068 2068 75.9 484016 497119 0.410 0.00540 7.59e+07 23 2196 2196 63.4 421439 501348 0.920 0.01451 6.34e+07 24 2264 2264 82.4 477861 478778 0.810 0.00983 8.24e+07 25 2269 2269 62.3 442000 476361 0.696 0.01117 6.23e+07 26 2280 2280 58.4 412619 458520 1.017 0.01740 5.84e+07 27 2284 2284 52.2 428161 458036 0.667 0.01279 5.22e+07 28 2364 2364 23.5 406824 470502 0.477 0.02032 2.35e+07 29 2384 2384 48.0 431866 497177 0.588 0.01225 4.80e+07 30 2387 2387 127.9 417931 497463 1.300 0.01017 1.28e+08 replaceNumber sim1 sim2 sim3 sim4 sim5 sim6 sim7 1 1 0.01339 0.01285 0.01036 0.01525 0.00879 0.01083 0.01117 2 2 0.01308 0.01459 0.01208 0.01243 0.01307 0.01249 0.01378 3 3 0.01215 0.01265 0.01323 0.01154 0.01184 0.00444 0.01167 4 4 0.01324 0.01059 0.01083 0.00873 0.01441 0.01255 0.01028 5 5 0.01069 0.01001 0.01108 0.01304 0.01149 0.01277 0.01238 6 6 0.01076 0.00914 0.00910 0.00655 0.00942 0.00769 0.00884 7 7 0.00958 0.00985 0.00976 0.01027 0.00835 0.00522 0.00957 8 8 0.01070 0.00791 0.00964 0.01217 0.00739 0.01053 0.00650 9 9 0.01142 0.00951 0.01067 0.01215 0.01054 0.00895 0.01000 10 10 0.01038 0.00954 0.01092 0.01153 0.00933 0.00977 0.00839 11 11 0.01148 0.00980 0.01054 0.01213 0.01192 0.01091 0.01174 12 12 0.01844 0.01517 0.01265 0.01561 0.01995 0.01557 0.01743 13 13 0.00773 0.00991 0.01181 0.00856 0.00658 0.01020 0.00949 14 14 0.01621 0.01723 0.01995 0.02314 0.01963 0.02078 0.02059 15 15 0.01086 0.00939 0.01036 0.01124 0.00953 0.01076 0.01092 16 16 0.00977 0.00792 0.01075 0.00940 0.00916 0.01099 0.00995 17 17 0.00934 0.01326 0.01278 0.01177 0.00936 0.00897 0.01438 18 18 0.00990 0.00895 0.01035 0.00855 0.01014 0.00939 0.01038 19 19 0.01049 0.01087 0.00936 0.01066 0.00767 0.00867 0.01081 20 20 0.01179 0.01055 0.01104 0.00930 0.00835 0.00947 0.01035 21 21 0.00781 0.01019 0.00855 0.00872 0.00642 0.00643 0.00706 22 22 0.00819 0.00882 0.00675 0.00403 0.00665 0.00671 0.00923 23 23 0.01233 0.01366 0.01484 0.01159 0.01052 0.01108 0.01196 24 24 0.00923 0.00914 0.01118 0.01119 0.01111 0.01036 0.01345 25 25 0.00954 0.00898 0.01096 0.01021 0.00843 0.01012 0.01017 26 26 0.01856 0.01691 0.01766 0.01248 0.01687 0.01539 0.01446 27 27 0.01251 0.01291 0.01216 0.01385 0.01317 0.01470 0.01271 28 28 0.01935 0.01488 0.01535 0.01387 0.01877 0.01604 0.01800 29 29 0.01506 0.01130 0.01191 0.01109 0.01028 0.00900 0.01288 30 30 0.00986 0.00939 0.01119 0.01003 0.00852 0.01101 0.01044 sim8 sim9 sim10 1 0.01495 0.00948 0.01526 2 0.01787 0.01606 0.01145 3 0.00863 0.01312 0.01280 4 0.00878 0.00787 0.01165 5 0.01263 0.01070 0.01175 6 0.00627 0.01197 0.01278 7 0.01182 0.01123 0.00954 8 0.01013 0.01084 0.00870 9 0.01319 0.01050 0.01130 10 0.01028 0.00911 0.01124 11 0.01396 0.01395 0.01229 12 0.02081 0.01813 0.01697 13 0.01069 0.00682 0.00961 14 0.02033 0.01979 0.02098 15 0.01153 0.00977 0.01335 16 0.00973 0.01092 0.01103 17 0.01101 0.01013 0.01140 18 0.00939 0.00985 0.01057 19 0.00963 0.00995 0.00957 20 0.01129 0.01052 0.01211 21 0.00660 0.00899 0.00994 22 0.00726 0.00807 0.00686 23 0.01615 0.01404 0.01088 24 0.00972 0.01114 0.00948 25 0.01023 0.00971 0.01288 26 0.01917 0.01580 0.01769 27 0.01058 0.01186 0.01486 28 0.02334 0.01921 0.02014 29 0.01327 0.01246 0.01032 30 0.01140 0.01101 0.00949 > > > > proc.time() user system elapsed 28.06 0.96 29.04