Geographical Convergent Cross Mapping (GCCM)
Wenbo Lv
2025-02-20
Source:vignettes/GCCM.Rmd
GCCM.Rmd1.1 Install the spEDM package
Install the stable version from CRAN with:
install.packages("spEDM", dep = TRUE)Alternatively, you can install the development version from R-universe with:
install.packages("spEDM",
repos = c("https://stscl.r-universe.dev",
"https://cloud.r-project.org"),
dep = TRUE)1.2 An example of spatial lattice data
Load the spEDM package:
Load the county-level population density data from the
spEDM package:
popd_nb = spdep::read.gal(system.file("extdata/popdensity_nb.gal",
package = "spEDM"))
## Warning in spdep::read.gal(system.file("extdata/popdensity_nb.gal", package = "spEDM")): neighbour
## object has 4 sub-graphs
popd_nb
## Neighbour list object:
## Number of regions: 2806
## Number of nonzero links: 15942
## Percentage nonzero weights: 0.2024732
## Average number of links: 5.681397
## 4 disjoint connected subgraphs
popdensity = readr::read_csv(system.file("extdata/popdensity.csv",
package = "spEDM"))
## Rows: 2806 Columns: 7
## ── Column specification ────────────────────────────────────────────────────────────────────────────────
## Delimiter: ","
## dbl (7): x, y, popDensity, DEM, Tem, Pre, slop
##
## ℹ Use `spec()` to retrieve the full column specification for this data.
## ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.
popdensity
## # A tibble: 2,806 × 7
## x y popDensity DEM Tem Pre slop
## <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 117. 30.5 780. 8 17.4 1528. 0.452
## 2 117. 30.6 395. 48 17.2 1487. 0.842
## 3 117. 30.8 261. 49 16.0 1456. 3.56
## 4 116. 30.1 258. 23 17.4 1555. 0.932
## 5 116. 30.5 211. 101 16.3 1494. 3.34
## 6 117. 31.0 386. 10 16.6 1382. 1.65
## 7 117. 30.2 350. 23 17.5 1569. 0.346
## 8 117. 30.7 470. 22 17.1 1493. 1.88
## 9 117. 30.6 1226. 11 17.4 1526. 0.208
## 10 116. 30.9 137. 598 13.9 1458. 5.92
## # ℹ 2,796 more rows
popd_sf = sf::st_as_sf(popdensity, coords = c("x","y"), crs = 4326)
popd_sf
## Simple feature collection with 2806 features and 5 fields
## Geometry type: POINT
## Dimension: XY
## Bounding box: xmin: 74.9055 ymin: 18.2698 xmax: 134.269 ymax: 52.9346
## Geodetic CRS: WGS 84
## # A tibble: 2,806 × 6
## popDensity DEM Tem Pre slop geometry
## * <dbl> <dbl> <dbl> <dbl> <dbl> <POINT [°]>
## 1 780. 8 17.4 1528. 0.452 (116.912 30.4879)
## 2 395. 48 17.2 1487. 0.842 (116.755 30.5877)
## 3 261. 49 16.0 1456. 3.56 (116.541 30.7548)
## 4 258. 23 17.4 1555. 0.932 (116.241 30.104)
## 5 211. 101 16.3 1494. 3.34 (116.173 30.495)
## 6 386. 10 16.6 1382. 1.65 (116.935 30.9839)
## 7 350. 23 17.5 1569. 0.346 (116.677 30.2412)
## 8 470. 22 17.1 1493. 1.88 (117.066 30.6514)
## 9 1226. 11 17.4 1526. 0.208 (117.171 30.5558)
## 10 137. 598 13.9 1458. 5.92 (116.208 30.8983)
## # ℹ 2,796 more rowsSelect the appropriate embedding dimension E:
simplex(popd_sf,"Pre",lib = 1:2000,pred = 2001:nrow(popd_sf),k = 6,nb = popd_nb)
## The suggested embedding dimension E for variable Pre is 2
## E rho mae rmse
## [1,] 1 0.9782781 34.80551 54.00773
## [2,] 2 0.9802540 32.18977 53.33620
## [3,] 3 0.9775425 34.25069 58.26795
## [4,] 4 0.9749178 36.43196 61.47662
## [5,] 5 0.9732294 38.25935 63.20869
## [6,] 6 0.9723361 40.15206 65.04856
## [7,] 7 0.9697997 43.27062 69.01631
## [8,] 8 0.9658433 46.51147 73.37128
## [9,] 9 0.9629774 49.73771 76.06562
## [10,] 10 0.9615392 51.52285 77.57833
simplex(popd_sf,"popDensity",lib = 1:2000,pred = 2001:nrow(popd_sf),k = 6,nb = popd_nb)
## The suggested embedding dimension E for variable popDensity is 5
## E rho mae rmse
## [1,] 1 0.8142262 821.6575 2283.544
## [2,] 2 0.8799178 648.7316 1918.029
## [3,] 3 0.8828657 636.7622 1881.740
## [4,] 4 0.8868719 637.2017 1850.598
## [5,] 5 0.8956243 610.3778 1805.386
## [6,] 6 0.8939081 640.9628 1819.206
## [7,] 7 0.8926276 626.0359 1831.724
## [8,] 8 0.8894480 645.2201 1851.623
## [9,] 9 0.8846390 660.4611 1880.918
## [10,] 10 0.8846890 659.1725 1876.446We choose the E with the highest rho and the lowest MAE and RMSE as
the most suitable one. Under the selected lib and pred, the optimal
embedding dimension E for the variable Pre is 2, and for
the variable popDensity, it is 5.
Then, run GCCM:
startTime = Sys.time()
pd_res = gccm(data = popd_sf,
cause = "Pre",
effect = "popDensity",
libsizes = seq(10, 2800, by = 100),
E = c(2,5),
k = 6,
nb = popd_nb,
progressbar = FALSE)
endTime = Sys.time()
print(difftime(endTime,startTime, units ="mins"))
## Time difference of 23.05816 mins
pd_res
## libsizes Pre->popDensity popDensity->Pre
## 1 10 0.05589905 0.01004636
## 2 110 0.18109636 0.03478543
## 3 210 0.24325566 0.04932428
## 4 310 0.27657692 0.06417637
## 5 410 0.30630281 0.07634983
## 6 510 0.34116169 0.08815848
## 7 610 0.38078810 0.09890353
## 8 710 0.41612716 0.10867846
## 9 810 0.44420437 0.11714337
## 10 910 0.46815429 0.12485229
## 11 1010 0.48876449 0.13194552
## 12 1110 0.50538462 0.13802611
## 13 1210 0.52300459 0.14335037
## 14 1310 0.53902632 0.14813342
## 15 1410 0.55309048 0.15284706
## 16 1510 0.56754023 0.15777314
## 17 1610 0.58387624 0.16301596
## 18 1710 0.60029184 0.16768202
## 19 1810 0.61631190 0.17202960
## 20 1910 0.63213899 0.17601166
## 21 2010 0.64785899 0.17979878
## 22 2110 0.66289302 0.18304365
## 23 2210 0.67740971 0.18596497
## 24 2310 0.69156193 0.18906187
## 25 2410 0.70541793 0.19196548
## 26 2510 0.71925253 0.19503355
## 27 2610 0.73287448 0.19819076
## 28 2710 0.74435871 0.20116518Visualize the result:
Figure 1. The cross-mapping
prediction outputs between population density and county-level
Precipitation.
1.3 An example of spatial grid data
Load the spEDM package:
Load the farmland NPP data from the spEDM package:
npp = terra::rast(system.file("extdata/npp.tif", package = "spEDM"))
npp
## class : SpatRaster
## dimensions : 404, 483, 3 (nrow, ncol, nlyr)
## resolution : 10000, 10000 (x, y)
## extent : -2625763, 2204237, 1877078, 5917078 (xmin, xmax, ymin, ymax)
## coord. ref. : CGCS2000_Albers
## source : npp.tif
## names : npp, pre, tem
## min values : 164.00, 384.3409, -47.8194
## max values : 16606.33, 23878.3555, 263.6938
terra::plot(npp, nc = 3,
mar = rep(0.1,4),
oma = rep(0.1,4),
axes = FALSE,
legend = FALSE)
Figure 2. Maps of farmland NPP
and climate factors.
To save the computation time, we will aggregate the data by 3 times and select 1500 non-NA pixels to predict:
npp = terra::aggregate(npp, fact = 3, na.rm = TRUE)
terra::global(npp,"isNA")
## isNA
## npp 14815
## pre 14766
## tem 14766
terra::ncell(npp)
## [1] 21735
nnamat = terra::as.matrix(!is.na(npp[[1]]), wide = TRUE)
nnaindice = terra::rowColFromCell(npp,which(nnamat))
dim(nnaindice)
## [1] 6920 2
set.seed(42)
indices = sample(nrow(nnaindice), size = 1500, replace = FALSE)
lib = nnaindice[-indices,]
pred = nnaindice[indices,]Due to the high number of NA values in the npp raster data, we used all non-NA cell as the libraries when testing for the most suitable embedding dimension.
simplex(npp,"pre",nnaindice,pred,k = 5)
## The suggested embedding dimension E for variable pre is 2
## E rho mae rmse
## [1,] 1 0.9935452 192.2437 257.0223
## [2,] 2 0.9964007 138.1465 192.6528
## [3,] 3 0.9933468 152.9612 261.6461
## [4,] 4 0.9947055 152.8392 232.8033
## [5,] 5 0.9946261 159.5439 234.5342
## [6,] 6 0.9945737 161.3816 235.8791
## [7,] 7 0.9944200 166.2916 239.2797
## [8,] 8 0.9944050 169.3894 239.9678
## [9,] 9 0.9941899 170.0724 244.2513
## [10,] 10 0.9947184 167.4744 233.1708
simplex(npp,"npp",nnaindice,pred,k = 5)
## The suggested embedding dimension E for variable npp is 6
## E rho mae rmse
## [1,] 1 0.9309598 424.9317 608.9018
## [2,] 2 0.9361632 382.7654 583.7951
## [3,] 3 0.9376302 374.3546 575.0216
## [4,] 4 0.9406945 370.6129 561.0225
## [5,] 5 0.9466347 349.3912 532.1607
## [6,] 6 0.9497837 340.4782 516.3542
## [7,] 7 0.9488605 341.0529 520.8572
## [8,] 8 0.9461863 346.1658 533.9391
## [9,] 9 0.9448451 344.5612 540.3860
## [10,] 10 0.9458182 341.8827 535.7265Under the selected lib and pred, the optimal embedding dimension E
for the variable pre is 2, and for the variable
npp, it is 6.
Then, run GCCM:
startTime = Sys.time()
npp_res = gccm(data = npp,
cause = "pre",
effect = "npp",
libsizes = seq(10,100,5),
E = c(2,6),
k = 5,
pred = pred,
progressbar = FALSE)
endTime = Sys.time()
print(difftime(endTime,startTime, units ="mins"))
## Time difference of 21.55607 mins
npp_res
## libsizes pre->npp npp->pre
## 1 10 0.1034609 0.06800981
## 2 15 0.1249670 0.08183160
## 3 20 0.1568069 0.09673624
## 4 25 0.1850738 0.11198871
## 5 30 0.2077103 0.11475138
## 6 35 0.2364462 0.12451989
## 7 40 0.2688678 0.13374827
## 8 45 0.3029108 0.14820557
## 9 50 0.3436574 0.16371411
## 10 55 0.4011195 0.19156188
## 11 60 0.4471125 0.21612723
## 12 65 0.4741198 0.24044547
## 13 70 0.4871746 0.24450378
## 14 75 0.5123036 0.23758264
## 15 80 0.5634948 0.19366801
## 16 85 0.5849654 0.15820213
## 17 90 0.5941501 0.16207449
## 18 95 0.6568359 0.16705620
## 19 100 0.8210545 0.13783912Visualize the result:
plot(npp_res,xlimits = c(9, 101),ylimits = c(-0.05,1)) +
ggplot2::theme(legend.justification = c(0.95,1))
Figure 3. The cross-mapping
prediction outputs between farmland NPP and precipitation.