From 6e366f54707ee72ebbca501735a870a32d63bede Mon Sep 17 00:00:00 2001
From: Francesco Sabatini <francesco.sabatini@idiv.de>
Date: Mon, 30 Nov 2020 18:54:36 +0100
Subject: [PATCH] Fixed bug in A98 optimized

---
 code/A98_PredictorsExtract.R | 20 +++++++++++---------
 1 file changed, 11 insertions(+), 9 deletions(-)

diff --git a/code/A98_PredictorsExtract.R b/code/A98_PredictorsExtract.R
index 19e18c0..b630aac 100644
--- a/code/A98_PredictorsExtract.R
+++ b/code/A98_PredictorsExtract.R
@@ -25,11 +25,8 @@ PredExtr <- function(x.shp, myfunction=NA, output=NA,
   print(paste("output will be:", output))
   if(typp=="raster"){ mypredictor <- raster(toextract)} else
   mypredictor <- readOGR(toextract)
-  colnames(mypredictor@data)
   header.shp <- readOGR(x.shp)
   crs(mypredictor) <- crs(header.shp) #should be verified beforehand!
-  
-  
 
   ## Divide in chunks if requested
   if(chunkn>1 & !is.na(chunk.i)){
@@ -47,14 +44,16 @@ PredExtr <- function(x.shp, myfunction=NA, output=NA,
   cl <- makeForkCluster(ncores, outfile="" )
   registerDoParallel(cl)
 
-  print("start main foreach loop")
+  
   if(typp=="raster"){
-  out <- foreach(i=1:length(header.shp), .packages=c('raster'), .combine=rbind) %dopar% { 
+    print("start main foreach loop")
+    out <- foreach(i=1:length(header.shp), .packages=c('raster'), .combine=rbind) %dopar% { 
     tmp <- raster::extract(mypredictor, header.shp[i,], 
                          buffer=min(10000,  max(250, header.shp@data$lc_ncrt[i])), fun=myfunction) 
   }
   if(!is.na(output)) {write.csv(out, file = output)} else{return(out)}
   } else {
+    print("Match sequentially")
     out <- sp::over(x=header.shp, y=mypredictor) 
     i.toassign <- which(is.na(out[,1]))
     toassign <- header.shp[i.toassign,]
@@ -67,7 +66,7 @@ PredExtr <- function(x.shp, myfunction=NA, output=NA,
         as_Spatial()
       } else {
         mypredictor.expl <- mypredictor}
-    
+    print("start main foreach loop")
     nearest <- foreach(i=1:length(toassign), .packages=c('raster'), .combine=rbind) %dopar% { 
       print(i)
       ## create a subset of geoentities based on a 5° buffer radius around each target plot.
@@ -84,13 +83,16 @@ PredExtr <- function(x.shp, myfunction=NA, output=NA,
         }
       )
       # find nearest neighbour  
-      nearest.tmp <- tryCatch(geosphere::dist2Line(toassign[i,], tmp.mypredictor),
-                              error = function(e){data.frame(i=i, distance=NA, lon=NA, lat=NA,ID=NA)}
+      nearest.tmp <- tryCatch(tmp.mypredictor@data[geosphere::dist2Line(toassign[i,], tmp.mypredictor)[,"ID"],],
+                              error = function(e){
+                                ee <- myptedictor@data[1,]
+                                ee[1,] <- rep(NA, ncol(mypredictor))
+                                }
                               )
       return(nearest.tmp)
     }
 
-    out[i.toassign,] <- tmp.mypredictor@data[nearest[,"ID"],]
+    out[i.toassign,] <- nearest
    }
     if(!is.na(output)) {write.csv(out, file = output)} else{return(out)}
   }
-- 
GitLab