SSURGO to DrainMod via Rosetta

ssurgo2rosetta.R

## this script extracts soil component(s) data from SSURGO DB for use
## in DrainMod model. Rosetta should be used to derive VGM parameters
library(RODBC)
## library(RPostgreSQL)
library(plyr)

file <- 'ssurgo2rosetta.txt' # output file name to import into Rosetta

con <- odbcConnectAccess("C:/webdata/soils/soil_mn061/soildb_MN_2003.mdb")
## con <- odbcDriverConnect("Driver={Microsoft Access Driver (*.mdb, *.accdb)};DBQ=C:/webdata/soils/soil_mn061/soildb_MN_2003.mdb")
## con <- dbConnect(PostgreSQL(), dbname='lccmr')

## SSURGO has ks [um/s], we want [cm/hr] for DrainMod
## res <- dbGetQuery(con, '
## SELECT row_number() over (ORDER BY mukey, compname, hzdepb_r) as "order",
##   concat_ws(\'_\', musym, compname, hzname) AS "desc",
##   concat_ws(\'_\', musym, compname, comppct_r) as comp,
##   hzdepb_r,
##   sandtotal_r AS PSAND,
##   silttotal_r AS PSILT,
##   claytotal_r AS PCLAY,
##   dbthirdbar_r AS bd,
##   wthirdbar_r/100 AS theta33,
##   wfifteenbar_r/100 AS theta1500,
##   ksat_r * 0.36 AS Ks, awc_r AS theta, om_r AS ORGMAT
## FROM ssurgo.mapunit
##   JOIN ssurgo.component USING (mukey)
##   JOIN ssurgo.chorizon USING (cokey)
## WHERE (majcompflag!=\'No\' or majcompflag is null)
##   and mukey in (select distinct mukey from mike.ssurgo)
## ORDER BY mukey, compname, hzdepb_r
## ')             # extract preselected soils from SSURGO converted to PG

## dbDisconnect(con)

musym <- c("158B", "167b")        # substitute manually as appropriate

sql <- sprintf("
SELECT musym,
  musym & '_' & compname & '_' & hzname AS `desc`,
  hzdepb_r - hzdept_r AS HSUBLAY,
  sandtotal_r AS PSAND,
  silttotal_r AS PSILT,
  claytotal_r AS PCLAY,
  dbthirdbar_r AS bd,
  wthirdbar_r/100 AS theta33,
  wfifteenbar_r/100 AS theta1500,
  ksat_r AS Ks,
  awc_r AS theta,
  om_r AS ORGMAT
FROM (((mapunit
  INNER JOIN component ON mapunit.mukey = component.mukey)
  INNER JOIN chorizon ON component.cokey = chorizon.cokey))
WHERE musym IN ('%s')
  AND ((component.majcompflag='Yes' Or IsNull(component.majcompflag))=True)
ORDER BY musym,chorizon.hzdepb_r
", paste0(musym, collapse="','"))

res <- sqlQuery(con, sql)
odbcClose(con)

# sometimes several components are returned instead of major one
# in this case manual clean up is necessary like res <- res[5:7,]
head(res)

res$order <- unlist(by(res$musym, res$musym, seq))
rosetta <- subset(res, select=c(order, desc, PSAND, PSILT, PCLAY, bd, theta33, theta1500))

write("SSURGO export", file)
write.table(rosetta, file, append=TRUE, row.names=FALSE, quote=FALSE)

############################# P A U S E   H E R E ##############################
## 1. Run ROSETTA                                                              #
## 2. Create new DB named jd4-rosetta.mdb                                      #
## 3. Import ssurgo2rosetta.txt                                                #
## 4. Estimate parameters for all soils                                        #
## 5. Quit Rosetta and continue                                                #
################################################################################

con2 <- odbcConnectAccess("jd4-rosetta.mdb")

ros <- sqlQuery(con2, "
SELECT Description AS `desc`,
  Theta_r AS ORES, Theta_s AS OSAT,
  10^Alpha AS ALFA, 10^PredictedVG.Npar AS NPAR,
  10^Ko/24 AS KoSAT, 10^Ks/24 AS KSAT, Lpar AS LEXP,
  ALFA AS ALFAW, '0.0' AS H_ENPR
FROM ((PredictedVG
  INNER JOIN Rawdata ON PredictedVG.Code=Rawdata.Code)
  INNER JOIN UnsatMVG ON PredictedVG.Code=UnsatMVG.Code)
  INNER JOIN PredictedKs ON PredictedVG.Code=PredictedKs.Code
", as.is=TRUE)

odbcClose(con2)

out <- merge(res, ros, by="desc")

ret <- function(soil) {
    ## h <- seq(0, 15000, by=1000)
    h <- c(0, 50, 100, 200, 330, 500, 1000, 5000, 15000)
    data.frame(sat = with(soil,
               ORES + (OSAT - ORES) / (1 + abs(ALFA * h)^NPAR )^(1 - 1/NPAR)
               ), h = h)
}

## comp <- '102B_Clarion'
## x <- subset(out, comp==comp)

dir.create('soils')

dummy <- dlply(out, 'comp', function(x) {
    comp <- x[1, 'comp']
    fname <- sprintf("soils/%s.soi", comp)
    header <- sprintf("%d %s", nrow(x), comp)
    write(header, fname)
    write(sprintf("%d 0", nrow(x)), fname, append=TRUE)
    dummy <- dlply(x, 'hzdepb_r', function(soil) {
        ret.tbl <- ret(soil)
        header <- sprintf("%d %.2f %.1f ", nrow(ret.tbl), soil$KSAT, soil$hzdepb_r)
        write(header, fname, append=TRUE)
        ret.tbl <- format(ret.tbl, nsmall=1)
        write.table(ret.tbl, fname, append=TRUE, row.name=FALSE, col.names=FALSE, quote=FALSE)
    })
    footer <- sprintf("%g, %g, 1", 45, 30)
    write(footer, fname, append=TRUE)
})