#!/opt/local/bin/perl
#
#  Name: get_kfact.pl -- Determine the KFACT for each STATSGO Mapunit
#
#  Purpose:
#
#	This script uses the Soil Erodibility Factor adjusted for the
#	effects of rock fragments (KFACT) and the rock fragment free
#	version (KFFACT) for the top layer of each STATSGO map unit
#	component and supporting data from the STATSGO COMP and LAYER
#	tables to determine average surface erodibility values for each
#	mapunit in the STATSGO MAPUNIT table.  The output is a decimal
#	fraction, with values between 0.00 and 1.00.
#
#  Usage:
#	perl get_kfact.pl workspace outfile rptfile
#
#  Options and Arguments:
#	workspace	Arc/Info Workspace location of INFO tables.
#	outfile		File to receive output table.
#	rptfile		File to receive report of processing statistics,
#			problems, etc.
#
#  External associations:
#	This script is part of a package for processing the STATSGO soil
#	properties database. In particular, it was written to access the
#	Layer and Component STATSGO tables to determine the soil
#	erodibility on a mapunit basis.  Output from this script is in
#	an ASCII format; script pack_kfact.pl may be used to convert it
#	to a format suitable for ingest into an Arc/Info INFO table
#	defined by script def_kfact_items.aml.
#
#  Algorithmic details:
#
#	The STATSGO LAYER table gives values of the soil erodibility
#	factor adjusted for the effects of rock fragments (KFACT) for
#	each layer of each component.  For most layers, a
#	rock-fragment-free erodibility factor, KFFACT, is also given;
#	this is the value normally used in the Universal Soil Loss
#	Equation to calculate soil loss by water.
#
#	This script extracts the values of KFACT and KFFACT for the
#	topmost layer for each component, and then averages the values
#	for all components of a given mapunit using the COMPPCT values
#	in the COMP table as weights.  If a non-zero value value of
#	KFACT is specified for a layer and KFFACT is missing or zero,
#	KFFACT is assumed to have the same value as KFACT.  Otherwise,
#	if KFACT is greater than KFFACT, a warning message is written.
#
#	If a mapunit has one or more components which are all water, and
#	also has other non-water components, the average KFACT and
#	KFFACT are computed for the non-water components only.  If the
#	mapunit is entirely water, KFACT and KFFACT are set to zero.  If
#	a component is all rock and there are no layers specified for
#	it, KFACT and KFFACT are also set to zero.
#
#  External files accessed:
#	<workspace>	Arc/Info Workspace location of INFO tables being 
#			accessed for processing.
#	<outfile>	WRITE to record KFACT/KFFACT for each mapunit.
#			Each record contains data for all components of
#			one mapunit.
#	<rptfile>	WRITE to record processing arguments, statistics
#			and other problems, etc., encountered during
#			processing.
#
#  Environment variables referenced:
#	PERLINC		Libraries to be searched by require() ahead of
#			those specified by -I and the default libraries.
#			This variable should be set if openinfo.pl is in
#			neither the directory from which get_kfact.pl is
#			invoked nor in the default libraries.
#
#  External subroutines included:
#	 openinfo.pl	Open an INFO table, returning unpack() template for
#			reading a specified list of items.
#
#  Internal subroutines defined:
#	help		Output help text.
#	add_comp	Subroutine to compute KFACT for one mapunit
#			component for the layers defined by successive
#			segments of the specified column lengths.
#	output_mu	Generate output values for each component of one
#			mapunit, and reinitialize for next component.
#	set_allrock	Subroutine to build single rock layer when no soil.
#	usage		Output terse usage text.
#
#  Internal variables:
#	COMP		File handle for reading COMP table listing.
#	LAYER		File handle for reading LAYER table listing.	 
#	RPT		File handle for output report file.
#	bufc		Array of values from COMP table.
#	bufl		Array of values from LAYER table.
#	ckfact		KFACT value for top layer of component.
#	ckffact		KFFACT value for top layer of component.
#	cl		Current component layer.
#	cmd		Name of this script; used for help text.
#	compname	Component name from COMP table.
#	comppct		Percent of component in map unit from COMP table.
#	comprec		Number of current COMP record.
#	curmuid		ID of mapunit currently being processed.
#	errcnt		Count of severe errors.
#	kfact1		KFACT value for last layer 1.
#	kffact1		KFFACT value for last layer 1.
#	kfsum		Sum of area-weighted KFACT for mapunit.
#	kffsum		Sum of area-weighted KFFACT for mapunit.
#	layer		Number of current output layer.
#	laynum		Value of LAYERNUM item.
#	layerrec	Number of current LAYER record
#	lmuid		Mapunit ID of current layer.
#	lseqnum		Component sequence number of current layer.
#	mubig		Flag set if MU comp has KFACT > KFFACT > 0.
#	muid		Mapunit ID.
#	muzero		Flag set if MU comp has KFACT > 0, KFFACT = 0.
#	nonwtotpct	Sum of percents for non water components in map
#			unit.
#	opt		Name of current option.
#	outfile		Name of file to receive relayered output.
#	rdeph		Maximum depth to bedrock for component.
#	seqnum		Component sequence number.
#	stats		Associative array of statistics for anomalies.
#	template_comp	Template for reading COMP table, returned by
#			openinfo.
#	template_lay	Template for reading LAYER table, returned by
#			openinfo.
#	time		Array of values returned by localtime().
#	totpct		Sum of percents for all components in map unit.
#	value		Value of one array element.
#
#  Script history:
#
#	1998-05		Initial version, R. A. White.
#	1999-10 	Changed treatment of KFFACT = 0, raw.
#
############   End of Prolog for get_kfact.pl   ############

#  Load external library modules, adding $PERLINC to search path

unshift (@INC, split (/ /,$ENV{"PERLINC"}));
require "openinfo.pl";

#  Define usage and help text

@cmd =	(split ('/', $0));
$cmd = pop (@cmd);

$usage = <<"End of Usage";

Usage:  $cmd workspace outfile rptfile
        $cmd -u|h
End of Usage

$help = <<"End of Help";

$cmd -- determine average surface soil erodabilities
$usage
Options and arguments are
     -h   	Write this help text to stdout and exit.
     -u   	Display terse usage message and exit.
   workspace	Arc/Info workspace location of INFO tables.
    outfile	File to receive output table.
    rptfile	File to receive report of processing statistics,
		problems, etc.
End of Help

#
#  Start of main routine
#

#  Get options

while (@ARGV && ($ARGV[0] =~ /^\w*-(\w+)($|\s+)/))
    {
    push (@opt, $1);
    shift;
    }	
(@opt) && ( (($opt[0] eq 'u') && (print $usage) && exit)
	 || (($opt[0] eq 'h') && (print $help) && exit) );

#  Preprocess arguments

(@ARGV != 3) && (print $usage) && exit;
$workspace = shift;
$outfile = shift;
$rptfile = shift;

#  Open	 output files

open (OUT, "> $outfile") || die "Error opening output file, $outfile,";
open (RPT, "> $rptfile") || die "Error opening report file, $rptfile,";

#  Create report header

@time = localtime (time);
$time[5] += 1900;
printf RPT "	%s, %02d/%02d/%02d %02d:%02d:%02d\n\n",
  "Summary of get_kfact.pl processing", $time[5], $time[4]+1,
  $time[3], $time[2], $time[1];

print RPT "Arguments are:\n";
printf RPT "    %-16s%s\n", "workspace", $workspace;
printf RPT "    %-16s%s\n", "outfile", $outfile;
printf RPT "    %-16s%s\n\n", "rptfile", $rptfile;

#  Call openinfo to open and unpack the INFO formatted COMP table

(($reclen_comp, $template_comp) = &openinfo ("COMP", $workspace,
    "COMP","MUID","SEQNUM","COMPNAME","COMPPCT", "ROCKDEPH"))
  || die "Openinfo returned error";

#  Call openinfo to open and unpack the INFO formatted LAYER table

(($reclen_layer, $template_lay) = &openinfo ("LAYER", $workspace,
   "LAYER", "MUID", "SEQNUM", "LAYERNUM", "KFACT", "KFFACT"))
 || die "Openinfo returned error";
   
#  Process components of each mapunit in turn

$curmuid = "";
$totpct = 0;
$nonwtotpct = 100;
$kfsum = 0;
$kfact1 = -1.0;
$kffsum = 0;
$kffact1 = -1.0;
$muzero = 0;
$mubig = 0;
$comprec = 0;
$layerrec = 0;
$errcnt = 0;

while (read (COMP, $bufc, $reclen_comp))
    {
    $comprec++;
    (($muid, $seqnum, $compname, $comppct, $rdeph)
	=  unpack ($template_comp, $bufc))
      || die "Bad COMP record $comprec";
    $curmuid = &output_mu if ($muid ne $curmuid);

 #  Check for missing LAYER entry for component

    if (defined $laynum)
	{
	print RPT "Bad LAYERNUM \"$laynum\" for layer 1,",
	  " $lmuid $lseqnum\n"
	  if (($laynum != 1) && ( ! eof (LAYER)));
	if (($lmuid ne $curmuid) || ($lseqnum != $seqnum))
	    {
	    &add_comp if (&set_allrock);
	    next;
	    }
	}
    while (read (LAYER, $bufl, $reclen_layer))
	{
	$layerrec++;
	(($lmuid, $lseqnum, $laynum, $kfact, $kffact)
	    = unpack ($template_lay, $bufl))
	  || die "Bad LAYER record $layerrec";

#  Get erodibilities for top layer of component

	if ($laynum == 1)
	    {
	    $ckfact = $kfact1;
	    $ckffact = $kffact1;
	    if ($kfact eq "")
		{
		print RPT "KFACT missing, $lmuid $lseqnum, layer",
		"$laynum\n";
		$stats{"nokfact"}++;
		$kfact = 0;
		}
	    if ($kffact eq "")
		{
		print RPT "KFFACT missing, $lmuid $lseqnum, layer",
		"$laynum\n";
		$stats{"nokffact"}++;
		$kffact = 0;
		}
	    $kfact1 = $kfact;
	    $kffact1 = $kffact;
	    }
	last if (($lseqnum != $seqnum) || ($lmuid ne $curmuid));
	}
    &add_comp;
    }

&output_mu;		# Output results for last mapunit

close COMP;
close LAYER;
close OUT;

#  Add statistics to processing report.

print RPT "\nProcessing statistics\n\n";
$value = 0 unless ($value = $stats{"totalcomp"});
print RPT "    Total components in input compfile:	$value\n"; 
print RPT "    Number of top layers with missing KFACT:   $value\n"
  if (($value = $stats{"nokfact"}) > 0);
print RPT "    Number of top layers with missing KFFACT:   $value\n"
  if (($value = $stats{"nokffact"}) > 0);
print RPT "    Number of top layers with KFACT > 0, KFFACT = 0:",
  "  $value in $stats{'mukffzero'} mapunits\n"
  if (($value = $stats{"kffzero"}) > 0);
print RPT "    Number of top layers with KFACT > KFFACT > 0:",
  "  $value in $stats{'mukfbig'} mapunits\n"
  if (($value = $stats{"kfbig"}) > 0);

close RPT;

exit 0;		# From Script get_kfact.pl

#
#  Subroutine to accumulate area-weighted values for current mapunit
#

sub add_comp
    {
    $stats{"totalcomp"}++;	# Total components
    return 1 if (($ckfact < 0.0) && ( ! &set_allrock));
    if ($ckfact > $ckffact)
	{
	if ($ckffact == 0.0)
	    {
	    $stats{"kffzero"}++;
	    $stats{"mukffzero"}++ unless ($muzero++);
	    $ckffact = $ckfact;
	    }
	else
	    {
	    $stats{"kfbig"}++;
	    $stats{"mukfbig"}++ unless ($mubig++);
	    print RPT "KFACT > KFFACT > 0.00 for layer 1,",
	      " $curmuid $seqnum\n";
	    }
	}
    $kfsum += $ckfact * $comppct;
    $kffsum += $ckffact * $comppct;
    $totpct += $comppct;

    return 1;

    }		# End of subroutine add_comp

#
#  Subroutine to output layer data for one mapunit
#

sub output_mu
    {
    return $muid if ($curmuid eq "");

#  Report problems and anomalies

    print RPT "Total COMPPCT for MUID $curmuid = $totpct\n"
      if ((($totpct < 99) || ($totpct > 101 ))
	&& (++$stats{"badtotpct"}));
#    print RPT "$curmuid has $muzero component(s) with KFACT > 0.00,",
#      " KFFACT = 0.00\n" if ($muzero);

#  Compute mean KFACT and KFFACT for entire mapunit

    if ($nonwtotpct > 0)
	{
	$ckfact =  $kfsum / $nonwtotpct;
	$ckffact =  $kffsum / $nonwtotpct;
	}
    else
	{
	print RPT "Non-water total COMPPCT for MUID $curmuid",
	" = $nonwtotpct -- setting KFACT and KFFACT to 0.00\n";
	$ckfact = 0.0;
	$ckffact = 0.0;
	}


#  Write record for this mapunit

    printf OUT "%-8s%4.2f%6.2f\n", $curmuid, $ckfact, $ckffact;

#  Reinitialize accumlators

    $totpct = 0;
    $nonwtotpct = 100;
    $kfsum = 0;
    $kffsum = 0;
    $muzero = 0;
    $mubig = 0;

    return $muid;
    }			#  End of subroutine output_mu

#
#  Subroutine to handle all-rock or all-water component
#

sub set_allrock
    {
    if (($compname =~ /^\s*WATER\s*$/) || ($curmuid =~ /^[A-Z]{2}W/))
	{
	if (($compname =~ /^\s*WATER\s*$/) && ($comppct > 0) &&
	  ($comppct < 100))
	    {
	    $ckfact = 0.0;
	    $ckffact = 0.0;
	    $nonwtotpct -= $comppct;
	    }
	else
	    {
	    $ckfact = 0.0;
	    $ckffact = 0.0;
	    $comppct = 100;
	    }
	}
    elsif ($rdeph == 0)
	{
	print RPT "No layers for $curmuid  $seqnum -- assuming",
	  " solid rock\n";
	$ckfact = 0.0;
	$ckffact = 0.0;
	}
    else
	{
	print RPT "MUID/SEQNUM mismatch -- ",
	  "compfile:  $curmuid $seqnum,	 ",
	  "layerfile:  $lmuid $lseqnum\n";
	die "Too many errors -- see report file $rptfile: "
	  if ($errcnt++ > 10);
	return 0;
	}
    return 1;
    }			#  End of subroutine set_allrock

#  End of script get_kfact.pl
