Skip to content
Snippets Groups Projects
Commit 55c2bfc4 authored by Neha Hunka's avatar Neha Hunka
Browse files

prep for DPS

parent 475d96e5
No related branches found
No related tags found
No related merge requests found
############# 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
#!/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
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
#!/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
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
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment