diff --git a/NASA_CMS_2023/Mexico/DPS_algorithm/GMB_Mexico.R b/NASA_CMS_2023/Mexico/DPS_algorithm/GMB_Mexico.R new file mode 100644 index 0000000000000000000000000000000000000000..9a46f3d3aead266d93bd946ac2a67d48502217a0 --- /dev/null +++ b/NASA_CMS_2023/Mexico/DPS_algorithm/GMB_Mexico.R @@ -0,0 +1,97 @@ +############# OPEN R, THEN RUN THE FOLLOWING COMMAND ################### + +options(repos=c(CRAN="https://cran.r-project.org")) +install.packages("INLA",repos=c(getOption("repos"),INLA="https://inla.r-inla-download.org/R/stable"), dep=TRUE) + +############## LOAD PACKAGES ########### + +library("fmesher") +library(MatrixModels) +library(Matrix) +library(INLA) +library(inlabru) +library(sf) +library(terra) +library(dplyr) +library(spdep) +library(exactextractr) +library(sn) + +######################################## +args <- commandArgs(trailingOnly = TRUE) + +#### LOAD THE DATA USED TO MAKE THE MODEL ###### +mexico = st_read(args[2], quiet = TRUE) %>% st_union() %>% st_transform(crs = 6933) # +cci.rast = rast(args[3]) +hei.rast = rast(args[4]) + +DATA <- read.csv(args[5]) #Read the saved data +cci.plot <- DATA$cci.plot +hei.plot <- DATA$hei.plot +loc.plot <- data.frame(matrix(ncol = 2, nrow = length(hei.plot))) +loc.plot[,1] <- DATA$X +loc.plot[,2] <- DATA$X.1 +agbd.plot <- sqrt(DATA$agbd.plot) +loc.plot <- data.matrix(loc.plot) + +###### REMAKE THE MESH IN CASE IT IS NOT STORED IN MEMORY ###### +max.edge = 10*10^3 +mexico.buffer = st_buffer(mexico, dist = max.edge*5) +mesh = inla.mesh.2d(boundary = list(as(mexico, "Spatial"), as(mexico.buffer, "Spatial")), max.edge = c(max.edge, 3*max.edge), cutoff = 2*max.edge/3, offset = c(max.edge, 5*max.edge)) +k = mesh$n +A.plot = inla.spde.make.A(mesh = mesh, loc = loc.plot) + +##### LOAD THE MODEL RESULTS #### +load(args[6]) +samples = inla.posterior.sample(n = 250, result = model_fit_v2) + +#### LOAD AN URBAN MASK AND PROBABILITY OF FOREST MASK FOR MEXICO #### +URBAN = st_read(args[7],quiet=TRUE) +FNF = rast(args[8]) + +##### MAKE PREDICTIONS ########### + +PROJECTS = st_read(args[9], quiet = TRUE) %>% st_transform(crs = 6933) +PROJECTS$ID <- seq.int(nrow(PROJECTS)) +PROJECTS$GMB_mean_AGBD_dense = NA +PROJECTS$GMB_sd_AGBD_dense = NA + +for (f in (1:nrow(PROJECTS))) { + PROJECTS_ID = PROJECTS[f,]$ID + grid_to_predict = st_make_grid(PROJECTS[f,], cellsize = c(100,100), what = "centers") %>% st_as_sf() %>% st_filter(PROJECTS[f,]) + grid_to_predict = grid_to_predict[URBAN, ,op=st_disjoint] + + if (nrow(grid_to_predict[URBAN, ,op=st_disjoint])>0){ + grid = grid_to_predict %>% st_coordinates() + cci.pred = extract(cci.rast, grid)[,1] + hei.pred = extract(hei.rast, grid)[,1] + + if (length(cci.pred) > 0 && length(hei.pred) > 0){ + A.pred = inla.spde.make.A(mesh = mesh, loc = grid) + pred_fun = function(...){ + drop(intercept + + cci.pred*cci + + hei.pred*hei + + A.pred%*%alpha.spat[1:k] + + Diagonal(x = cci.pred)%*%A.pred%*%beta.spat[1:k] + + Diagonal(x = hei.pred)%*%A.pred%*%eta.spat[1:k]) + + rnorm(nrow(A.pred), sd = sqrt(1/theta[1])) + } + + # Generate prediction samples + pred.samples = (inla.posterior.sample.eval(fun = pred_fun, samples = samples)) + pred.samples[pred.samples < 0] = 0 + + # Model expectations and SD's at the grid locations + pred.mu_dense = Matrix::rowMeans(pred.samples^2, na.rm = TRUE) + pred.sd_dense = apply(pred.samples^2, 1, sd, na.rm = TRUE) + + # Enter information to CSV file + PROJECTS$GMB_mean_AGBD_dense[PROJECTS$ID == PROJECTS_ID] <- round(mean(pred.mu_dense,na.rm=TRUE),digits=2) + PROJECTS$GMB_sd_AGBD_dense[PROJECTS$ID == PROJECTS_ID] <- round(sd(colMeans(pred.samples^2, na.rm = T),na.rm=TRUE),digits=2) + } + } +} + +#### SAVE OUR RESULTS +write.csv(st_drop_geometry(PROJECTS),paste0(args[1])) \ No newline at end of file diff --git a/NASA_CMS_2023/Mexico/DPS_algorithm/build_command_main_ADE.sh b/NASA_CMS_2023/Mexico/DPS_algorithm/build_command_main_ADE.sh new file mode 100644 index 0000000000000000000000000000000000000000..bb7586250950208bf2ac39f7c5cc3abc38ac7195 --- /dev/null +++ b/NASA_CMS_2023/Mexico/DPS_algorithm/build_command_main_ADE.sh @@ -0,0 +1,16 @@ +#!/bin/bash +set -x +basedir=$( cd "$(dirname "$0")" ; pwd -P ) +conda env update -n base --solver=libmamba -f ${basedir}/env_main_ADE.yaml + +conda install conda-forge::r-fmesher --yes +conda install conda-forge::r-exactextractr --yes +conda install conda-forge::r-sn --yes +conda install conda-forge::r-MatrixModels --yes +conda install conda-forge::r-inlabru --yes +conda install conda-forge::r-terra --yes +conda install conda-forge::r-dplyr --yes +conda install conda-forge::r-spdep --yes +conda install conda-forge::r-sf --yes +conda install conda-forge::r-fields --yes +conda install conda-forge::r-Matrix --yes \ No newline at end of file diff --git a/NASA_CMS_2023/Mexico/DPS_algorithm/env_main_ADE.yaml b/NASA_CMS_2023/Mexico/DPS_algorithm/env_main_ADE.yaml new file mode 100644 index 0000000000000000000000000000000000000000..75474a878b0d65a3277bab71f00845166a63fcb3 --- /dev/null +++ b/NASA_CMS_2023/Mexico/DPS_algorithm/env_main_ADE.yaml @@ -0,0 +1,16 @@ +name: Packages_for_GMB +channels: + - conda-forge + - defaults +dependencies: + - r-fmesher + - r-exactextractr + - r-sn + - r-MatrixModels + - r-inlabru + - r-terra + - r-dplyr + - r-spdep + - r-sf + - r-fields + - r-Matrix \ No newline at end of file diff --git a/NASA_CMS_2023/Mexico/DPS_algorithm/run_GMBpredictions_MEX.sh b/NASA_CMS_2023/Mexico/DPS_algorithm/run_GMBpredictions_MEX.sh new file mode 100644 index 0000000000000000000000000000000000000000..7b09c0195d8fa26ed1cf0638e01798c082ef7948 --- /dev/null +++ b/NASA_CMS_2023/Mexico/DPS_algorithm/run_GMBpredictions_MEX.sh @@ -0,0 +1,19 @@ +#!/bin/bash + +basedir=$( cd "$(dirname "$0")" ; pwd -P ) + +mkdir -p output + +output_csv=${1} +mexico="input/ecort08gw_DESECON1_DISS.gpkg" +ccirast="input/cci_mexico_6933.tif" +heirast="input/GLAD_FH_mexico_UINT16_C_6399.tif" +data="input/NFI_CCI_GEDIheights.csv" +INLA_model_fit_v2="input/INLA_model_fit_v2.RData" +MGN2020_INEGI_Urban_Mex_mask_DISS="input/MGN2020_INEGI_Urban_Mex_mask_DISS.gpkg" +All_Products_Comp_Over30_Binary_6933="input/All_Products_Comp_Over30_Binary_6933.tif" +Projects="input/Projects.gpkg" + +OUTPUTDIR="${PWD}/output" + +Rscript ${basedir}/GMB_Mexico.R ${output_csv} ${mexico} ${ccirast} ${heirast} ${data} ${INLA_model_fit_v2} ${MGN2020_INEGI_Urban_Mex_mask_DISS} ${All_Products_Comp_Over30_Binary_6933} ${Projects} \ No newline at end of file diff --git a/NASA_CMS_2023/Mexico/DPS_algorithm/run_GMBpredictions_MEX.yaml b/NASA_CMS_2023/Mexico/DPS_algorithm/run_GMBpredictions_MEX.yaml new file mode 100644 index 0000000000000000000000000000000000000000..e1b33f9fb5b5a043e2c1f700740d0f63740a7b14 --- /dev/null +++ b/NASA_CMS_2023/Mexico/DPS_algorithm/run_GMBpredictions_MEX.yaml @@ -0,0 +1,50 @@ +algorithm_description: Runs GMB predictions for Mexico +algorithm_name: run_GMBpredictions_MEX +algorithm_version: GMBpredictions_v1 +repository_url: https://repo.maap-project.org/lduncanson/biomass_harmonization +docker_container_url: mas.maap-project.org/root/maap-workspaces/base_images/vanilla:v3.1.3 +queue: maap-dps-worker-16gb +build_command: biomass_harmonization/NASA_CMS_2023/Mexico/DPS_algorithm/build_command_main_ADE.sh +run_command: biomass_harmonization/NASA_CMS_2023/Mexico/DPS_algorithm/run_GMBpredictions_MEX.sh +disk_space: 20GB + +inputs: + config: [] + file: + - default: 'mexico' + description: mexico + name: mexico + required: true + - default: 'ccirast' + description: ccirast + name: ccirast + required: true + - default: 'heirast' + description: heirast + name: heirast + required: true + - default: 'data' + description: data + name: data + required: true + - default: 'INLA_model_fit_v2' + description: INLA_model_fit_v2 + name: INLA_model_fit_v2 + required: true + - default: 'MGN2020_INEGI_Urban_Mex_mask_DISS' + description: MGN2020_INEGI_Urban_Mex_mask_DISS + name: MGN2020_INEGI_Urban_Mex_mask_DISS + required: true + - default: 'All_Products_Comp_Over30_Binary_6933' + description: All_Products_Comp_Over30_Binary_6933 + name: All_Products_Comp_Over30_Binary_6933 + required: true + - default: 'Projects' + description: Projects + name: Projects + required: true + positional: + - default: 'output_csv' + description: output_csv + name: output_csv + required: true \ No newline at end of file