# R script for batch parsing and loading GHCN daily station files

# (*.dly) into a PostgreSQL database. The appropriate GHCN files are

# assumed to have been downloaded to the location specified by

# 'ghcndir', with the daily files themselves in a "ghcnd_all"

# subdirectory exactly as generated by unpacking "ghcnd_all.tar.gz";

# note that for the purposes of storage efficiency, we're not currently

# keeping these uncompressed files on disk, so this tarball would need

# to be unpacked again if for some reason this script needs to be

# re-run.

#

# The station data (names, locations) are now also imported as a

# separate 'stations' table in the database.

#

# As currently written, the script assumes that the 'ghcn' database

# already exists locally but has no tables, and that it can be accessed

# via 'ident' authentication (i.e., as the user executing the script,

# without password prompt).

#

# At the moment, PRCP, TMIN, and TMAX records are loaded.

#

# Jim Regetz

# NCEAS

# Created on 10-May-2012

require(RPostgreSQL)

#-------------#

# "constants" #

#-------------#

# location of ghcn daily data (on atlas)

ghcndir <- "/home/layers/data/climate/ghcn/v2.92-upd-2012052822"

# output file

logfile <- "~/ghcn-psql-load.log"

# name of target db

db.name <- "ghcn"

# variables to keep

VARS <- c("PRCP", "TMIN", "TMAX")

# pre-generate several command statements to avoid doing the work

# repeatedly when processing each of the 75k daily files

AWK.CMD <- paste(

    "awk -v FIELDWIDTHS='",

    paste(c(11, 4, 2, 4, rep(c(5,1,1,1), times=31)), collapse=" "),

    "' -v OFS=',' '{ $1=$1 \"\"; print }'", sep="")

SQL.UNION <- paste(

    sprintf(

        "(select id, year, month, %d as day, element,

            value%d as value, sflag%d as sflag,

            mflag%d as mflag, qflag%d as qflag from ghcn_wide)",

        1:31, 1:31, 1:31, 1:31, 1:31),

    collapse=" union all ")

#------------------#

# helper functions #

#------------------#

# Function invoking OS to preprocess data into CSV format with

# grep/awk/tr, then pipe directly into psql to invoke COPY statement for

# bulk load. Returns FALSE if filtering yields zero rows (or fails for

# some other reason), hence no data staged into the database.

loadAsCSV <- function(dly, patt=NULL) {

    tmpfile <- tempfile(tmpdir="/tmp")

    tr <- "tr -d ' '"

    if (is.null(patt)) {

        cmd <- paste(AWK.CMD, dly, "|", tr)

    } else {

        patt <- shQuote(paste(patt, collapse="\\|"))

        cmd <- paste("grep -e", patt, dly, "|", AWK.CMD, "|", tr)

    }

    cmd <- paste(cmd, tmpfile, sep=" > ")

    if (system(cmd)==0 & 0<file.info(tmpfile)$size) {

        dbGetQuery(db,

            sprintf("COPY ghcn_wide FROM '%s' WITH CSV", tmpfile))

        file.remove(tmpfile)

        TRUE

    } else {

        FALSE

    }

}

# Function to make postgresl move data from the wide-form table (days

# across columns) into the long form table (unique row for each

# day*element).

wideToLong <- function(db) {

   dbGetQuery(db, paste("insert into ghcn", SQL.UNION))

}

#-----------------#

# procedural code #

#-----------------#

# establish database connection

drv <- dbDriver("PostgreSQL")

db <- dbConnect(drv, dbname=db.name)

# create ghcn staging table

dbGetQuery(db, paste(

    "CREATE TABLE ghcn_wide (

        id char(11),

        year integer,

        month integer,

        element char(4), ",

        paste(sprintf("value%d integer, mflag%d char(1),

            qflag%d char(1), sflag%d char(1)", 1:31, 1:31, 1:31, 1:31),

            collapse=", "),

        ")", sep=""))

# create main ghcn table

dbGetQuery(db, paste(

    "CREATE TABLE ghcn (

        id char(11),

        year integer,

        month integer,

        day integer,

        element char(4),

        value integer,

        mflag char(1),

        qflag char(1),

        sflag char(1)

        )"))

# process and insert daily data

dailies <- list.files(file.path(ghcndir, "ghcnd_all"), pattern="*.dly")

for (file in dailies) {

    cat(date(), "\t", file=logfile, append=TRUE)

    if (loadAsCSV(file.path(ghcndir, "ghcnd_all", file), VARS)) {

        wideToLong(db)

        dbGetQuery(db, 'delete from ghcn_wide')

    }

    cat(file, "\n", file=logfile, append=TRUE)

}

# create ghcn stations table

dbGetQuery(db, paste(

    "CREATE TABLE stations (

        id char(11) primary key,

        latitude numeric,

        longitude numeric,

        elevation numeric,

        state char(2),

        name char(30),

        gsnflag char(3),

        hcnflag char(3),

        wmoid char(5)

        )"))

# process and insert station data

FMT.COLS <- strsplit(paste(c("A11", "F8", "F9", "F6", "A2", "A30", "A3",

    "A3", "A5"), collapse=",X1,"), ",")[[1]]

stations <- read.fortran(file.path(ghcndir, "ghcnd-stations.txt"),

    FMT.COLS, comment.char="", strip.white=TRUE, na="")

dbWriteTable(db, "stations", stations, row.names=FALSE, append=TRUE)

# add foreign key constraint to ghcn table

dbGetQuery(db,

    "ALTER TABLE ghcn

        ADD FOREIGN KEY (station)

        REFERENCES stations (id)")

dbDisconnect(db)

# alternative approach: load some number (>1) of station data files into

# the wide staging table before applying the wideToLong step; initial

# testing suggests that this approach is a little faster than

# one-at-a-time, but I don't know what the optimal number would be,

# particularly given that the amount of data per file is variable

#BATCH.SIZE <- 10

#counter <- 1

#for (file in dailies) {

#    loadAsCSV(file.path(ghcndir, "ghcnd_all", file), VARS)

#    if (counter %% BATCH.SIZE == 0) {

#        wideToLong(db)

#        dbGetQuery(db, 'delete from ghcn_wide')

#    }

#    counter <- counter + 1

#}