-
Notifications
You must be signed in to change notification settings - Fork 2
/
spatial_indices.R
812 lines (720 loc) · 27.5 KB
/
spatial_indices.R
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
598
599
600
601
602
603
604
605
606
607
608
609
610
611
612
613
614
615
616
617
618
619
620
621
622
623
624
625
626
627
628
629
630
631
632
633
634
635
636
637
638
639
640
641
642
643
644
645
646
647
648
649
650
651
652
653
654
655
656
657
658
659
660
661
662
663
664
665
666
667
668
669
670
671
672
673
674
675
676
677
678
679
680
681
682
683
684
685
686
687
688
689
690
691
692
693
694
695
696
697
698
699
700
701
702
703
704
705
706
707
708
709
710
711
712
713
714
715
716
717
718
719
720
721
722
723
724
725
726
727
728
729
730
731
732
733
734
735
736
737
738
739
740
741
742
743
744
745
746
747
748
749
750
751
752
753
754
755
756
757
758
759
760
761
762
763
764
765
766
767
768
769
770
771
772
773
774
775
776
777
778
779
780
781
782
783
784
785
786
787
788
789
790
791
792
793
794
795
796
797
798
799
800
801
802
803
804
805
806
807
808
809
810
811
812
#%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
#### spatial consistency index ####
#%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
#' @title Spatial consistency index
#'
#' @description Calculate a spatial consistency index
#'
#' @details This index is experimental, it aims to measure how much a clustering solution
#' is spatially consistent. A classification is spatially inconsistent if
#' neighbouring observation do not belong to the same group. See detail for
#' a description of its calculation
#'
#' The total spatial inconsistency (*Scr*) is calculated as follow
#'
#' \deqn{isp = \sum_{i}\sum_{j}\sum_{k} (u_{ik} - u_{jk})^{2} * W_{ij}}
#'
#' With U the membership matrix, i an observation, k the neighbours of i and W
#' the spatial weight matrix This represents the total spatial inconsistency of
#' the solution (true inconsistency) We propose to compare this total with
#' simulated values obtained by permutations (simulated inconsistency). The
#' values obtained by permutation are an approximation of the spatial
#' inconsistency obtained in a random context Ratios between the true
#' inconsistency and simulated inconsistencies are calculated A value of 0
#' depict a situation where all observations are identical to their neighbours
#' A value of 1 depict a situation where all observations are as much different
#' as their neighbours that what randomness can produce A classification
#' solution able to reduce this index has a better spatial consistency
#'
#' @template FCMresobj-arg
#' @template nblistw2-arg
#' @param window if rasters were used for the classification, the window must be
#' specified instead of a list.w object. Can also be NULL if object is a FCMres object.
#' @param nrep An integer indicating the number of permutation to do to simulate
#' spatial randomness. Note that if rasters are used, each permutation can be very long.
#' @param adj A boolean indicating if the adjusted version of the indicator must be
#' calculated when working with rasters (globally standardized). When working with vectors, see the function
#' adjustSpatialWeights to modify the list.w object.
#' @param mindist When adj is true, a minimum value for distance between two observations. If two
#' neighbours have exactly the same values, then the euclidean distance between
#' them is 0, leading to an infinite spatial weight. In that case, the minimum
#' distance is used instead of 0.
#' @return A named list with
#' \itemize{
#' \item Mean : the mean of the spatial consistency index
#' \item prt05 : the 5th percentile of the spatial consistency index
#' \item prt95 : the 95th percentile of the spatial consistency index
#' \item samples : all the value of the spatial consistency index
#' \item sum_diff : the total sum of squarred difference between observations and their neighbours
#' }
#' @export
#' @examples
#' data(LyonIris)
#' AnalysisFields <-c("Lden","NO2","PM25","VegHautPrt","Pct0_14","Pct_65","Pct_Img",
#' "TxChom1564","Pct_brevet","NivVieMed")
#' dataset <- sf::st_drop_geometry(LyonIris[AnalysisFields])
#' queen <- spdep::poly2nb(LyonIris,queen=TRUE)
#' Wqueen <- spdep::nb2listw(queen,style="W")
#' result <- SFCMeans(dataset, Wqueen,k = 5, m = 1.5, alpha = 1.5, standardize = TRUE)
#' # NOTE : more replications are needed for proper inference
#' spConsistency(result$Belongings, nblistw = Wqueen, nrep=25)
spConsistency <- function(object, nblistw = NULL, window = NULL, nrep = 999, adj = FALSE, mindist = 1e-11) {
if(inherits(object, "FCMres")){
belongmat <- as.matrix(object$Belongings)
if(object$isRaster & is.null(window)){
window <- object$window
if(is.null(window)){
stop("impossible to find a window in the given object, please
specify one by hand.")
}
}
if(object$isRaster == FALSE & is.null(nblistw)){
nblistw <- object$nblistw
}
}else{
belongmat <- as.matrix(object)
}
# if we are not in raster mode
if(is.null(window)){
if(is.null(nblistw)){
stop("The nblistw must be provided if spatial vector data is used")
}
weights <- nblistw$weights
neighbours <- nblistw$neighbours
## calcul de l'inconsistence spatiale actuelle
obsdev <- sapply(1:nrow(belongmat), function(i) {
row <- belongmat[i, ]
idneighbour <- neighbours[[i]]
neighbour <- belongmat[idneighbour, ]
if (length(idneighbour) == 1){
neighbour <- t(as.matrix(neighbour))
}
W <- weights[[i]]
diff <- (neighbour-row[col(neighbour)])**2 * W
tot <- sum(rowSums(diff))
return(tot)
})
totalcons <- sum(obsdev)
## simulation de l'inconsistance spatiale
belongmat <- t(belongmat)
n <- ncol(belongmat)
simulated <- vapply(1:nrep, function(d) {
belong2 <- belongmat[,sample(n)]
simvalues <- vapply(1:ncol(belong2), function(i) {
row <- belong2[,i]
idneighbour <- neighbours[[i]]
neighbour <- belong2[,neighbours[[i]]]
if (length(idneighbour) == 1){
neighbour <- t(as.matrix(neighbour))
}
W <- weights[[i]]
diff <- (neighbour-row)
tot <- sum(diff^2 * W)
return(tot)
}, FUN.VALUE = 1)
return(sum(simvalues))
},FUN.VALUE = 1)
ratio <- totalcons / simulated
# if we are using a raster mode.
}else{
# we must calculate for each pixel its distance to its neighbours
# on the membership matrix. So we will calculate the distance for each
# raster in object$rasters and then sum them all group
rastnames <- names(object$rasters)
ok_names <- rastnames[grepl("group",rastnames, fixed = TRUE)]
rasters <- object$rasters[ok_names]
matrices <- lapply(rasters, terra::as.matrix, wide = TRUE)
mat_dim <- dim(matrices[[1]])
if(adj){
dataset <- lapply(1:ncol(object$Data), function(ic){
vec1 <- object$Data[,ic]
vec2 <- rep(NA,length(object$missing))
vec2[object$missing] <- vec1
rast <- object$rasters[[1]]
terra::values(rast) <- vec2
mat <- terra::as.matrix(rast, wide = TRUE)
return(mat)
})
totalcons <- calc_raster_spinconsistency(matrices, window, adj, dataset, mindist = mindist)
}else{
totalcons <- calc_raster_spinconsistency(matrices,window)
}
# we must now do the same thing but with resampled values
warning("Calculating the permutation for the spatial inconsistency
when using raster can be long, depending on the raster size.
Note that the high number of cell in a raster reduces the need of
a great number of replications.")
# creating a vector of ids for each cell in raster
all_ids <- 1:terra::ncell(rasters[[1]])
# converting the matrices (columns of membership matrix) into 1d vectors
mem_vecs <- lapply(rasters, function(rast){
mat <- terra::as.matrix(rast, wide = TRUE)
dim(mat) <- NULL
return(mat)
})
# if necessary, doing the same with the original data
if(adj){
data_vecs <- lapply(dataset, function(mat){
vec <- mat
dim(vec) <- NULL
return(vec)
})
}
# extracting the dimension of the raster
#dim(terra::as.matrix(rasters[[1]]))
# starting the simulations
simulated <- sapply(1:nrep, function(i){
# resampling the ids
Ids <- sample(all_ids)
# resampling the matrices of memberships
new_matrices <- lapply(mem_vecs, function(vec){
new_vec <- vec[Ids]
# swapping the NA at their original place
val_na <- new_vec[!object$missing]
loc_na <- is.na(new_vec)
new_vec[!object$missing] <- NA
new_vec[loc_na] <- val_na
dim(new_vec) <- mat_dim
return(new_vec)
})
if(adj){
# resampling the matrices of the data
new_dataset <- lapply(data_vecs, function(vec){
new_vec <- vec[Ids];
# swapping the NA at their original place
val_na <- new_vec[!object$missing]
loc_na <- is.na(new_vec)
new_vec[!object$missing] <- NA
new_vec[loc_na] <- val_na
dim(new_vec) <- mat_dim
return(new_vec)
})
# calculating the index value
inconsist <- calc_raster_spinconsistency(new_matrices, window, adj, new_dataset, mindist = mindist)
}else{
# calculating the index value
inconsist <- calc_raster_spinconsistency(new_matrices, window)
}
return(inconsist)
})
ratio <- totalcons / simulated
}
return(list(Mean = mean(ratio), Median = quantile(ratio, probs = c(0.5)),
prt05 = quantile(ratio, probs = c(0.05)),
prt95 = quantile(ratio, probs = c(0.95)),
samples = ratio,
sum_diff = totalcons))
}
#' @title calculate spatial inconsistency for raster
#'
#' @description Calculate the spatial inconsistency sum for a set of rasters
#'
#' @param matrices A list of matrices
#' @param window The window to use to define spatial neighbouring
#' @param adj A boolean indicating if the adjusted version of the algorithm must be
#' calculated
#' @param dataset A list of matrices with the original data (if adj = TRUE)
#' @param mindist When adj is true, a minimum value for distance between two observations. If two
#' neighbours have exactly the same values, then the euclidean distance between
#' them is 0, leading to an infinite spatial weight. In that case, the minimum
#' distance is used instead of 0.
#' @return A float: the sum of spatial inconsistency
#' @keywords internal
#' @examples
#' # this is an internal function, no example provided
calc_raster_spinconsistency <- function(matrices, window, adj = FALSE, dataset = NULL, mindist = 1e-11){
if(adj & is.null(dataset)){
stop("When the adjusted version of spinconsistency is required, dataset must be given")
}
arr <- array(do.call(c,matrices), c(nrow(matrices[[1]]), ncol(matrices[[1]]), length(matrices)))
if(adj){
arr2 <- array(do.call(c,dataset), c(nrow(dataset[[1]]), ncol(dataset[[1]]), length(dataset)))
totalcons <- adj_spconsist_arr_window_globstd(arr2, arr, window, mindist = mindist)
}else{
totalcons <- sum(focal_euclidean_arr_window(arr,window), na.rm = TRUE)
}
return(totalcons)
}
#%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
#### ELSA for HARD PARTITION ####
#%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
#' @title Check validity of a dissimilarity matrix
#'
#' @description Check the validity of a dissimilarity matrix
#'
#' @param matdist A dissimilarity matrix
#' @keywords internal
#' @examples
#' # this is an internal function, no example provided
check_matdist <- function(matdist){
if(!isSymmetric(matdist)){
stop("matdist must be a symetric matrix")
}
if(sum(diag(matdist)) != 0){
stop("matdist must have a diagonal filled with zeros")
}
}
#' @title calculate ELSA statistic for a hard partition
#'
#' @description Calculate ELSA statistic for a hard partition. This local indicator of
#' spatial autocorrelation can be used to determine where observations belong to different
#' clusters.
#'
#' @details The ELSA index \insertCite{naimi2019elsa}{geocmeans} can be used to measure
#' local autocorrelation for a categorical variable. It varies between 0 and 1, 0 indicating
#' a perfect positive spatial autocorrelation and 1 a perfect heterogeneity. It is based on
#' the Shanon entropy index, and uses a measure of difference between categories. Thus it
#' can reflect that proximity of two similar categories is still a form of positive
#' autocorelation. The authors suggest to calculate the mean of the index at several lag
#' distance to create an entrogram which quantifies global spatial structure and can be
#' represented as a variogram-like graph.
#'
#' @param object A FCMres object, typically obtained from functions CMeans,
#' GCMeans, SFCMeans, SGFCMeans. Can also be a vector of categories. This vector must
#' be filled with integers starting from 1. -1 can be used to indicate missing categories.
#' @template nblistw-arg
#' @template window2-arg
#' @template matdist-arg
#' @return A depending of the input, a vector of ELSA values or a raster with the ELSA values.
#' @examples
#' data(LyonIris)
#' AnalysisFields <-c("Lden","NO2","PM25","VegHautPrt","Pct0_14","Pct_65","Pct_Img",
#' "TxChom1564","Pct_brevet","NivVieMed")
#' dataset <- sf::st_drop_geometry(LyonIris[AnalysisFields])
#' queen <- spdep::poly2nb(LyonIris,queen=TRUE)
#' Wqueen <- spdep::nb2listw(queen,style="W")
#' result <- SFCMeans(dataset, Wqueen,k = 5, m = 1.5, alpha = 1.5, standardize = TRUE)
#' elsa_valus <- calcELSA(result)
calcELSA <- function(object, nblistw = NULL, window = NULL, matdist = NULL){
# testing if we have all the required parameters
if(inherits(object, "FCMres") == FALSE){
if(inherits(object, "numeric") == FALSE){
stop("if object is not a FCMres object, it must be a numeric vector")
}
vec <- object[object != -1]
if(min(vec) != 1){
stop("if object is a numeric vector, its lower value must be 1")
}
if(length(unique(vec)) != length(1:max(vec))){
stop("if object is a numeric vector, its values must be like 1,2,3,4,...m, -1 can be used to indicate missing values")
}
if(is.null(matdist)){
stop("if object is not a FCMres object, matdist must be provided")
}else{
check_matdist(matdist)
}
if(nrow(matdist) != length(unique(vec))){
stop("the dimension of matdist must equal the number of categories in object.")
}
if(is.null(nblistw)){
stop("if object is not a FCMres object, nblistw must be provided.")
}
}else{
if(object$isRaster){
if(is.null(object$window) & is.null(window)){
stop("impossible to extract window from object, window must be given.")
}
}else{
if(is.null(object$nblistw) & is.null(nblistw)){
stop("impossible to extract nblistw from object, nblistw must be given.")
}
}
}
# case 1 : object is a simple vector of categories
if(inherits(object,'FCMres') == FALSE){
vals <- elsa_vector(object, nblistw, matdist)
}else{
# case 2 : the object is a FCMres object
if(is.null(matdist)){
matdist <- as.matrix(stats::dist(object$Centers))
}
if(object$isRaster){
if(is.null(window)){
window <- object$window
}
vals <- elsa_raster(object$rasters$Groups, window, matdist)
}else{
vec <- as.numeric(gsub("V","",object$Groups,fixed = TRUE))
if(is.null(nblistw)){
vals <- elsa_vector(vec, object$nblistw, matdist)
}else{
vals <- elsa_vector(vec, nblistw, matdist)
}
}
}
return(vals)
}
#' @title calculate ELSA spatial statistic for vector dataset
#'
#' @description calculate ELSA spatial statistic for vector dataset
#'
#' @param categories An integer vector representing the m categories (1,2,3,..., m),
#' -1 is used to indicate missing values.
#' @param nblistw A listw object from spdep representing neighbour relations
#' @param dist A numeric matrix (m*m) representing the distances between categories
#' @return A vector: the local values of ELSA
#' @keywords internal
#' @examples
#' # this is an internal function, no example provided
elsa_vector <- function(categories, nblistw, dist){
d <- max(dist)
m <- length(unique(categories))
values <- sapply(1:length(categories), function(i){
xi <- categories[[i]]
if(xi == -1){
return(-1)
}else{
neighbours <- nblistw$neighbours[[i]]
w <- nblistw$weights[[i]]
xjs <- categories[neighbours]
Eai <- sum(dist[xi,xjs] * w) / (d * sum(w))
nn <- length(w) # the number of observations in the area...
if(nn > m){
mi <- m
}else{
mi <- nn
}
probs <- table(c(xjs,xi)) / (length(xjs)+1)
probs <- probs[probs > 0]
if(mi == 1){
return(0)
}else{
Eci <- -1 * (sum(probs * log2(probs)) / log2(mi))
return(Eai * Eci)
}
}
})
return(values)
}
#' @title calculate ELSA spatial statistic for raster dataset
#'
#' @description calculate ELSA spatial statistic for vector dataset
#'
#' @param rast An integer raster or matrix representing the m categories (0,1,2,..., m)
#' @template window2-arg
#' @template matdist-arg
#' @return A raster or a matrix: the local values of ELSA
#' @keywords internal
#' @examples
#' # this is an internal function, no example provided
elsa_raster <- function(rast, window, matdist){
if(inherits(window, "matrix")){
u <- unique(window)
if(length(union(c(0,1),u)) > 2){
stop("The provided matrix to calculate ELSA is not a binary matrix (0,1)")
}
fun <- Elsa_categorical_matrix_window
}else{
stop("for calculating ELSA on raster, window must be a binary matrix")
}
isRaster <- inherits(rast, "SpatRaster")
if(isRaster){
mat <- terra::as.matrix(rast, wide = TRUE)
}else{
mat <- rast
}
mat <- ifelse(is.na(mat),-1,mat)
# let us check if the lowest value beside -1 is 0
vec <- c(mat)
m <- min(vec[vec!=-1])
if(m > 0){
mat<- ifelse(mat > 0, mat-1,mat)
}
cats <- unique(c(mat))
refcats <- 0:max(cats)
cats <- cats[cats != -1]
cats <- cats[order(cats)]
if(sum(refcats - cats)!=0){
stop(paste("the values of the raster used for ELSA calculation must be integers starting from 0 (or 1).
There must be no jumps between categories. The categories in the actual raster are : ", paste(cats,collapse = ","),sep=""))
}
mat2 <- fun(mat, window, matdist)
mat2 <- ifelse(mat2 == -1, NA, mat2)
if(isRaster){
terra::values(rast) <- mat2
return(rast)
}else{
return(mat2)
}
}
#%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
#### ELSA for fuzzy PARTITION ####
#%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
#' @title calculate ELSA statistic for a fuzzy partition
#'
#' @description Calculate ELSA statistic for a fuzzy partition. This local indicator of
#' spatial autocorrelation can be used to identify areas where close observations tend to
#' belong to different clusters.
#'
#' @details The fuzzy ELSA index is a generalization of the ELSA index \insertCite{naimi2019elsa}{geocmeans}. It can be used to measure
#' local autocorrelation for a membership matrix. It varies between 0 and 1, 0 indicating
#' a perfect positive spatial autocorrelation and 1 a perfect heterogeneity. It is based on
#' the Shannon entropy index, and uses a measure of dissimilarity between categories.
#'
#' @param object A FCMres object, typically obtained from functions CMeans,
#' GCMeans, SFCMeans, SGFCMeans. Can also be a membership matrix. Each row of this matrix
#' must sum up to 1. Can also be a list of rasters, in which case each raster must represent
#' the membership values for one cluster and the sum of all the rasters must be a raster filled
#' with ones.
#' @template nblistw-arg
#' @template window2-arg
#' @template matdist-arg
#' @return either a vector or a raster with the ELSA values.
#' @examples
#' data(LyonIris)
#' AnalysisFields <-c("Lden","NO2","PM25","VegHautPrt","Pct0_14","Pct_65","Pct_Img",
#' "TxChom1564","Pct_brevet","NivVieMed")
#' dataset <- sf::st_drop_geometry(LyonIris[AnalysisFields])
#' queen <- spdep::poly2nb(LyonIris,queen=TRUE)
#' Wqueen <- spdep::nb2listw(queen,style="W")
#' result <- SFCMeans(dataset, Wqueen,k = 5, m = 1.5, alpha = 1.5, standardize = TRUE)
#' elsa_valus <- calcFuzzyELSA(result)
calcFuzzyELSA <- function(object, nblistw = NULL, window = NULL, matdist = NULL){
# testing if we have all the required parameters
if(inherits(object, c("FCMres","matrix","list")) == FALSE){
stop("object must be a FCMres object, a matrix or a list of rasters")
}
if(inherits(object,"matrix")){
sums <- rowSums(object) != 1
if(any(sums)){
stop("if object is a matrix, the sum of each row must be 1")
}
}
if(inherits(object,"FCMres") == FALSE & is.null(matdist)){
stop("if object is not a FCMre, matdist must be specified")
}
if(inherits(object,"matrix") & is.null(nblistw)){
stop("if object is a matrix, nblistw must be specified")
}
if(inherits(object,"list") & inherits(object,"list") == FALSE & is.null(window)){
stop("if object is a list, window must be specified")
}
if(is.null(matdist)==FALSE){
if(!isSymmetric(matdist)){
stop("matdist must be a symetric matrix")
}
}
if(inherits(object,"matrix")){
if(nrow(matdist) != ncol(object)){
stop("the number of columns in object (matrix) must match the number of rows in matdist")
}
if(is.null(nblistw)){
stop("if object is a matrix, nblistw must be provided")
}
}
if(inherits(object,"list") & (inherits(object,"FCMres") == FALSE)){
if(nrow(matdist) != length(object)){
stop("the number of rasters in object (list) must match the number of rows in matdist")
}
if(is.null(window)){
stop("if object is a list of rasters, window must be provided")
}
}
if(inherits(object,"FCMres")){
if(object$isRaster){
if(is.null(object$window) & is.null(window)){
stop("impossible to extract window from object, window must be provided")
}
}else{
if(is.null(object$nblistw) & is.null(nblistw)){
stop("impossible to extract nblistw from object, nblistw must be provided")
}
}
}
if(inherits(object,"FCMres")){
if(object$isRaster==FALSE){
#case 1 : object is a FCMres and vector type
mat <- object$Belongings
if(is.null(matdist)){
matdist <- as.matrix(stats::dist(object$Centers))
}
if(is.null(nblistw)){
nblistw <- object$nblistw
}
return(elsa_fuzzy_vector(mat, nblistw, matdist))
}else{
#case 2 : object is a FCMres and raster type
if(is.null(matdist)){
matdist <- as.matrix(stats::dist(object$Centers))
}
if(is.null(window)){
window <- object$window
}
mats <- lapply(object$rasters[1:object$k], terra::as.matrix, wide = TRUE)
arr <- do.call(c,mats)
arr <- array(arr, dim = c(nrow(mats[[1]]), ncol(mats[[1]]),object$k))
values <- Elsa_fuzzy_matrix_window(arr, window, matdist)
rast <- object$rasters[[1]]
terra::values(rast) <- values
return(rast)
}
}else if (inherits(object,"matrix")){
#case 3 : object is a matrix
values <- elsa_fuzzy_vector(object, nblistw, matdist)
return(values)
}else if (inherits(object,"list")){
#case 4 : object is a list of rasters
mats <- lapply(object, terra::as.matrix, wide = TRUE)
arr <- do.call(c,mats)
arr <- array(arr, dim = c(nrow(mats[[1]]), ncol(mats[[1]]), length(object)))
values <- Elsa_fuzzy_matrix_window(arr, window, matdist)
rast <- object[[1]]
terra::values(rast) <- values
return(rast)
}else{
stop("invalid arguments passed to calcFuzzyELSA")
}
return(values)
}
#' @title Local Fuzzy ELSA statistic for vector
#'
#' @description Calculate the Local Fuzzy ELSA statistic using a nblistw object
#'
#' @param memberships A membership matrix
#' @param nblistw The spatial weight matrix (nblistw object from spdep)
#' @template matdist-arg
#' @return A vector of local ELSA values
#' @keywords internal
#' @examples
#' # this is an internal function, no example provided
elsa_fuzzy_vector <- function(memberships, nblistw, matdist){
d <- max(matdist)
m <- ncol(memberships)
values <- sapply(1:nrow(memberships), function(i){
xi <- memberships[i,]
neighbours <- nblistw$neighbours[[i]]
xj <- memberships[neighbours,]
diffs <- abs(xi - t(xj))
sis <- apply(diffs, 2, function(x){
out <- outer(x,x)
return(sum(out * matdist))
})
s1 <- sum(sis)/2.0
nn <- length(neighbours)
eai <- s1 / (d*nn)
pks <- (colSums2(xj) + xi) / (nn+1)
pks <- pks[pks>0]
if(nn > m){
mi <- m
}else{
mi <- nn
}
eci <- -1 * (sum(pks * log2(pks)) / log2(mi))
return(eai * eci)
})
return(values)
}
#' @title Local Fuzzy ELSA statistic for raster
#'
#' @description Calculate the Local Fuzzy ELSA statistic for a numeric raster
#'
#' @param rasters A List of SpatRaster or a List of matrices, or an array
#' @template window2-arg
#' @template matdist-arg
#' @return A raster or a matrix (depending on the input): the values of
#' local fuzzy ELSA statistic
#' @keywords internal
#' @examples
#' # this is an internal function, no example provided
calcFuzzyElsa_raster <- function(rasters,window,matdist){
if(inherits(rasters[[1]], "matrix")){
mats <- rasters
isRaster <- FALSE
vals <- do.call(c,mats)
arr <- array(vals, c(nrow(mats[[1]]),ncol(mats[[1]]),length(rasters)))
}else if (inherits(rasters[[1]], "SpatRaster")){
mats <- lapply(rasters, terra::as.matrix, wide = TRUE)
isRaster <- TRUE
vals <- do.call(c,mats)
arr <- array(vals, c(nrow(mats[[1]]),ncol(mats[[1]]),length(rasters)))
}else if (inherits(rasters,"array")){
arr <- rasters
isRaster <- FALSE
}else{
stop("rasters must be a list of matrix or a list of SpatRaster or an array")
}
elsa <- Elsa_fuzzy_matrix_window(arr, window, matdist)
if(isRaster){
rast <- rasters[[1]]
terra::values(rast) <- elsa
}else{
dims <- dim(arr)
rast <- matrix(elsa, nrow = dims[[1]], ncol = dims[[2]])
}
return(rast)
}
#%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
#### Moran for rasters ####
#%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
#' @title Global Moran I for raster
#'
#' @description Calculate the global Moran I for a numeric raster
#'
#' @param rast A SpatRaster or a matrix
#' @param window The window defining the neighbour weights
#' @return A float: the global Moran I
#' @examples
#' Arcachon <- terra::rast(system.file("extdata/Littoral4_2154.tif", package = "geocmeans"))
#' names(Arcachon) <- c("blue", "green", "red", "infrared", "SWIR1", "SWIR2")
#' rast <- Arcachon[[1]]
#' w <- matrix(1, nrow = 3, ncol = 3)
#' calc_moran_raster(rast, w)
calc_moran_raster <- function(rast, window){
if(inherits(rast,"SpatRaster")){
mat <- terra::as.matrix(rast, wide = TRUE)
}else if (inherits(rast,"matrix")){
mat <- rast
}else{
stop("rast parameter must be on of matrix or SpatRaster")
}
if(inherits(window,"matrix")){
fun <- moranI_matrix_window
}else if (inherits(window, "numeric")){
#fun <- moranI_matrix
fun <- moranI_matrix_window
size <- n+1+n
window <- matrix(1, nrow = size, ncol = size)
}else{
stop("window parameter must be an integer or a matrix")
}
val <- fun(mat, window)
return(val)
}
#' @title Local Moran I for raster
#'
#' @description Calculate the Local Moran I for a numeric raster
#'
#' @param rast A SpatRaster or a matrix
#' @param window The window defining the neighbour weights
#' @return A SpatRaster or a matrix depending on the input with the local Moran I values
#' @examples
#' Arcachon <- terra::rast(system.file("extdata/Littoral4_2154.tif", package = "geocmeans"))
#' names(Arcachon) <- c("blue", "green", "red", "infrared", "SWIR1", "SWIR2")
#' rast <- Arcachon[[1]]
#' w <- matrix(1, nrow = 3, ncol = 3)
#' calc_local_moran_raster(rast, w)
calc_local_moran_raster <- function(rast, window){
if(inherits(rast,"SpatRaster")){
mat <- terra::as.matrix(rast, wide = TRUE)
}else if (inherits(rast,"matrix")){
mat <- rast
}else{
stop("rast parameter must be on of matrix or SpatRaster")
}
if(inherits(window,"matrix")){
window <- window
}else if (inherits(window,"numeric")){
w <- 1+2*inherits
window <- matrix(1, nrow = w, ncol = w)
window[ceiling(w/2),ceiling(w/2)] <- 0
}else{
stop("window parameter must be an integer or a matrix")
}
vals <- local_moranI_matrix_window(mat, window)
if(inherits(rast,"SpatRaster")){
terra::values(rast) <- vals
return(rast)
}else{
return(matrix(vals, ncol = ncol(mat), nrow = nrow(mat)))
}
}