# This R-script calculates the C stock for each mapping unit of the 
# Harmonized World Soil Database (HWSD, v.1.2, 2012) and produces a recoding 
# table that can be used by GRASS GIS to convert a map of the mapping 
# units to a map of C stocks and the areal share of histosols
# and gelic phases (using the command r.recode).
# The script was written for R 2.13.0 
#
# R Development Core Team (2011). R: A language and environment for statistical
#   computing. R Foundation for Statistical Computing, Vienna, Austria. ISBN
#   3-900051-07-0, URL http://www.R-project.org/.
# GRASS Development Team, 2011. Geographic Resources Analysis Support System 
#   (GRASS) Software, Version 6.4.1. Open Source Geospatial Foundation, Available 
#   from http://grass.osgeo.org.
# FAO, IIASA, ISRIC, ISSCAS, JRC (2012) Harmonized World Soil Database (version 1.2). 
#   FAO and IIASA, Rome, Italy, and Laxenburg, Austria.
#
# The last lines of the script as it is published are removed here, because they were not 
#  a necessary part and required too much memory of R to be processed.

# Gobal soil mass after each step (HWSD 1.1):
# --- 2469.464
# -0- 2469.464
# -1- 2469.464
# -2- 2469.464
# -3- 2470.638
# -6- 2471.133
# -7- 2471.322
# -8c 1230.218 >3% OC function
# -8a 1113.296 only histosol BSD = 0.1
# -8b 1060.874 >3% OC function & histosol BSD = 0.1

### prepare data
# 1. obtain the HWSD from http://webarchive.iiasa.ac.at/Research/LUC/External-World-soil-database/HWSD_Viewer/HWSD_viewer_setup.exe
# 2. manually export the tables of the HWSD as tab-delimited text
# 3. the next steps assume that the table are stored in the local directory ~/HWSD

setwd("~/HWSD") # set the path to the directory containing the exported tables
hwsd.att=read.table("HWSD_DATA.txt", header=TRUE, sep=";")

hwsd.att$T_BSDadj = hwsd.att$T_REF_BULK_DENSITY # gap filling and correction to soils with OC>3%
hwsd.att$S_BSDadj = hwsd.att$S_REF_BULK_DENSITY
hwsd.att$T_BSDadj2 = hwsd.att$T_REF_BULK_DENSITY # gap filling and BD correction to soils with OC>3% and HIST: BD=0.1 
hwsd.att$S_BSDadj2 = hwsd.att$S_REF_BULK_DENSITY
hwsd.att$T_BSDadj3 = hwsd.att$T_REF_BULK_DENSITY # gap filling and using bulk density when provided, else use ref bulk density 
hwsd.att$S_BSDadj3 = hwsd.att$S_REF_BULK_DENSITY
hwsd.att$T_BSDadj4 = hwsd.att$T_REF_BULK_DENSITY # gap filling and use Hiederer-correction and and HIST: BD=0.1
hwsd.att$S_BSDadj4 = hwsd.att$S_REF_BULK_DENSITY
hwsd.att$T_OCadj = hwsd.att$T_OC
hwsd.att$S_OCadj = hwsd.att$S_OC

hwsd.att$SEQ[hwsd.att$MU_GLOBAL==31651]<-1:3 # replace double entries
hwsd.att$SEQ[hwsd.att$MU_GLOBAL==7000]<-1 # add missing value

### Adjustment of gaps and inconsistencies.
# The adjustments have been suggested by R. Hiederer, JRC.
###

# 0. Replace ISSOIL where value is 0 although HSf has soil characteristics:
# MU 7004, 12972, 16608, 16995, 26851

hwsd.att$ISSOILadj = ifelse(hwsd.att$MU_GLOBAL==7004 | hwsd.att$MU_GLOBAL==12972 | hwsd.att$MU_GLOBAL==16608 | hwsd.att$MU_GLOBAL==16995 | hwsd.att$MU_GLOBAL==26851, 1, hwsd.att$ISSOIL)

# 1. Assign OC to FAO90 TOP_OC without value:
hwsd.att$T_OCadj = ifelse(hwsd.att$SU_SYM90=='CMx' & (hwsd.att$T_OC==0 | is.na(hwsd.att$T_OC)), 0.93, hwsd.att$T_OCadj) # 2 cases
hwsd.att$T_OCadj = ifelse(hwsd.att$SU_SYM90=='ARo' & (hwsd.att$T_OC==0 | is.na(hwsd.att$T_OC)), 0.27, hwsd.att$T_OCadj) # 1 case, median of remaining
hwsd.att$T_OCadj = ifelse(hwsd.att$SU_SYM90=='HSf' & (hwsd.att$T_OC==0 | is.na(hwsd.att$T_OC)), 33.63, hwsd.att$T_OCadj) # 1 case

# 2. Assign OC to FAO90 SUB_OC without value:
hwsd.att$S_OCadj = ifelse(hwsd.att$SU_SYM90=='ARg' & (hwsd.att$S_OC==0 | is.na(hwsd.att$S_OC)), 0.12, hwsd.att$S_OCadj) # 4 cases, median of remaining
hwsd.att$S_OCadj = ifelse(hwsd.att$SU_SYM90=='ARo' & (hwsd.att$S_OC==0 | is.na(hwsd.att$S_OC)), 0.11, hwsd.att$S_OCadj) # 4 cases
hwsd.att$S_OCadj = ifelse(hwsd.att$SU_SYM90=='CMx' & (hwsd.att$S_OC==0 | is.na(hwsd.att$S_OC)), 0.38, hwsd.att$S_OCadj) # 6 cases
hwsd.att$S_OCadj = ifelse(hwsd.att$SU_SYM90=='LXf' & (hwsd.att$S_OC==0 | is.na(hwsd.att$S_OC)), 0.32, hwsd.att$S_OCadj) # 2 cases
hwsd.att$S_OCadj = ifelse(hwsd.att$SU_SYM90=='HSf' & (hwsd.att$S_OC==0 | is.na(hwsd.att$S_OC)), 32.89, hwsd.att$S_OCadj) # 1 case

# 3. Assign OC values to blank FAO90 SUB_OC values (107 cases), 
# R. Hiederer used different replacement values
hwsd.att$S_OCadj = ifelse(hwsd.att$SU_SYM90=='CLl' & is.na(hwsd.att$S_OC) & is.na(hwsd.att$S_REF_BULK_DENSITY), 0.2, hwsd.att$S_OCadj)#  nearest value
hwsd.att$S_BSDadj = ifelse(hwsd.att$SU_SYM90=='CLl' & is.na(hwsd.att$S_OC) & is.na(hwsd.att$S_REF_BULK_DENSITY), 1.36, hwsd.att$S_BSDadj)
hwsd.att$S_OCadj = ifelse(hwsd.att$SU_SYM90=='ACg' & is.na(hwsd.att$S_OC) & is.na(hwsd.att$S_REF_BULK_DENSITY), 0.5, hwsd.att$S_OCadj) # 1 case, nearest value
hwsd.att$S_BSDadj = ifelse(hwsd.att$SU_SYM90=='ACg' & is.na(hwsd.att$S_OC) & is.na(hwsd.att$S_REF_BULK_DENSITY), 1.42, hwsd.att$S_BSDadj)
hwsd.att$S_OCadj = ifelse(hwsd.att$SU_SYM90=='ALg' & is.na(hwsd.att$S_OC) & is.na(hwsd.att$S_REF_BULK_DENSITY), 0.34, hwsd.att$S_OCadj)# 1 case, nearest value
hwsd.att$S_BSDadj = ifelse(hwsd.att$SU_SYM90=='ALg' & is.na(hwsd.att$S_OC) & is.na(hwsd.att$S_REF_BULK_DENSITY), 1.38, hwsd.att$S_BSDadj)
hwsd.att$S_OCadj = ifelse(hwsd.att$SU_SYM90=='ANz' & is.na(hwsd.att$S_OC) & is.na(hwsd.att$S_REF_BULK_DENSITY), 3.76, hwsd.att$S_OCadj)# 5 cases, nearest value
hwsd.att$S_BSDadj = ifelse(hwsd.att$SU_SYM90=='ANz' & is.na(hwsd.att$S_OC) & is.na(hwsd.att$S_REF_BULK_DENSITY), 1.63, hwsd.att$S_BSDadj)
hwsd.att$S_BSDadj = ifelse(hwsd.att$SU_SYM90=='ANz' & hwsd.att$S_OC==5.2 & is.na(hwsd.att$S_REF_BULK_DENSITY), 1.63, hwsd.att$S_BSDadj)
hwsd.att$S_OCadj = ifelse(hwsd.att$SU_SYM90=='CLh' & is.na(hwsd.att$S_OC) & is.na(hwsd.att$S_REF_BULK_DENSITY), 0.22, hwsd.att$S_OCadj)# 1 case, nearest value
hwsd.att$S_BSDadj = ifelse(hwsd.att$SU_SYM90=='CLh' & is.na(hwsd.att$S_OC) & is.na(hwsd.att$S_REF_BULK_DENSITY), 1.36, hwsd.att$S_BSDadj)
hwsd.att$S_OCadj = ifelse(hwsd.att$SU_SYM90=='CLp' & is.na(hwsd.att$S_OC) & is.na(hwsd.att$S_REF_BULK_DENSITY), 0.34, hwsd.att$S_OCadj)# 9 cases, nearest value
hwsd.att$S_BSDadj = ifelse(hwsd.att$SU_SYM90=='CLp' & is.na(hwsd.att$S_OC) & is.na(hwsd.att$S_REF_BULK_DENSITY), 1.42, hwsd.att$S_BSDadj)
hwsd.att$S_OCadj = ifelse(hwsd.att$MU_GLOBAL==18631, 0.23, hwsd.att$S_OCadj)# replaces one of the above values, nearest value
hwsd.att$S_BSDadj = ifelse(hwsd.att$MU_GLOBAL==18631, 1.44, hwsd.att$S_BSDadj)
hwsd.att$S_OCadj = ifelse(hwsd.att$SU_SYM90=='CMd' & is.na(hwsd.att$S_OC) & is.na(hwsd.att$S_REF_BULK_DENSITY), 0.45, hwsd.att$S_OCadj)# 1 case, nearest value
hwsd.att$S_BSDadj = ifelse(hwsd.att$SU_SYM90=='CMd' & is.na(hwsd.att$S_OC) & is.na(hwsd.att$S_REF_BULK_DENSITY), 1.29, hwsd.att$S_BSDadj)
hwsd.att$S_OCadj = ifelse(hwsd.att$SU_SYM90=='CMo' & is.na(hwsd.att$S_OC) & is.na(hwsd.att$S_REF_BULK_DENSITY), 0.2, hwsd.att$S_OCadj)# 1 case, nearest value
hwsd.att$S_BSDadj = ifelse(hwsd.att$SU_SYM90=='CMo' & is.na(hwsd.att$S_OC) & is.na(hwsd.att$S_REF_BULK_DENSITY), 1.44, hwsd.att$S_BSDadj)
hwsd.att$S_OCadj = ifelse(hwsd.att$SU_SYM90=='CMv' & is.na(hwsd.att$S_OC) & is.na(hwsd.att$S_REF_BULK_DENSITY), 0.24, hwsd.att$S_OCadj)# 1 case, nearest value
hwsd.att$S_BSDadj = ifelse(hwsd.att$SU_SYM90=='CMv' & is.na(hwsd.att$S_OC) & is.na(hwsd.att$S_REF_BULK_DENSITY), 1.31, hwsd.att$S_BSDadj)
hwsd.att$S_OCadj = ifelse(hwsd.att$SU_SYM90=='CMx' & is.na(hwsd.att$S_OC) & is.na(hwsd.att$S_REF_BULK_DENSITY), 0.41, hwsd.att$S_OCadj)# 4 cases, nearest value
hwsd.att$S_BSDadj = ifelse(hwsd.att$SU_SYM90=='CMx' & is.na(hwsd.att$S_OC) & is.na(hwsd.att$S_REF_BULK_DENSITY), 1.42, hwsd.att$S_BSDadj)
hwsd.att$S_OCadj = ifelse(hwsd.att$SU_SYM90=='FLd' & is.na(hwsd.att$S_OC) & is.na(hwsd.att$S_REF_BULK_DENSITY), 0.19, hwsd.att$S_OCadj)# 1 case, nearest value
hwsd.att$S_BSDadj = ifelse(hwsd.att$SU_SYM90=='FLd' & is.na(hwsd.att$S_OC) & is.na(hwsd.att$S_REF_BULK_DENSITY), 1.46, hwsd.att$S_BSDadj)
hwsd.att$S_OCadj = ifelse(hwsd.att$SU_SYM90=='FLe' & is.na(hwsd.att$S_OC) & is.na(hwsd.att$S_REF_BULK_DENSITY), 0.23, hwsd.att$S_OCadj)# 6 cases, top value
hwsd.att$S_BSDadj = ifelse(hwsd.att$SU_SYM90=='FLe' & is.na(hwsd.att$S_OC) & is.na(hwsd.att$S_REF_BULK_DENSITY), 1.26, hwsd.att$S_BSDadj)
hwsd.att$S_OCadj = ifelse(hwsd.att$SU_SYM90=='FRh' & is.na(hwsd.att$S_OC) & is.na(hwsd.att$S_REF_BULK_DENSITY), 0.3, hwsd.att$S_OCadj)# 1 case, nearest value
hwsd.att$S_BSDadj = ifelse(hwsd.att$SU_SYM90=='FRh' & is.na(hwsd.att$S_OC) & is.na(hwsd.att$S_REF_BULK_DENSITY), 1.48, hwsd.att$S_BSDadj)
hwsd.att$S_OCadj = ifelse(hwsd.att$SU_SYM90=='GLe' & is.na(hwsd.att$S_OC) & is.na(hwsd.att$S_REF_BULK_DENSITY), 0.42, hwsd.att$S_OCadj)# 3 cases, nearest value
hwsd.att$S_BSDadj = ifelse(hwsd.att$SU_SYM90=='GLe' & is.na(hwsd.att$S_OC) & is.na(hwsd.att$S_REF_BULK_DENSITY), 1.29, hwsd.att$S_BSDadj)
hwsd.att$S_OCadj = ifelse(hwsd.att$SU_SYM90=='GLm' & is.na(hwsd.att$S_OC) & is.na(hwsd.att$S_REF_BULK_DENSITY), 0.55, hwsd.att$S_OCadj)# 3 cases, nearest value
hwsd.att$S_BSDadj = ifelse(hwsd.att$SU_SYM90=='GLm' & is.na(hwsd.att$S_OC) & is.na(hwsd.att$S_REF_BULK_DENSITY), 1.31, hwsd.att$S_BSDadj)
hwsd.att$S_OCadj = ifelse(hwsd.att$SU_SYM90=='GYk' & is.na(hwsd.att$S_OC) & is.na(hwsd.att$S_REF_BULK_DENSITY), 0.18, hwsd.att$S_OCadj)# 1 case, nearest value
hwsd.att$S_BSDadj = ifelse(hwsd.att$SU_SYM90=='GYk' & is.na(hwsd.att$S_OC) & is.na(hwsd.att$S_REF_BULK_DENSITY), 1.43, hwsd.att$S_BSDadj)
hwsd.att$S_OCadj = ifelse(hwsd.att$SU_SYM90=='GYl' & is.na(hwsd.att$S_OC) & is.na(hwsd.att$S_REF_BULK_DENSITY), 0.13, hwsd.att$S_OCadj)# 1 case, nearest value
hwsd.att$S_BSDadj = ifelse(hwsd.att$SU_SYM90=='GYl' & is.na(hwsd.att$S_OC) & is.na(hwsd.att$S_REF_BULK_DENSITY), 1.36, hwsd.att$S_BSDadj)
hwsd.att$S_OCadj = ifelse(hwsd.att$SU_SYM90=='KSh' & is.na(hwsd.att$S_OC) & is.na(hwsd.att$S_REF_BULK_DENSITY), 0.46, hwsd.att$S_OCadj)# 6 cases, nearest values
hwsd.att$S_BSDadj = ifelse(hwsd.att$SU_SYM90=='KSh' & is.na(hwsd.att$S_OC) & is.na(hwsd.att$S_REF_BULK_DENSITY), 1.32, hwsd.att$S_BSDadj)
hwsd.att$S_BSDadj = ifelse(hwsd.att$MU_GLOBAL==12257, 1.32, hwsd.att$S_BSDadj)
hwsd.att$S_OCadj = ifelse(hwsd.att$SU_SYM90=='KSk' & is.na(hwsd.att$S_OC) & is.na(hwsd.att$S_REF_BULK_DENSITY), 0.52, hwsd.att$S_OCadj)# 1 case, nearest value
hwsd.att$S_BSDadj = ifelse(hwsd.att$SU_SYM90=='KSk' & is.na(hwsd.att$S_OC) & is.na(hwsd.att$S_REF_BULK_DENSITY), 1.32, hwsd.att$S_BSDadj)
hwsd.att$S_OCadj = ifelse(hwsd.att$SU_SYM90=='LVh' & is.na(hwsd.att$S_OC) & is.na(hwsd.att$S_REF_BULK_DENSITY), 0.71, hwsd.att$S_OCadj)# 1 case, nearest value
hwsd.att$S_BSDadj = ifelse(hwsd.att$SU_SYM90=='LVh' & is.na(hwsd.att$S_OC) & is.na(hwsd.att$S_REF_BULK_DENSITY), 1.32, hwsd.att$S_BSDadj)
hwsd.att$S_OCadj = ifelse(hwsd.att$SU_SYM90=='LVx' & is.na(hwsd.att$S_OC) & is.na(hwsd.att$S_REF_BULK_DENSITY), 0.78, hwsd.att$S_OCadj)# 18 cases, nearest value
hwsd.att$S_BSDadj = ifelse(hwsd.att$SU_SYM90=='LVx' & is.na(hwsd.att$S_OC) & is.na(hwsd.att$S_REF_BULK_DENSITY), 1.29, hwsd.att$S_BSDadj)
hwsd.att$S_OCadj = ifelse(hwsd.att$MU_GLOBAL==18631, 0.39, hwsd.att$S_OCadj)# replaces one of the above values, nearest value
hwsd.att$S_BSDadj = ifelse(hwsd.att$MU_GLOBAL==18631, 1.34, hwsd.att$S_BSDadj)
hwsd.att$S_OCadj = ifelse(hwsd.att$MU_GLOBAL==29645 | hwsd.att$MU_GLOBAL==29646, 0.33, hwsd.att$S_OCadj)
hwsd.att$S_BSDadj = ifelse(hwsd.att$MU_GLOBAL==29645| hwsd.att$MU_GLOBAL==29646, 1.32, hwsd.att$S_BSDadj)
hwsd.att$S_OCadj = ifelse(hwsd.att$SU_SYM90=='PHh' & is.na(hwsd.att$S_OC) & is.na(hwsd.att$S_REF_BULK_DENSITY), 0.11, hwsd.att$S_OCadj)# 1 case, nearest value
hwsd.att$S_BSDadj = ifelse(hwsd.att$SU_SYM90=='PHh' & is.na(hwsd.att$S_OC) & is.na(hwsd.att$S_REF_BULK_DENSITY), 1.41, hwsd.att$S_BSDadj)
hwsd.att$S_OCadj = ifelse(hwsd.att$SU_SYM90=='PHl' & is.na(hwsd.att$S_OC) & is.na(hwsd.att$S_REF_BULK_DENSITY), 0.76, hwsd.att$S_OCadj)# 1 case, nearest value
hwsd.att$S_BSDadj = ifelse(hwsd.att$SU_SYM90=='PHl' & is.na(hwsd.att$S_OC) & is.na(hwsd.att$S_REF_BULK_DENSITY), 1.19, hwsd.att$S_BSDadj)
hwsd.att$S_OCadj = ifelse(hwsd.att$SU_SYM90=='PLd' & is.na(hwsd.att$S_OC) & is.na(hwsd.att$S_REF_BULK_DENSITY), 0.34, hwsd.att$S_OCadj)# 5 cases, nearest value
hwsd.att$S_BSDadj = ifelse(hwsd.att$SU_SYM90=='PLd' & is.na(hwsd.att$S_OC) & is.na(hwsd.att$S_REF_BULK_DENSITY), 1.32, hwsd.att$S_BSDadj)
hwsd.att$S_OCadj = ifelse(hwsd.att$SU_SYM90=='RGd' & is.na(hwsd.att$S_OC) & is.na(hwsd.att$S_REF_BULK_DENSITY), 1.08, hwsd.att$S_OCadj)# 2 cases, nearest value
hwsd.att$S_BSDadj = ifelse(hwsd.att$SU_SYM90=='RGd' & is.na(hwsd.att$S_OC) & is.na(hwsd.att$S_REF_BULK_DENSITY), 1.27, hwsd.att$S_BSDadj)
hwsd.att$S_OCadj = ifelse(hwsd.att$MU_GLOBAL==18325, 1.87, hwsd.att$S_OCadj)
hwsd.att$S_BSDadj = ifelse(hwsd.att$MU_GLOBAL==18325, 1.45, hwsd.att$S_BSDadj)
hwsd.att$S_OCadj = ifelse(hwsd.att$SU_SYM90=='SCh' & is.na(hwsd.att$S_OC) & is.na(hwsd.att$S_REF_BULK_DENSITY), 0.19, hwsd.att$S_OCadj)# 3 cases, nearest value
hwsd.att$S_BSDadj = ifelse(hwsd.att$SU_SYM90=='SCh' & is.na(hwsd.att$S_OC) & is.na(hwsd.att$S_REF_BULK_DENSITY), 1.42, hwsd.att$S_BSDadj)
hwsd.att$S_OCadj = ifelse(hwsd.att$SU_SYM90=='SCm' & is.na(hwsd.att$S_OC) & is.na(hwsd.att$S_REF_BULK_DENSITY), 0.5, hwsd.att$S_OCadj)# 2 cases, top value
hwsd.att$S_BSDadj = ifelse(hwsd.att$SU_SYM90=='SCm' & is.na(hwsd.att$S_OC) & is.na(hwsd.att$S_REF_BULK_DENSITY), 1.47, hwsd.att$S_BSDadj)
hwsd.att$S_OCadj = ifelse(hwsd.att$SU_SYM90=='SNh' & is.na(hwsd.att$S_OC) & is.na(hwsd.att$S_REF_BULK_DENSITY), 0.22, hwsd.att$S_OCadj)# 2 cases, top value
hwsd.att$S_BSDadj = ifelse(hwsd.att$SU_SYM90=='SNh' & is.na(hwsd.att$S_OC) & is.na(hwsd.att$S_REF_BULK_DENSITY), 1.34, hwsd.att$S_BSDadj)

#steps 4-5 by R. Hiederer skipped

# 6. Assign values to TOP_BD where blank and not histsols (23 cases)
hwsd.att$T_BSDadj = ifelse(substr(hwsd.att$SU_SYM90,1,2)=='AR' & is.na(hwsd.att$T_REF_BULK_DENSITY), 1.8, hwsd.att$T_REF_BULK_DENSITY) # 17 cases, value from "Lecture notes on the major soils of the world"
hwsd.att$T_BSDadj = ifelse(hwsd.att$SU_SYM90=='KSh' & is.na(hwsd.att$T_REF_BULK_DENSITY), 1.31, hwsd.att$T_BSDadj) # 1 case, median
hwsd.att$T_BSDadj = ifelse(hwsd.att$SU_SYM90=='LPm' & is.na(hwsd.att$T_REF_BULK_DENSITY), 1.5, hwsd.att$T_BSDadj) # 2 cases, regression
hwsd.att$T_BSDadj = ifelse(hwsd.att$SU_SYM90=='PZh' & is.na(hwsd.att$T_REF_BULK_DENSITY), 1.78, hwsd.att$T_BSDadj) # 1 case, nearest value
hwsd.att$T_BSDadj = ifelse(hwsd.att$SU_SYM90=='RGd' & is.na(hwsd.att$T_REF_BULK_DENSITY), 1.44, hwsd.att$T_BSDadj) # 1 case, nearest value
hwsd.att$T_BSDadj = ifelse(hwsd.att$SU_SYM90=='RGe' & is.na(hwsd.att$T_REF_BULK_DENSITY), 1.33, hwsd.att$T_BSDadj) # 1 case, nearest value

# 7. Assign values to SUB_BD where blank from TOP_BD (4 cases): MU_GLOBAL FAO90: # 12256	RGe, 12737	PZc, 12790	PZc, 13420	ANz, 31357	ARa
hwsd.att$S_BSDadj = ifelse(is.na(hwsd.att$S_BSDadj) & is.na(hwsd.att$S_OCadj)==FALSE & is.na(hwsd.att$T_BSDadj)==FALSE, hwsd.att$T_BSDadj, hwsd.att$S_BSDadj) # 5 cases

# 7a For use with HWSD 1.2: Use BULK_DENSITY from profiles if provided, else use corrected REF_BULK_DENSITY from previous step.
hwsd.att$T_BSDadj3 = ifelse(hwsd.att$T_BULK_DENSITY>0, hwsd.att$T_BULK_DENSITY, hwsd.att$T_BSDadj)
hwsd.att$S_BSDadj3 = ifelse(hwsd.att$S_BULK_DENSITY>0 | hwsd.att$S_BULK_DENSITY<10, hwsd.att$S_BULK_DENSITY, hwsd.att$S_BSDadj)

#8. Adjust BD for SOC > 3% to BD = f(OC); try different adjustments

hwsd.att$T_BSDadj = ifelse(hwsd.att$T_OCadj>3, -0.31*log(hwsd.att$T_OCadj)+1.38, hwsd.att$T_BSDadj) # kg/dm3, R. Hiederer, pers. comm.
hwsd.att$S_BSDadj = ifelse(hwsd.att$S_OCadj>3, -0.32*log(hwsd.att$S_OCadj)+1.38, hwsd.att$S_BSDadj) # kg/dm3, R. Hiederer, pers. comm.

hwsd.att$hist = ifelse(substr(hwsd.att$SU_SYM74,1,1)=='O'|substr(hwsd.att$SU_SYM90,1,2)=='HS', 1, 0) # kg/dm3
hwsd.att$T_BSDadj2 = ifelse(hwsd.att$hist, 0.1, hwsd.att$T_BSDadj) # kg/dm3
hwsd.att$S_BSDadj2 = ifelse(hwsd.att$hist, 0.1, hwsd.att$S_BSDadj) # kg/dm3



# calculate top layer C stock
# (if T_REF_BULK-DENSITY is missing (751 cases), no C stocks are calculated)
hwsd.att$T_OCG.stock = ifelse(is.na(hwsd.att$T_OC),0,hwsd.att$T_OC/100)*hwsd.att$T_REF_BULK_DENSITY*pmin(3, hwsd.att$REF_DEPTH/10)*(1-ifelse(is.na(hwsd.att$T_GRAVEL),0,hwsd.att$T_GRAVEL/100)) # result is in kg/(dm2*m)
hwsd.att$T_OCG.stock.BSDadj = ifelse(is.na(hwsd.att$T_OCadj),0,hwsd.att$T_OCadj/100)*hwsd.att$T_BSDadj*pmin(3, hwsd.att$REF_DEPTH/10)*(1-ifelse(is.na(hwsd.att$T_GRAVEL),0,hwsd.att$T_GRAVEL/100)) # result is in kg/(dm2*m)
hwsd.att$T_OCG.stock.BSDadj2 = ifelse(is.na(hwsd.att$T_OCadj),0,hwsd.att$T_OCadj/100)*hwsd.att$T_BSDadj2*pmin(3, hwsd.att$REF_DEPTH/10)*(1-ifelse(is.na(hwsd.att$T_GRAVEL),0,hwsd.att$T_GRAVEL/100)) # result is in kg/(dm2*m)
hwsd.att$T_OCG.stock.BSDadj3 = ifelse(is.na(hwsd.att$T_OCadj),0,hwsd.att$T_OCadj/100)*hwsd.att$T_BSDadj3*pmin(3, hwsd.att$REF_DEPTH/10)*(1-ifelse(is.na(hwsd.att$T_GRAVEL),0,hwsd.att$T_GRAVEL/100)) # result is in kg/(dm2*m)
hwsd.att$T_OCG.stock.BSDadj4 = ifelse(is.na(hwsd.att$T_OCadj),0,hwsd.att$T_OCadj/100)*hwsd.att$T_BSDadj4*pmin(3, hwsd.att$REF_DEPTH/10)*(1-ifelse(is.na(hwsd.att$T_GRAVEL),0,hwsd.att$T_GRAVEL/100)) # result is in kg/(dm2*m)

# calculate bottom layer C stock
S_OCG.stock=hwsd.att$S_OC/100 * hwsd.att$S_REF_BULK_DENSITY * pmax(0,(hwsd.att$REF_DEPTH-30)/10)*(1-ifelse(is.na(hwsd.att$S_GRAVEL),0,hwsd.att$S_GRAVEL/100))
S_OCG.stock.BSDadj=hwsd.att$S_OCadj/100 * hwsd.att$S_BSDadj * pmax(0,(hwsd.att$REF_DEPTH-30)/10)*(1-ifelse(is.na(hwsd.att$S_GRAVEL),0,hwsd.att$S_GRAVEL/100))
S_OCG.stock.BSDadj2=hwsd.att$S_OCadj/100 * hwsd.att$S_BSDadj2 * pmax(0,(hwsd.att$REF_DEPTH-30)/10)*(1-ifelse(is.na(hwsd.att$S_GRAVEL),0,hwsd.att$S_GRAVEL/100))
S_OCG.stock.BSDadj3=hwsd.att$S_OCadj/100 * hwsd.att$S_BSDadj3 * pmax(0,(hwsd.att$REF_DEPTH-30)/10)*(1-ifelse(is.na(hwsd.att$S_GRAVEL),0,hwsd.att$S_GRAVEL/100))
S_OCG.stock.BSDadj4=hwsd.att$S_OCadj/100 * hwsd.att$S_BSDadj4 * pmax(0,(hwsd.att$REF_DEPTH-30)/10)*(1-ifelse(is.na(hwsd.att$S_GRAVEL),0,hwsd.att$S_GRAVEL/100))

hwsd.att$S_OCG.stock = ifelse(is.na(S_OCG.stock),0,S_OCG.stock) # result is in kg/(dm2*m) 
hwsd.att$S_OCG.stock.BSDadj = ifelse(is.na(S_OCG.stock.BSDadj),0,S_OCG.stock.BSDadj) # result is in kg/(dm2*m) 
hwsd.att$S_OCG.stock.BSDadj2 = ifelse(is.na(S_OCG.stock.BSDadj2),0,S_OCG.stock.BSDadj2) # result is in kg/(dm2*m) 
hwsd.att$S_OCG.stock.BSDadj3 = ifelse(is.na(S_OCG.stock.BSDadj3),0,S_OCG.stock.BSDadj3) # result is in kg/(dm2*m) 
hwsd.att$S_OCG.stock.BSDadj4 = ifelse(is.na(S_OCG.stock.BSDadj4),0,S_OCG.stock.BSDadj4) # result is in kg/(dm2*m) 

# add top and bottom layer C stocks
hwsd.att$TS_OCG.stock = hwsd.att$T_OCG.stock + hwsd.att$S_OCG.stock
hwsd.att$TS_OCG.stock.BSDadj = hwsd.att$T_OCG.stock.BSDadj + hwsd.att$S_OCG.stock.BSDadj
hwsd.att$TS_OCG.stock.BSDadj2 = hwsd.att$T_OCG.stock.BSDadj2 + hwsd.att$S_OCG.stock.BSDadj2
hwsd.att$TS_OCG.stock.BSDadj3 = hwsd.att$T_OCG.stock.BSDadj3 + hwsd.att$S_OCG.stock.BSDadj3
hwsd.att$TS_OCG.stock.BSDadj4 = hwsd.att$T_OCG.stock.BSDadj4 + hwsd.att$S_OCG.stock.BSDadj4

# calculate the fractional share of soils with a gelic phase
hwsd.att$gelic = (ifelse(hwsd.att$ADD_PROP==2,1,0)*hwsd.att$SHARE/100)

# calculate the fraction of soil per pixel area
hwsd.soil.share = aggregate(hwsd.att$SHARE/100*hwsd.att$ISSOIL, by=list(MU_GLOBAL=hwsd.att$MU_GLOBAL), function(x){sum(x, na.rm=TRUE)}) # this accounts for regions where some portion of the soil is non-soil (ISSOIL=0)
write.table(hwsd.soil.share[hwsd.soil.share$x>0,c(1,1,2,2)], "HWSD_soilshare.txt", sep=":", row.names=FALSE, col.names=FALSE, quote=FALSE, na="0")
hwsd.soil.share.adj = aggregate(hwsd.att$SHARE/100*hwsd.att$ISSOILadj, by=list(MU_GLOBAL=hwsd.att$MU_GLOBAL), function(x){sum(x, na.rm=TRUE)}) # this accounts for regions where some portion of the soil is non-soil (ISSOIL=0)
write.table(hwsd.soil.share.adj[hwsd.soil.share.adj$x>0,c(1,1,2,2)], "HWSD_soilshareadj.txt", sep=":", row.names=FALSE, col.names=FALSE, quote=FALSE, na="0")

# sum the properties over each pixel fraction
hwsd.ocg.kg.m2 = aggregate((hwsd.att$TS_OCG.stock * hwsd.att$SHARE/100), by=list(MU_GLOBAL=hwsd.att$MU_GLOBAL), function(x){sum(x*100, na.rm=TRUE)}) # result is in kg/(m2*m)
hwsd.ocg.BSDadj.kg.m2 = aggregate((hwsd.att$TS_OCG.stock.BSDadj * hwsd.att$SHARE/100), by=list(MU_GLOBAL=hwsd.att$MU_GLOBAL), function(x){sum(x*100, na.rm=TRUE)}) # result is in kg/(m2*m)
hwsd.ocg.BSDadj2.kg.m2 = aggregate((hwsd.att$TS_OCG.stock.BSDadj2 * hwsd.att$SHARE/100), by=list(MU_GLOBAL=hwsd.att$MU_GLOBAL), function(x){sum(x*100, na.rm=TRUE)}) # result is in kg/(m2*m)
hwsd.ocg.BSDadj3.kg.m2 = aggregate((hwsd.att$TS_OCG.stock.BSDadj3 * hwsd.att$SHARE/100), by=list(MU_GLOBAL=hwsd.att$MU_GLOBAL), function(x){sum(x*100, na.rm=TRUE)}) # result is in kg/(m2*m)
hwsd.ocg.BSDadj4.kg.m2 = aggregate((hwsd.att$TS_OCG.stock.BSDadj4 * hwsd.att$SHARE/100), by=list(MU_GLOBAL=hwsd.att$MU_GLOBAL), function(x){sum(x*100, na.rm=TRUE)}) # result is in kg/(m2*m)
#### stop here to include proportion of bare soil ####

# The C stock as calculated here assumes that the whole pixel is covered by soil. 
# This is intended for maps of C stock, but for maps of soil mass per pixel 
#  C stock must be multiplied by pixel area and total soil share in a GIS.
hwsd.ocg.kg.m2$x=round(hwsd.ocg.kg.m2$x/hwsd.soil.share$x,3) # this overrepresents bare soil
hwsd.ocg.BSDadj.kg.m2$x=round(hwsd.ocg.BSDadj.kg.m2$x/hwsd.soil.share.adj$x,3) # this overrepresents bare soil
hwsd.ocg.BSDadj2.kg.m2$x=round(hwsd.ocg.BSDadj2.kg.m2$x/hwsd.soil.share.adj$x,3) # this overrepresents bare soil
hwsd.ocg.BSDadj3.kg.m2$x=round(hwsd.ocg.BSDadj3.kg.m2$x/hwsd.soil.share.adj$x,3) # this eliminates bare soil
hwsd.ocg.BSDadj4.kg.m2$x=round(hwsd.ocg.BSDadj4.kg.m2$x/hwsd.soil.share.adj$x,3) # this eliminates bare soil
write.table(hwsd.ocg.kg.m2[hwsd.soil.share$x>0,c(1,1,2,2)], "HWSD_TSOCG_kgm2.txt", sep=":", row.names=FALSE, col.names=FALSE, quote=FALSE, na="0")
write.table(hwsd.ocg.BSDadj.kg.m2[hwsd.soil.share.adj$x>0,c(1,1,2,2)], "HWSD_TSOCG_BSDadj1_kgm2.txt", sep=":", row.names=FALSE, col.names=FALSE, quote=FALSE, na="0")
write.table(hwsd.ocg.BSDadj2.kg.m2[hwsd.soil.share.adj$x>0,c(1,1,2,2)], "HWSD_TSOCG_BSDadj_kgm2.txt", sep=":", row.names=FALSE, col.names=FALSE, quote=FALSE, na="0")

# The fraction of soil per pixel being a histosol
hwsd.hist = aggregate((hwsd.att$hist * hwsd.att$SHARE/100), by=list(MU_GLOBAL=hwsd.att$MU_GLOBAL), function(x){sum(x, na.rm=TRUE)})
hwsd.hist$x=round(hwsd.hist$x/hwsd.soil.share.adj$x,3)
write.table(hwsd.hist[hwsd.soil.share.adj$x>0,c(1,1,2,2)], "HWSD_HISTadj.txt", sep=":", row.names=FALSE, col.names=FALSE, quote=FALSE, na="0")

# The fraction of soil per pixel having a gelic phase
hwsd.gelic = aggregate(hwsd.att$gelic, by=list(MU_GLOBAL=hwsd.att$MU_GLOBAL), function(x){sum(x, na.rm=TRUE)})
hwsd.gelic$x = round(hwsd.gelic$x/hwsd.soil.share.adj$x,2)
write.table(hwsd.gelic[hwsd.soil.share.adj$x>0,c(1,1,2,2)], "HWSD_GELICadj.txt", sep=":", row.names=FALSE, col.names=FALSE, quote=FALSE, na="0")