agreement between workflows, areaproportion according to maxAgreement…

… levels and LCZ mode.
MGousseff · Oct 1, 2024 · f15831e · f15831e
1 parent 330c8c8
commit f15831e
Show file tree

Hide file tree

Showing 3 changed files with 77 additions and 37 deletions.
diff --git a/R/LCZmodeByAgreementLevel.R b/R/LCZmodeByAgreementLevel.R
@@ -0,0 +1,57 @@
+LCZmodeByAgreementLevel <- 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))
+
+  LCZwfsNames<-grep( pattern = "LCZ*",  x = names(intersec_sf), value = TRUE)
+  intersec_sf$LCZmode<-apply(intersec_sf[,LCZwfsNames], 1, Mode)
+
+  modeLCZSurfbyAgreement <- intersec_sf %>% group_by(maxAgree, LCZmode) %>% summarize(modeLCZsurf = sum(area)) %>% mutate(modeLCZSurfPerc = modeLCZsurf/sum(modeLCZsurf)*100)
+
+  generalProp<-intersec_sf %>%select(area, LCZmode) %>% mutate(totalArea=sum(area)) %>% 
+    group_by(LCZmode) %>% 
+    summarize(modeLCZGenSurfPerc = sum(area), totalArea = mean(totalArea)) %>% 
+    mutate(modeLCZGenSurfPerc = modeLCZGenSurfPerc / totalArea *100 ) %>% 
+    select(LCZmode, modeLCZGenSurfPerc)
+
+  modeLCZSurfbyAgreement<-left_join(modeLCZSurfbyAgreement, generalProp, by = "LCZmode") %>% 
+    arrange(desc(maxAgree),desc(modeLCZSurfPerc))
+
+  # 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(modeLCZSurfbyAgreement )
+}
+
+Mode <- function(x) {
+  ux <- unique(x)
+  unlist(ux[which.max(tabulate(match(x, ux)))])
+}
+
+
+LCZmodeTest <- LCZmodeByAgreementLevel(multicompare_test$intersec_sf)
+LCZmodeTest[601:610,c(LCZwfsNamesTest,"LCZmode")]
+
+test1<-multicompare_test$intersec_sf
+LCZwfsNamesTest<-grep( pattern = "LCZ*",  x = names(test1), value = TRUE) 
+apply(test1[,LCZwfsNamesTest], 1, Mode)
+
+
+
+grep( pattern = "LCZ*", names(multicompare_test$intersec_sf), value = TRUE)
diff --git a/R/compareMultipleLCZ.R b/R/compareMultipleLCZ.R
@@ -36,7 +36,7 @@ compareMultipleLCZ<-function(sfList, LCZcolumns, refCrs=NULL, sfWf=NULL, trimPer
     X = intersec_sfnogeom[,1:length(sfList)], MARGIN = 1, function(x) max(table(x) ))
   intersec_sf<-cbind(intersec_sfnogeom,intersec_sf$geometry)  %>% st_as_sf()
   intersec_sf
-  intersec_sfLong<-pivot_longer(st_drop_geometry(intersec_sf),cols=rangeCol, names_to = "whichWfs", values_to = "agree")
+  intersec_sfLong<-pivot_longer(intersec_sfnogeom,cols=rangeCol, names_to = "whichWfs", values_to = "agree")
   intersec_sfLong$LCZref<-substr(intersec_sfLong$whichWfs,start = 1, stop=1 )
   print(head(intersec_sfLong[,c(1,2,9:10)]))
   whichLCZagree <- names(intersec_sfLong)[as.numeric(intersec_sfLong$LCZref)]
@@ -49,48 +49,30 @@ 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"))
-
-# # 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_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)
-
-# 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))
 

diff --git a/R/buildWorkflowAgreement.R → R/workflowsAgreement.R b/R/buildWorkflowAgreement.R → R/workflowsAgreement.R
@@ -1,4 +1,4 @@
-buildWorkflowAgreement <- function(intersec_sf, sfWfs = NULL){ 
+workflowsAgreement <- function(intersec_sf, sfWfs = NULL){ 
   if ( !is.null(intersec_sf$geometry)) {
     intersec_sf<-st_drop_geometry(intersec_sf)
   }
@@ -26,5 +26,6 @@ buildWorkflowAgreement <- function(intersec_sf, sfWfs = NULL){
  return(sort(agreement_by_pair[,1], decreasing = TRUE))
 }
 
-tetest<-buildWorkflowAgreement(intersec_sf = multicompare_test$intersec_sf, 
+
+tetest<-workflowsAgreement(intersec_sf = multicompare_test$intersec_sf, 
   sfWfs = c("BDT11", "BDT22", "OSM11", "OSM22", "WUDAPT"))