diff --git a/src/prediction/abundance_prediction.R b/src/prediction/abundance_prediction.R
index a2fa2204287677e912386f96dbb8bce9856e833d..a7c43fe9a89d37b1c9615ecc3de0cd54f5a0b71a 100644
--- a/src/prediction/abundance_prediction.R
+++ b/src/prediction/abundance_prediction.R
@@ -116,7 +116,7 @@ exclude_year <- as.numeric(options$exclude_year)
if(is_verbose){print(paste("exclude year", exclude_year))}
focal_change_predictor <- as.numeric(options$focal_change_predictor)
-if(is_verbose){print(paste("focal_change_predictor", focal_change_predictor))}
+if(is_verbose){print(paste("focal_change_predictor", focal_change_predictor))}
if(is_verbose){print(paste(Sys.time(), "0. start run"))}
#------------------------#
@@ -162,17 +162,9 @@ coeffs[is.na(coeffs) == T] <- 0
# Load estimates for observation and grid
# these are the predictors that will be used in the prediction
-# THEY MUST BE IN THE ORDER IN WHICH THEY APPEAR IN THE COEFFICIENTS
-# CODE THIS
# here only separate after first _
predictor_names_coeffs <- gsub("coeff_","", names(coeffs))
-#UNDERSTAND HERE WHAT IS HAPPENING
-#interaction_terms_names <- predictor_names_coeffs[predictor_names_coeffs %in%
-# paste0("year:", predictor_names_coeffs)]
-#interaction_terms_names <- gsub("year:", "", interaction_terms_names)
-#quadratic_terms_names <- predictor_names_coeffs[predictor_names_coeffs %in%
-# paste0("I(", predictor_names_coeffs, "^2)")]
-#quadratic_terms_names <- gsub("I\\(|\\^2\\)", "", quadratic_terms_names )
+
quadratic_terms_names <- c("rain_dry")
predictor_names_coeffs <- predictor_names_coeffs[predictor_names_coeffs != "(Intercept)"]
@@ -287,7 +279,6 @@ names(predictor_estimates) <- c("intercept", predictor_names,
if(is_verbose){print(paste("1. start pred_per_cell", Sys.time()))}
pred_per_cell <- foreach(i = 1:nrow(predictor_estimates), .combine = c) %dopar% {
- # pred_per_cell <- foreach(i = 1:100, .combine = c) %dopar% {
t_predictor_estimates <- t( predictor_estimates[i, ])
pred_estimates_wcoeffs <- data.frame(mapply(`*`, coeffs, t_predictor_estimates, SIMPLIFY = F))
pred_estimates_sum <- apply(pred_estimates_wcoeffs, 1, sum)
diff --git a/src/resource_use/resource_use_table_no_protected_country_split_MYS.py b/src/resource_use/resource_use_table_no_protected_country_split_MYS.py
index 4b99bd883111c8854ffa81149c9c78353224070a..0872ec669ed85bcac1c24df51789a4b3fbfe9a1a 100644
--- a/src/resource_use/resource_use_table_no_protected_country_split_MYS.py
+++ b/src/resource_use/resource_use_table_no_protected_country_split_MYS.py
@@ -622,19 +622,47 @@ abundance_data = mp.tiff.read_tif(abundance_path, 1)
resource_use = mp.tiff.read_tif(grid_layer_path, 1)
-abundance_data = np.where(absence == 0, abundance_data, 0)
-#abundance_data = np.where(abundance_data > 0.1, abundance_data, 0)
+abundance_data = np.where(populations != 0, abundance_data, 0)
+
# there are some pixels in the outside shape of borneo,
# where we have abundance but no grid, this we will clip
-abundance_data = np.where((resource_use == 0) &
- (abundance_data > 0), 0, abundance_data)
+
-abundance_data_MYS = np.where(borneo == 136, abundance_data, 0 )
+#abundance_data_MYS = np.where(borneo == 136, abundance_data, 0 )
abundance_data_SAB = np.where(borneo_province == 10, abundance_data, 0 )
abundance_data_SAW = np.where(borneo_province == 11, abundance_data, 0 )
abundance_data_IDN = np.where(borneo == 106, abundance_data, 0 )
+
+
+# checks
+"""
+mp.tiff.write_tif(file_with_srid = grid_layer_path,
+ full_output_name = "/homes/mv39zilo/work/Borneo/analysis/model_prep_and_running/results/resource_use/SAB.tif",
+ data = abundance_data_SAB ,
+ dtype = 4)
+
+mp.tiff.write_tif(file_with_srid = grid_layer_path,
+ full_output_name = "/homes/mv39zilo/work/Borneo/analysis/model_prep_and_running/results/resource_use/SAW.tif",
+ data = abundance_data_SAW ,
+ dtype = 4)
+
+mp.tiff.write_tif(file_with_srid = grid_layer_path,
+ full_output_name = "/homes/mv39zilo/work/Borneo/analysis/model_prep_and_running/results/resource_use/Kal.tif",
+ data = abundance_data_IDN ,
+ dtype = 4)
+
+"""
+np.sum(abundance_data) - (np.sum(abundance_data_SAB) + np.sum(abundance_data_SAW) + np.sum(abundance_data_IDN) )
+# 39 individuals in the whole abundance_data that are not in the part
+# so I need to get rid of those, otherwise slight inaccuracy
+
+abundance_data = abundance_data_IDN + abundance_data_SAB + abundance_data_SAW
+
+np.sum(abundance_data) - (np.sum(abundance_data_SAB) + np.sum(abundance_data_SAW) + np.sum(abundance_data_IDN) )
+
+
file_path_table_affected = "/homes/mv39zilo/work/Borneo/analysis/model_prep_and_running/results/resource_use/OU_area_affected_" + str(year) + ".csv"
affected = open(file_path_table_affected, 'w')
@@ -644,6 +672,10 @@ affected.write("type, country, ou_affected, perc_ou_affected, area_affected, per
# IT DOESNT MAKE A DIFFERENCE IF INCLUDING ONLY LARGER THAN 0.1
+# check different approaches to count > 0 pixels
+np.sum(total_ae_SAW)
+np.count_nonzero(abundance_data_SAW)
+
# how many OU have been affected
gone_ou_all = np.where(#(abundance_data > 0.1) &
@@ -653,8 +685,7 @@ gone_ou_all = np.where(#(abundance_data > 0.1) &
# # how much "area" was affected
total_ae = np.where(abundance_data > 0, 1, 0)
-gone_ae_all = np.where(
- # (abundance_data > 0.1) &
+gone_ae_all = np.where( (total_ae == 1) &
((resource_use == 1) |
(resource_use == 2) |
(resource_use == 3)), 1, 0)
@@ -665,13 +696,13 @@ str(np.sum(gone_ou_all)*100/np.sum(abundance_data)) + "," +
str(np.sum(gone_ae_all) * 0.0001) + "," + str(np.sum(gone_ae_all)*100/ np.sum(total_ae)) + "\n")
# how many OU have been affected
-gone_ou_defor = np.where(#(abundance_data > 0.1) &
+gone_ou_defor = np.where(
((resource_use == 1) |
(resource_use == 2) ), abundance_data, 0)
# # how much "area" was affected
gone_ae_defor = np.where(
- # (abundance_data > 0.1) &
+(total_ae == 1) &
((resource_use == 1) |
(resource_use == 2) ), 1, 0)
@@ -684,10 +715,8 @@ affected.write( "cover_change" + "," + "borneo" + "," +
-gone_ou_logging = np.where(#(abundance_data > 0.1) &
- (resource_use == 3), abundance_data, 0)
-gone_ae_logging = np.where(#(abundance_data > 0.1) &
- (resource_use == 3), 1, 0)
+gone_ou_logging = np.where((resource_use == 3), abundance_data, 0)
+gone_ae_logging = np.where((total_ae == 1) & (resource_use == 3), 1, 0)
affected.write( "logging" + "," + "borneo" + "," +
@@ -697,7 +726,7 @@ affected.write( "logging" + "," + "borneo" + "," +
str(np.sum(gone_ae_logging)*100/ np.sum(total_ae))+ "\n")
-gone_ou_IDN_all = np.where(#(abundance_data_IDN > 0.1) &
+gone_ou_IDN_all = np.where(
((resource_use_IDN == 1) |
(resource_use_IDN == 2) |
(resource_use_IDN == 3)), abundance_data_IDN, 0)
@@ -705,7 +734,7 @@ gone_ou_IDN_all = np.where(#(abundance_data_IDN > 0.1) &
# # how much "area" was affected
total_ae_IDN = np.where(abundance_data_IDN > 0, 1, 0)
gone_ae_IDN_all = np.where(
- # (abundance_data_IDN > 0.1) &
+ (total_ae_IDN == 1) &
((resource_use_IDN == 1) |
(resource_use_IDN == 2) |
(resource_use_IDN == 3)), 1, 0)
@@ -722,7 +751,7 @@ gone_ou_IDN_defor = np.where(#(abundance_data_IDN > 0.1) &
# # how much "area" was affected
gone_ae_IDN_defor = np.where(
- # (abundance_data_IDN > 0.1) &
+ (total_ae_IDN == 1) &
((resource_use_IDN == 1) |
(resource_use_IDN == 2) ), 1, 0)
@@ -735,9 +764,9 @@ affected.write( "cover_change" + "," + "IDN" + "," +
-gone_ou_IDN_logging = np.where(#(abundance_data_IDN > 0.1) &
+gone_ou_IDN_logging = np.where(
(resource_use_IDN == 3), abundance_data_IDN, 0)
-gone_ae_IDN_logging = np.where(#(abundance_data_IDN > 0.1) &
+gone_ae_IDN_logging = np.where((total_ae_IDN == 1) &
(resource_use_IDN == 3), 1, 0)
@@ -747,7 +776,7 @@ affected.write( "logging" + "," + "IDN" + "," +
str(np.sum(gone_ae_IDN_logging)* 0.0001) + "," +
str(np.sum(gone_ae_IDN_logging)*100/ np.sum(total_ae_IDN))+ "\n")
-
+"""
gone_ou_MYS_all = np.where(#(abundance_data_MYS > 0.1) &
((resource_use_MYS == 1) |
(resource_use_MYS == 2) |
@@ -799,20 +828,26 @@ affected.write( "logging" + "," + "MYS" + "," +
str(np.sum(gone_ae_MYS_logging) * 0.0001) + "," +
str(np.sum(gone_ae_MYS_logging)*100/ np.sum(total_ae_MYS)) + "\n")
+"""
+
# SAB
-gone_ou_SAB_all = np.where(#(abundance_data_SAB > 0.1) &
+gone_ou_SAB_all = np.where(
((resource_use_SAB == 1) |
(resource_use_SAB == 2) |
(resource_use_SAB == 3)), abundance_data_SAB, 0)
# # how much "area" was affected
total_ae_SAB = np.where(abundance_data_SAB > 0, 1, 0)
+
+
gone_ae_SAB_all = np.where(
- # (abundance_data_SAB > 0.1) &
+ (total_ae_SAB == 1) &
((resource_use_SAB == 1) |
(resource_use_SAB == 2) |
- (resource_use_SAB == 3)), 1, 0)
+ (resource_use_SAB == 3) ), 1, 0)
+
+
affected.write( "all" + "," + "SAB" + "," +
str(np.sum(gone_ou_SAB_all) ) + "," +
@@ -821,13 +856,13 @@ affected.write( "all" + "," + "SAB" + "," +
str(np.sum(gone_ae_SAB_all)*100/ np.sum(total_ae_SAB)) + "\n")
# how many OU have been affected
-gone_ou_SAB_defor = np.where(#(abundance_data_SAB > 0.1) &
+gone_ou_SAB_defor = np.where(
((resource_use_SAB == 1) |
(resource_use_SAB == 2)), abundance_data_SAB, 0)
# # how much "area" was affected
gone_ae_SAB_defor = np.where(
- # (abundance_data_SAB > 0.1) &
+ (total_ae_SAB == 1) &
((resource_use_SAB == 1) |
(resource_use_SAB == 2)), 1, 0)
@@ -840,9 +875,8 @@ affected.write( "cover_change" + "," + "SAB" + "," +
-gone_ou_SAB_logging = np.where(#(abundance_data_SAB > 0.1) &
- (resource_use_SAB == 3), abundance_data_SAB, 0)
-gone_ae_SAB_logging = np.where(#(abundance_data_SAB > 0.1) &
+gone_ou_SAB_logging = np.where( (resource_use_SAB == 3), abundance_data_SAB, 0)
+gone_ae_SAB_logging = np.where( (total_ae_SAB == 1) &
(resource_use_SAB == 3), 1, 0)
@@ -852,18 +886,24 @@ affected.write( "logging" + "," + "SAB" + "," +
str(np.sum(gone_ae_SAB_logging) * 0.0001) + "," +
str(np.sum(gone_ae_SAB_logging)*100/ np.sum(total_ae_SAB)) + "\n")
-gone_ou_SAW_all = np.where(#(abundance_data_SAW > 0.1) &
+gone_ou_SAW_all = np.where(
((resource_use_SAW == 1) |
(resource_use_SAW == 2) |
(resource_use_SAW == 3)), abundance_data_SAW, 0)
# # how much "area" was affected
total_ae_SAW = np.where(abundance_data_SAW > 0, 1, 0)
+
+
+
gone_ae_SAW_all = np.where(
- # (abundance_data_SAW > 0.1) &
+ (total_ae_SAW == 1) &
((resource_use_SAW == 1) |
(resource_use_SAW == 2) |
(resource_use_SAW == 3)), 1, 0)
+
+
+
affected.write( "all" + "," + "SAW" + "," +
str(np.sum(gone_ou_SAW_all) ) + "," +
@@ -872,13 +912,13 @@ affected.write( "all" + "," + "SAW" + "," +
str(np.sum(gone_ae_SAW_all)*100/ np.sum(total_ae_SAW)) + "\n")
# how many OU have been affected
-gone_ou_SAW_defor = np.where(#(abundance_data_SAW > 0.1) &
+gone_ou_SAW_defor = np.where(
((resource_use_SAW == 1) |
(resource_use_SAW == 2)), abundance_data_SAW, 0)
# # how much "area" was affected
gone_ae_SAW_defor = np.where(
- # (abundance_data_SAW > 0.1) &
+ (total_ae_SAW == 1) &
((resource_use_SAW == 1) |
(resource_use_SAW == 2)), 1, 0)
@@ -891,9 +931,8 @@ affected.write( "cover_change" + "," + "SAW" + "," +
-gone_ou_SAW_logging = np.where(#(abundance_data_SAW > 0.1) &
- (resource_use_SAW == 3), abundance_data_SAW, 0)
-gone_ae_SAW_logging = np.where(#(abundance_data_SAW > 0.1) &
+gone_ou_SAW_logging = np.where((resource_use_SAW == 3), abundance_data_SAW, 0)
+gone_ae_SAW_logging = np.where((total_ae_SAW == 1) &
(resource_use_SAW == 3), 1, 0)