# this script for HWSD data only.

cut_HWSD <- function(fn.r, 
                     fn.buf,
                     fn.dbf,
                     fout.mask,
                     fa.soil, 
                     fa.geol,
                     crs.out, 
                     crs.gcs){
  fun.gdalcut(f.in = fn.r, f.mask = fn.buf, f.out = fout.mask, s_srs = crs.gcs, t_srs = crs.out)
  r = raster::raster(fout.mask)
  ur = sort(raster::unique(r))
  x = foreign::read.dbf(fn.dbf)
  
  cn =  c('SILT', 'CLAY', 'OC', 'BULK_DEN')
  cn.t = paste0('T_', cn )
  cn.s=paste0('S_', cn )
  idx = x[, 1] %in% ur
  y.soil = x[idx,  c('ID', cn.t)]
  y.geol = x[idx,   c('ID', cn.s)]
  fn1 = fn2 = fout.mask
  raster::extension(fn1) ='.Soil.csv'
  raster::extension(fn2) ='.Geol.csv'
  write.table(y.soil, file = fn1, quote = FALSE, row.names = FALSE, col.names = TRUE)
  write.table(y.geol, file = fn2, quote = FALSE, row.names = FALSE, col.names = TRUE)
  # foreign::write.dbf(y.geol, file = fn2)
  # write.table(y, quote = FALSE, row.names = FALSE, col.names = TRUE)
  return(list('Soil'=y.soil, 'Geol'=y.geol))
}

fun.soil.hwsd <-function(xfg){
  indir = xfg$dir$soil
  outdir = xfg$dir$predata
  fn.buf = xfg$pd.gcs$wbd.buf
  
  fn.r = file.path(indir, 'hwsd.bil')
  fn.dbf = file.path(indir, 'hwsd.dbf')
  
  fout.mask = file.path(xfg$dir$predata, 'hwsd.tif')
  fa.soil = file.path(xfg$dir$predata, 'hwsd.Soil.csv')
  fa.geol = file.path(xfg$dir$predata, 'hwsd.Geol.csv')
  
  tmp = cut_HWSD(fn.r, fn.buf, fn.dbf, fout.mask,  fa.soil, fa.geol,  crs.out = xfg$crs.pcs, crs.gcs = xfg$crs.gcs)
  
  xfg = c(xfg, 
          list(fn.soil = fout.mask,
               fn.geol = fout.mask,
               fa.soil = fa.soil,
               fa.geol = fa.geol) )
  dat.soil = fun.Soil_Geol(xfg=xfg, TOP = TRUE)
  dat.geol = fun.Soil_Geol(xfg=xfg, TOP = FALSE)
  return(xfg)
}