buildWorkFlowAgreementOK. TODO : build LCZ agreement and then all sam…

…ple indicators
MGousseff · Sep 30, 2024 · 330c8c8 · 330c8c8
1 parent 2d6a4c5
commit 330c8c8
Show file tree

Hide file tree

Showing 10 changed files with 194 additions and 132 deletions.
diff --git a/NAMESPACE b/NAMESPACE
@@ -4,6 +4,7 @@ export(CohenKappa)
 export(LCZareas)
 export(areColors)
 export(compareLCZ)
+export(compareMultipleLCZ)
 export(confidSensib)
 export(fetchLCZ)
 export(groupLCZ)

diff --git a/R/buildWorkflowAgreement.R b/R/buildWorkflowAgreement.R
@@ -0,0 +1,30 @@
+buildWorkflowAgreement <- function(intersec_sf, sfWfs = NULL){ 
+  if ( !is.null(intersec_sf$geometry)) {
+    intersec_sf<-st_drop_geometry(intersec_sf)
+  }
+  columnNames<-names(intersec_sf)
+  pairNames<-grep(pattern = "[1-9]_[1-9]", x = columnNames)
+  agreement_by_pair<- t(
+    intersec_sf[,pairNames]) %*% as.matrix(drop_units(intersec_sf$area)) / 
+    sum(drop_units(intersec_sf$area))
+
+
+  if (!is.null(sfWfs)) {
+    lengthSfWfs<-length(sfWfs)
+    testLengthSfWfs<-factorial(lengthSfWfs)/(2*factorial(lengthSfWfs-2))
+    compNames<-NULL
+    if ( nrow(agreement_by_pair)==testLengthSfWfs ) {
+      for (firstWfIndice in 1:(length(sfWfs)-1)) {
+        for(secondWfIndice in (firstWfIndice + 1):length(sfWfs)){
+          compNames<-c(compNames,paste0(sfWfs[firstWfIndice],"_",sfWfs[secondWfIndice]))    
+          }
+        }
+      }
+    row.names(agreement_by_pair) <- compNames
+
+     }
+ return(sort(agreement_by_pair[,1], decreasing = TRUE))
+}
+
+tetest<-buildWorkflowAgreement(intersec_sf = multicompare_test$intersec_sf, 
+  sfWfs = c("BDT11", "BDT22", "OSM11", "OSM22", "WUDAPT"))
diff --git a/R/compareLCZ.R b/R/compareLCZ.R
@@ -314,19 +314,19 @@ compareLCZ<-function(sf1, geomID1="", column1 = "LCZ_PRIMARY", confid1="", wf1="
    # # Intersect geometries of both files
   ######################################################
     #intersection of geometries
-  echInt<-st_intersection(x=sf1[,nom1],y=sf2[,nom2]) %>% st_buffer(0) %>% 
+  intersec_sf<-st_intersection(x=sf1[,nom1],y=sf2[,nom2]) %>% st_buffer(0) %>% 
     mutate(area=drop_units(st_area(geometry)))
-  nbNoSurf<-nrow(subset(echInt,area==0))
+  nbNoSurf<-nrow(subset(intersec_sf,area==0))
   if (nbNoSurf>0){
   message(paste0(
     "The intersection of the two data set geometries return ",
     nbNoSurf, " geometries with an null area. They will be discarded."))
-  echInt<-subset(echInt,area!=0)
+  intersec_sf<-subset(intersec_sf,area!=0)
   }
 
 
     # checks if the two LCZ classifications agree
-  echInt$agree<-subset(echInt,select=column1,drop=T)==subset(echInt,select=column2,drop=T)
+  intersec_sf$agree<-subset(intersec_sf,select=column1,drop=T)==subset(intersec_sf,select=column2,drop=T)
 
 
 ######################################################
@@ -337,11 +337,11 @@ compareLCZ<-function(sf1, geomID1="", column1 = "LCZ_PRIMARY", confid1="", wf1="
 
   # Export of lcz and area for each geom for further analysis
 
-        #echInt<-echInt %>% mutate(area=st_area(geometry)) %>% drop_units
+        #intersec_sf<-intersec_sf %>% mutate(area=st_area(geometry)) %>% drop_units
     # Drop intersected geometries with area equal to zero
-        echInt<-subset(echInt,area!=0)
+        intersec_sf<-subset(intersec_sf,area!=0)
 
-        echIntExpo<-echInt %>% mutate(location=location,area=as.numeric(area)) %>%
+        intersec_sfExpo<-intersec_sf %>% mutate(location=location,area=as.numeric(area)) %>%
           st_set_geometry(NULL) %>% as.data.frame()
 
 
@@ -357,12 +357,12 @@ compareLCZ<-function(sf1, geomID1="", column1 = "LCZ_PRIMARY", confid1="", wf1="
                        getwd())
                 )
         if (!file.exists(filePath)){
-        write.table(x=echIntExpo, file =nom, append = TRUE, quote = TRUE, sep = ";",
+        write.table(x=intersec_sfExpo, file =nom, append = TRUE, quote = TRUE, sep = ";",
                     eol = "\n", na = "NA", dec = ".",
                    qmethod = c("escape", "double"),
                     fileEncoding = "", col.names=TRUE,row.names=F)
         }else{
-          write.table(x=echIntExpo, file =nom, append = TRUE, quote = TRUE, sep = ";",
+          write.table(x=intersec_sfExpo, file =nom, append = TRUE, quote = TRUE, sep = ";",
                       eol = "\n", na = "NA", dec = ".",
                       qmethod = c("escape", "double"),
                       fileEncoding = "", col.names=FALSE,row.names=F)
@@ -374,7 +374,7 @@ compareLCZ<-function(sf1, geomID1="", column1 = "LCZ_PRIMARY", confid1="", wf1="
 
 matConfOut<-matConfLCZ(sf1=sf1, column1=column1, sf2=sf2, column2=column2,
                        repr=repr, typeLevels=LCZlevels, plot=FALSE)
-matConfOut$data<-echIntExpo
+matConfOut$data<-intersec_sfExpo
 matConfLong<-matConfOut$matConf
 matConfLarge<-pivot_wider(matConfLong,names_from = column2,values_from = agree)
 matConfLarge<-matConfLarge %>% as.data.frame()
@@ -440,7 +440,7 @@ if (plot == TRUE){
 
   nbgeom1<-nrow(sf1)
   nbgeom2<-nrow(sf2)
-  nbgeomInter<-nrow(echInt)
+  nbgeomInter<-nrow(intersec_sf)
 
   # Plot the first classification
      l1Plot<- ggplot2::ggplot(boundary)+ # les données
@@ -463,7 +463,7 @@ if (plot == TRUE){
   # Plot areas where classifications agree
      agreePlot<-ggplot(boundary)+
        geom_sf(data=boundary, fill=NA,lty='blank')+
-       geom_sf(data=echInt,aes(fill=agree),lwd=0,colour=NA)+
+       geom_sf(data=intersec_sf,aes(fill=agree),lwd=0,colour=NA)+
        scale_fill_manual(values=c("red","green"),
                          name=paste0(
                            "The two classifications agree for \n ",percAgg, " % of the area Agreement"))+

diff --git a/R/compareMultipleLCZ.R b/R/compareMultipleLCZ.R
@@ -10,11 +10,7 @@
 #' @importFrom ggplot2 geom_sf guides ggtitle aes
 #' @import sf dplyr cowplot forcats units tidyr RColorBrewer utils grDevices rlang
 #' @return returns graphics of comparison and an object called matConfOut which contains :
-#' matConfLong, a confusion matrix in a longer form, 
-#' matConfPlot is a ggplot2 object showing the confusion matrix.
-#' percAgg is the general agreement between the two sets of LCZ, expressed as a percentage of the total area of the study zone
-#' pseudoK is a heuristic estimate of a Cohen's kappa coefficient of agreement between classifications
-#' If saveG is not an empty string, graphics are saved under "saveG.png"
+#'TO DO
 #' @export
 #' @examples
 #' 
@@ -53,88 +49,88 @@ compareMultipleLCZ<-function(sfList, LCZcolumns, refCrs=NULL, sfWf=NULL, trimPer
 }
 
 
-sfBDT_11_78030<-importLCZvect(dirPath="/home/gousseff/Documents/3_data/data_article_LCZ_diff_algos/GeoClimate/2011/bdtopo_2_78030/",
-                   file="rsu_lcz.fgb", column="LCZ_PRIMARY")
-class(sfBDT_11_78030)
-sfBDT_22_78030<-importLCZvect(dirPath="/home/gousseff/Documents/3_data/data_article_LCZ_diff_algos/GeoClimate/2022/bdtopo_3_78030/",
-                              file="rsu_lcz.fgb", column="LCZ_PRIMARY")
-sf_OSM_11_Auffargis<-importLCZvect(dirPath="//home/gousseff/Documents/3_data/data_article_LCZ_diff_algos/GeoClimate/2011/osm_Auffargis/",
-                                   file="rsu_lcz.fgb", column="LCZ_PRIMARY")
-sf_OSM_22_Auffargis<-importLCZvect(dirPath="/home/gousseff/Documents/3_data/data_article_LCZ_diff_algos/GeoClimate/2022/osm_Auffargis/",
-                                   file="rsu_lcz.fgb", column="LCZ_PRIMARY")
-sf_WUDAPT_78030<-importLCZvect("/home/gousseff/Documents/3_data/data_article_LCZ_diff_algos/WUDAPT/",
-                               file ="wudapt_Auffargis.fgb", column="lcz_primary")
+# sfBDT_11_78030<-importLCZvect(dirPath="/home/gousseff/Documents/3_data/data_article_LCZ_diff_algos/GeoClimate/2011/bdtopo_2_78030/",
+#                    file="rsu_lcz.fgb", column="LCZ_PRIMARY")
+# class(sfBDT_11_78030)
+# sfBDT_22_78030<-importLCZvect(dirPath="/home/gousseff/Documents/3_data/data_article_LCZ_diff_algos/GeoClimate/2022/bdtopo_3_78030/",
+#                               file="rsu_lcz.fgb", column="LCZ_PRIMARY")
+# sf_OSM_11_Auffargis<-importLCZvect(dirPath="//home/gousseff/Documents/3_data/data_article_LCZ_diff_algos/GeoClimate/2011/osm_Auffargis/",
+#                                    file="rsu_lcz.fgb", column="LCZ_PRIMARY")
+# sf_OSM_22_Auffargis<-importLCZvect(dirPath="/home/gousseff/Documents/3_data/data_article_LCZ_diff_algos/GeoClimate/2022/osm_Auffargis/",
+#                                    file="rsu_lcz.fgb", column="LCZ_PRIMARY")
+# sf_WUDAPT_78030<-importLCZvect("/home/gousseff/Documents/3_data/data_article_LCZ_diff_algos/WUDAPT/",
+#                                file ="wudapt_Auffargis.fgb", column="lcz_primary")
 
-sfList<-list(BDT11 = sfBDT_11_78030, BDT22 = sfBDT_22_78030, OSM11= sf_OSM_11_Auffargis, OSM22 = sf_OSM_22_Auffargis,
-             WUDAPT = sf_WUDAPT_78030)
+# sfList<-list(BDT11 = sfBDT_11_78030, BDT22 = sfBDT_22_78030, OSM11= sf_OSM_11_Auffargis, OSM22 = sf_OSM_22_Auffargis,
+#              WUDAPT = sf_WUDAPT_78030)
 
-intersected<-createIntersec(sfList = sfList, LCZcolumns = c(rep("LCZ_PRIMARY",4),"lcz_primary"),
-                            sfWf = c("BDT11","BDT22","OSM11","OSM22","WUDAPT"))
+# intersected<-createIntersec(sfList = sfList, LCZcolumns = c(rep("LCZ_PRIMARY",4),"lcz_primary"),
+#                             sfWf = c("BDT11","BDT22","OSM11","OSM22","WUDAPT"))
 
-# test_list<-list(a=c(1,2),b="top",c=TRUE)
-# length(test_list)
-# for (i in test_list[2:3]) print(str(i))
+# # test_list<-list(a=c(1,2),b="top",c=TRUE)
+# # length(test_list)
+# # for (i in test_list[2:3]) print(str(i))
 
 multicompare_test<-compareMultipleLCZ(sfList = sfList, LCZcolumns = c(rep("LCZ_PRIMARY",4),"lcz_primary"),
                                       sfWf = c("BDT11","BDT22","OSM11","OSM22","WUDAPT"),trimPerc = 0.5)
 multicompare_test
 
-test<-multicompare_test$intersec_sfLong
-test2<-test %>% subset(agree==TRUE) %>% group_by(LCZvalue) %>% summarize(agreementArea=sum(area)) %>% mutate(percAgreementArea=agreementArea/sum(agreementArea))
+# test<-multicompare_test$intersec_sfLong
+# test2<-test %>% subset(agree==TRUE) %>% group_by(LCZvalue) %>% summarize(agreementArea=sum(area)) %>% mutate(percAgreementArea=agreementArea/sum(agreementArea))
 
-test<-multicompare_test$intersec_sf[,1:5] %>% st_drop_geometry()
-prov1<-apply(X = test, MARGIN = 1, table )
-prov2<-apply(X = test, MARGIN = 1, function(x) max(table(x)) )
+# test<-multicompare_test$intersec_sf[,1:5] %>% st_drop_geometry()
+# prov1<-apply(X = test, MARGIN = 1, table )
+# prov2<-apply(X = test, MARGIN = 1, function(x) max(table(x)) )
 
-head(prov1)
-head(prov2)
+# head(prov1)
+# head(prov2)
 
-plot1<-showLCZ(sf = multicompare_test$intersec_sf, column="LCZBDT22", wf="22")
-plot2<-showLCZ(sf = multicompare_test$intersec_sf, column="LCZBDT11", wf="11")
+# plot1<-showLCZ(sf = multicompare_test$intersec_sf, column="LCZBDT22", wf="22")
+# plot2<-showLCZ(sf = multicompare_test$intersec_sf, column="LCZBDT11", wf="11")
 
-ggplot(data=multicompare_test$intersec_sf) +
-  geom_sf(aes(fill=maxAgree, color=after_scale(fill)))+
-  scale_fill_gradient(low = "red" , high = "green", na.value = NA)
+# ggplot(data=multicompare_test$intersec_sf) +
+#   geom_sf(aes(fill=maxAgree, color=after_scale(fill)))+
+#   scale_fill_gradient(low = "red" , high = "green", na.value = NA)
 
-hist(st_area(multicompare_test$intersec_sf$geometry))
+# hist(st_area(multicompare_test$intersec_sf$geometry))
 
-allSameList<-list(OSM11= sf_OSM_11_Auffargis, OSM11.2 = sf_OSM_11_Auffargis, 
-  OSM11.3 = sf_OSM_11_Auffargis, OSM11.4 = sf_OSM_11_Auffargis, OSM_11.5 = sf_OSM_11_Auffargis)
-showLCZ(sfList[[1]])
+# allSameList<-list(OSM11= sf_OSM_11_Auffargis, OSM11.2 = sf_OSM_11_Auffargis, 
+#   OSM11.3 = sf_OSM_11_Auffargis, OSM11.4 = sf_OSM_11_Auffargis, OSM_11.5 = sf_OSM_11_Auffargis)
+# showLCZ(sfList[[1]])
 
-sf_OSM_11_Auffargis[which.max(st_area(sf_OSM_11_Auffargis)),] 
+# sf_OSM_11_Auffargis[which.max(st_area(sf_OSM_11_Auffargis)),] 
 
-max(st_area(sf_OSM_11_Auffargis))
+# max(st_area(sf_OSM_11_Auffargis))
 
-multicompare_test_all_same<-compareMultipleLCZ(sfList = allSameList, 
-  LCZcolumns = rep("LCZ_PRIMARY",5),
-  sfWf = c("OSM1","OSM2", "OSM3", "OSM4", "OSM5"),
-  trimPerc = 0.5)
+# multicompare_test_all_same<-compareMultipleLCZ(sfList = allSameList, 
+#   LCZcolumns = rep("LCZ_PRIMARY",5),
+#   sfWf = c("OSM1","OSM2", "OSM3", "OSM4", "OSM5"),
+#   trimPerc = 0.5)
 
 
 
-areas_test<-st_area(multicompare_test_all_same$intersec_sf)
-hist(areas_test)
-hist(st_area(sf_OSM_11_Auffargis$geometry))
+# areas_test<-st_area(multicompare_test_all_same$intersec_sf)
+# hist(areas_test)
+# hist(st_area(sf_OSM_11_Auffargis$geometry))
 
 
-quantile(areas_test, prob = 0.5)
-test2<-multicompare_test_all_same$intersec_sf[
-  st_area(multicompare_test_all_same$intersec_sf) == 
-    max(st_area(multicompare_test_all_same$intersec_sf)),
-]
+# quantile(areas_test, prob = 0.5)
+# test2<-multicompare_test_all_same$intersec_sf[
+#   st_area(multicompare_test_all_same$intersec_sf) == 
+#     max(st_area(multicompare_test_all_same$intersec_sf)),
+# ]
 
 
-ggplot() + 
-    geom_sf(data=multicompare_test_all_same$intersec_sf, aes(fill=maxAgree))+
-    scale_fill_gradient(low = "red" , high = "green", na.value = NA)
+# ggplot() + 
+#     geom_sf(data=multicompare_test_all_same$intersec_sf, aes(fill=maxAgree))+
+#     scale_fill_gradient(low = "red" , high = "green", na.value = NA)
 
-ggplot() + 
-    geom_sf(data=test2, aes(color = "gray", fill=maxAgree)) +
-      scale_fill_gradient(low = "red" , high = "green", na.value = NA) + 
-    scale_linewidth(range=c(8))
+# ggplot() + 
+#     geom_sf(data=test2, aes(color = "gray", fill=maxAgree)) +
+#       scale_fill_gradient(low = "red" , high = "green", na.value = NA) + 
+#     scale_linewidth(range=c(8))
 
-ggplot() +
-  geom_sf(data = sf_OSM_11_Auffargis[which.max(st_area(sf_OSM_11_Auffargis)),], 
-aes(color = LCZ_PRIMARY) 
-)
+# ggplot() +
+#   geom_sf(data = sf_OSM_11_Auffargis[which.max(st_area(sf_OSM_11_Auffargis)),], 
+# aes(color = LCZ_PRIMARY) 
+# )