Skip to content
Snippets Groups Projects
functionssML.R 8.16 KiB
Newer Older
# recode CLC -------------------
Nicolas Saby's avatar
Nicolas Saby committed
occuprecode = function(attr) {
  recode(attr, 
         "10" = "0" , "11"= "0" , 
         "12"= "1" , "14"= "1" ,
         "15" = "2" ,"16"= "2" , "17"= "2" , 
         "18" = "3" ,"19"="3",  "20" = "3" ,
         "21" = "3" ,"22"  = "3","
                          23" = "4" ,  "24"  = "4" ,"25" = "4" ,
         "26" = "3" , "27"  = "3" ,"28" = "3" ,
         "19" = "3" , "36" = "3" ,
Nicolas Saby's avatar
Nicolas Saby committed
         "32" = "3",
         .default = "A")
}
# recode(sample2$FR_CLC, 
#        "10" = "0" , "11"= "0" , 
#        "12"= "1" , "14"= "1" ,
#        "15" = "2" ,"16"= "2" , "17"= "2" , 
#        "18" = "3" ,"19"="3",  "20" = "3" ,
#        "21" = "3" ,"22"  = "3","
#                           23" = "4" ,  "24"  = "4" ,"25" = "4" ,
#        "26" = "3" , "27"  = "3" ,"28" = "3" ,
#        "32" = "3",
#        .default = "A")
Nicolas Saby's avatar
Nicolas Saby committed


# optimisation of the SMN --------------------
Nicolas Saby's avatar
Nicolas Saby committed


Nicolas Saby's avatar
Nicolas Saby committed
optimisesampleSML <-function(
Nicolas Saby's avatar
Nicolas Saby committed
    framecens = NULL,
    X,
    Y,
    domainvalue,
    cv,
    eval= FALSE # should we run the ex ante evaluation
Nicolas Saby's avatar
Nicolas Saby committed
  ) {
  # in case of presence of census data  to be kept
Nicolas Saby's avatar
Nicolas Saby committed
  if ( ! is.null(framecens)) {
    
    print(" Include points from existing surveys (forcing to keep existing points) ")
    
    tab1 <- table(framesamp[,domainvalue])
    framesamp <- framesamp[framesamp[,domainvalue] %in% names(tab1[tab1 >= 10]),] 
    
    tab1 <- table(framecens[,domainvalue])
    framecens <- framecens[framecens[,domainvalue] %in% names(tab1[tab1 >= 10]),] 
    
    
    print(dim(framecens))
    
    sample2 <- rbind(framesamp, 
                     framecens)
    
    
Nicolas Saby's avatar
Nicolas Saby committed
    sample2$id <- c(1:nrow(sample2))
    

    frame1 <- buildFrameDF(df = sample2,
                           id = "id",
                           X = X,
                           Y = Y,
                           domainvalue = domainvalue)
    
    #Build initial strata
    strata1 <- buildStrataDF(frame1, progress=T)
    
    
    
    #Find initial solutions for each domain (NUTS2) by kmeans clustering
    
    init_sol3 <- KmeansSolution2(frame=frame1,
                                 errors=cv,
                                 maxclusters = 8, # why 8?
                                 showPlot = F) 
    
    nstrata3 <- tapply(init_sol3$suggestions,
                       init_sol3$domainvalue,
                       FUN=function(x) length(unique(x)))
    
    initial_solution3 <- prepareSuggestion(init_sol3,frame1,nstrata3)
    
    
    #Optimize initial solution by Genetic Algorithm
    print("optimStrata-------------------")
    
    ind_framecens <- c(rep(T, dim(framesamp)[1]),
                       rep(F, dim(framecens)[1]))

    frame1samp <- frame1[ ind_framecens , ]
    frame1cens <- frame1[ !ind_framecens , ]
    
    
    print(dim(frame1cens))
    
    solution3 <- optimStrata(method = "continuous",
                             errors = cv, 
                             framesamp = frame1samp ,
                             framecens =  frame1cens ,
                             iter = 20,
                             pops = 20,
                             nStrata = nstrata3,
                             suggestions = initial_solution3,
                             showPlot=FALSE,writeFiles=FALSE)
    
    framenew3 <- solution3$framenew
    outstrata3 <- solution3$aggr_strata
    
    strataStructure3 <- summaryStrata(framenew3,
                                      outstrata3,
                                      progress=T)
    
    #Total number of samples required
    
    
    print("selectSample-------------------")
    
    sample4 <- selectSample(framenew3, 
                            outstrata3,
                            writeFiles = F)
    
    
  } else {
    
    ## case when framecens is null
Nicolas Saby's avatar
Nicolas Saby committed
    sample2<- framesamp
    
    
    frame1 <- buildFrameDF(df = sample2,
                           id = "id",
                           X = X,
                           Y = Y,
                           domainvalue = domainvalue)
    
    #Build initial strata
    strata1 <- buildStrataDF(frame1, progress=T)
    
    
    
    #Find initial solutions for each domain (NUTS2) by kmeans clustering
    
    init_sol3 <- KmeansSolution2(frame=frame1,
                                 errors=cv,
                                 maxclusters = 8, # why 8?
                                 showPlot = F) 
    
    nstrata3 <- tapply(init_sol3$suggestions,
                       init_sol3$domainvalue,
                       FUN=function(x) length(unique(x)))
    
    initial_solution3 <- prepareSuggestion(init_sol3,frame1,nstrata3)
    
    
    #Optimize initial solution by Genetic Algorithm
    print("optimStrata-------------------")
    
    solution3 <- optimStrata(method = "continuous",
                             errors = cv, 
                             framesamp = frame1,
                             framecens = framecens ,
                             iter = 20,
                             pops = 20,
                             nStrata = nstrata3,
                             suggestions = initial_solution3,
                             showPlot=FALSE,writeFiles=FALSE)
    
    framenew3 <- solution3$framenew
    outstrata3 <- solution3$aggr_strata
    
    strataStructure3 <- summaryStrata(framenew3,
                                      outstrata3,
                                      progress=T)
    
    #Total number of samples required
    
    
    print("selectSample-------------------")
    
    sample4 <- selectSample(framenew3, 
                            outstrata3,
                            writeFiles = F)
    
  }# End else
    
Nicolas Saby's avatar
Nicolas Saby committed

  
  
print("evalSolution-------------------")
  
  res = eval3 = NA
  
  if(eval == TRUE)  {
    eval3 <- evalSolution(frame = solution3$framenew, 
                          outstrata = solution3$aggr_strata, 
                          nsampl = 200, # why not the optimised number of sample ?
                          progress = TRUE,
                          writeFiles = FALSE) 
    res = expected_CV(solution3$aggr_strata)
  } 
   
Nicolas Saby's avatar
Nicolas Saby committed
   return(list(
     smlsample = sample4,
    optimstrata=  solution3,
          tableau = framenew3,
          eval = eval3,
          NbSample = sum(strataStructure3$Allocation),
          ExpCV = res)
          )
  
}

## ----------- plot

plotResu <- function(
    res, # the output of the function optim , a list of results
    nom , # the name of the file to be saved
    crs,
    sample2 ,
    stack500 
    ) {
Nicolas Saby's avatar
Nicolas Saby committed
  
  points4 <- sample2[res$smlsample$ID,]
  
  # coordinates(points4)=~x+y
  points4 <- st_as_sf(points4,
                      coords=c("x","y"),
                      crs = crs
  )
  
  
  # plotStrata2d(res$optimstrata$framenew, 
  #              res$optimstrata$aggr_strata,
  #              domain = 5, 
  #              vars = c("X1","X2"),
  #              labels = c("OC","pH"))
  # 
  
  # dom = 1
  # hist(eval3$est$Y1[res$eval3$est$dom == dom], col = "grey", border = "white",
  #      xlab = expression(hat(Y)[1]),
  #      freq = FALSE, 
  #      main = paste("Variable Y1 Domain ",dom,sep=""))
  # abline(v = mean(eval3$est$Y1[eval3$est$dom == dom]), col = "blue", lwd = 2)
  # abline(v = mean(frame1$Y1[frame1$domainvalue==dom]), col = "red")
  # legend("topright", c("distribution mean", "true value"),
  #        lty = 1, col = c("blue", "red"), box.col = NA, cex = 0.8)
  
  

  
Nicolas Saby's avatar
Nicolas Saby committed
  plNUTS = tm_shape(stack500["FR_pH"]) + 
Nicolas Saby's avatar
Nicolas Saby committed
    tm_raster(style = "quantile",n = 10,
              palette = "RdYlGn",
              alpha = 0.5) +
    tm_shape(points4) +
    tm_dots(col = "black",
            size = 0.002)+
    tm_layout(legend.outside = T)
  
Nicolas Saby's avatar
Nicolas Saby committed
  tmap_save(plNUTS,
            filename = paste0("output/",nom,".jpeg" ) 
            
            
  )
  
Nicolas Saby's avatar
Nicolas Saby committed
  # points4 <- st_as_sf(sample2,
  #                     coords=c("x","y"),
  #                     crs = crs
  # )
  # 
  # plNUTS = tm_shape(stack500["FR_pH"]) + 
  #   tm_raster(style = "quantile",n = 10,
  #             palette = "RdYlGn",
  #             alpha = 0.5) +
  #   tm_shape(points4) +
  #   tm_dots(col = "black",
  #           size = 0.002)+
  #   tm_layout(legend.outside = T)
  # 
  # tmap_save(plNUTS,
  #           filename = paste0("output/",nom,"Start.jpeg" ) 
  #           
  #           
  # )
  # 
Nicolas Saby's avatar
Nicolas Saby committed
  
Nicolas Saby's avatar
Nicolas Saby committed
}