GeolUnitArea.sml

### GeolUnitArea.sml

## 26 January 2006
## R. Smith
## MicroImages, Inc.

## standalone script to process sample points and basins
#$ created by the SampleCatchments script.

## Input: Results vector object containing relocated sample points and
## USDbasins vector containing basins (catchments) associated with 
## sample points.

## Output: copy of the USDbasins vector with added GeolUnitArea table 
## in the polygon database with records attached to polygon elements.
## Fields in this table record summed areas and percentages of each
## rock unit for each basin (catchment area).

## Assumption 1:  Polygons have attached PERCENTAGE table created by 
## the Polygon Properties process.  Multiple records in this table are 
## attached to each basin polygon, one for each geologic map
## unit polygon that overlaps the basin polygon.  Each such record 
## includes the area (in square meters) of the current basin polygon that
## is covered by that rock unit polygon.  More than one rock unit polygon 
## of the same type can overlap a given basin polygon.

## Assumption 2: Basins can be composite, consisting
## of more than one polygon element.

## Assumption 3: The sample for a particular basin represents it and
## all upstream basins in the same watershed.  The PointToPoint table
## records the number of upstream samples for each sample, and the
## IDs of those that are immediately adjacent to that sample's basin.
## Compiling a list of all upstream basins requires a recursive
## procedure to exhaustively trace these upstream-adjacent relationships.

## Assumption 4: Geologic map units in the area are indicated by symbols 
## Cu, Du, Mu, Ou, Qal, Tg, Ts, Tv, pCgr, and pCgs in the PERCENTAGE table. 

## Global variables
class VECTOR SampleVect, BasinVectIn, BasinVectOut;
class STRINGLIST idList;
class DATABASE polyDB;
class DBTABLEINFO ptTable, basinInfoTbl, unitAreaTbl, percentTbl;
class DBFIELDINFO sampID, sampIDunit, totArea, checkArea; 
class DBFIELDINFO cuAreaFld, cuPctFld, duAreaFld, duPctFld, muAreaFld, muPctFld, ouAreaFld;
class DBFIELDINFO ouPctFld, qalAreaFld, qalPctFld, tgAreaFld, tgPctFld, tsAreaFld, tsPctFld;
class DBFIELDINFO pCgrAreaFld, pCgrPctFld, pCgsAreaFld, pCgsPctFld, tvAreaFld, tvPctFld;
class DBFIELDINFO totUnitPctFld;
numeric i;


##########################################################################
## define recursive function to get upstream sample IDs from given sample
proc getUpstreamIDs ( string sampleID$ ) 
	{
	local numeric j;
	local string fieldbase$ = "UpSample";		# base string for numbered UpSample fields

	# add current sample ID to global ID List;
	idList.AddToEnd(sampleID$);

	# get record number in PointToPoint table for current sample ID
	local numeric recNum = TableKeyFieldLookup(ptTable, "REC_NO", sampleID$);

	# get number of upstream samples immediately adjacent to the current basin
	local numeric numAdjUp = TableReadFieldNum(ptTable, "NumUpSamples", recNum);


	if ( numAdjUp > 0 )	# if there are adjacent upstream polygons
		{
		# get list of sample IDs for adjacent upstream samples from UpSample[num] fields 
		# in the PolyToPoly table
		local class STRINGLIST tempList;
		for j = 1 to numAdjUp
			{
			# create string for field name holding the appropriate UpSample basin ID
			# and add to local stringlist
			local string field$ = fieldbase$ + NumToStr(j);
			tempList.AddToEnd( TableReadFieldStr(ptTable, field$, recNum) );
			}

		# loop through local list of UpSample fields to get their adjacent upstream basin IDs
		# and call procedure to check for basins further upstream
		for j = 1 to numAdjUp
			{
			getUpstreamIDs( tempList[j-1] );
			}
		}
	}

############################################################
################### Main process ###########################
############################################################

clear();		# clear Console window

###############################################################
## Get input vector and handles for existing tables needed for processing

GetInputVector(SampleVect);
GetInputVector(BasinVectIn);

GetOutputVector(BasinVectOut);
VectorCopyElements(BasinVectIn, BasinVectOut);
CloseVector(BasinVectIn);

# get handles for existing point and polygon tables needed for processing
ptTable = TableGetInfo(SampleVect.point.PointToPoint);
basinInfoTbl = TableGetInfo(BasinVectOut.poly.PolyToPoly);
sampID = FieldGetInfoByName(basinInfoTbl, "REC_NO");		# get DBFIELDINFO for REC_NO field
totArea = FieldGetInfoByName(basinInfoTbl, "TotalArea_sq_km");
percentTbl = TableGetInfo(BasinVectOut.poly.PERCENTAGE);

#########################################################
## create new polygon table
polyDB = OpenVectorPolyDatabase(BasinVectOut);
unitAreaTbl = TableCreate(polyDB, "UnitArea", "Geologic unit areas per basin"); 
unitAreaTbl.RelatedOnly = 1;

# create REC_NO field in GraniteArea table and make it foreign
# key field pointing to same field in PolyToPoly table
sampIDunit = TableAddFieldString(unitAreaTbl, "REC_NO", 12); 
sampIDunit.IsPrimaryKey = 0;
sampIDunit.IsKey = 1;
sampIDunit.Key = sampID;

# add TotalArea_sq_km field in new table using DBFIELDINFO from same field in PolyToPoly
TableAddField(unitAreaTbl, totArea);

# add CheckArea field to verify if search of upstream IDs succeeded
checkArea = TableAddFieldFloat(unitAreaTbl, "CheckArea", 10, 4);
checkArea.UnitType = "Area";
checkArea.Units = "square kilometers";

# add fields for area and % for each rock unit
cuAreaFld = TableAddFieldFloat(unitAreaTbl, "Cu_Area", 10, 4); 
cuPctFld = TableAddFieldFloat(unitAreaTbl, "Cu_Pct", 10, 2);
duAreaFld = TableAddFieldFloat(unitAreaTbl, "Du_Area", 10, 4); 
duPctFld = TableAddFieldFloat(unitAreaTbl, "Du_Pct", 10, 2);
muAreaFld = TableAddFieldFloat(unitAreaTbl, "Mu_Area", 10, 4); 
muPctFld = TableAddFieldFloat(unitAreaTbl, "Mu_Pct", 10, 2);
ouAreaFld = TableAddFieldFloat(unitAreaTbl, "Ou_Area", 10, 4); 
ouPctFld = TableAddFieldFloat(unitAreaTbl, "Ou_Pct", 10, 2);
qalAreaFld = TableAddFieldFloat(unitAreaTbl, "Tg_Area", 10, 4); 
qalPctFld = TableAddFieldFloat(unitAreaTbl, "Tg_Pct", 10, 2);
tgAreaFld = TableAddFieldFloat(unitAreaTbl, "Tg_Area", 10, 4); 
tgPctFld = TableAddFieldFloat(unitAreaTbl, "Tg_Pct", 10, 2);
tsAreaFld = TableAddFieldFloat(unitAreaTbl, "Ts_Area", 10, 4); 
tsPctFld = TableAddFieldFloat(unitAreaTbl, "Ts_Pct", 10, 2);
tvAreaFld = TableAddFieldFloat(unitAreaTbl, "Tv_Area", 10, 4); 
tvPctFld = TableAddFieldFloat(unitAreaTbl, "Tv_Pct", 10, 2);
pCgrAreaFld = TableAddFieldFloat(unitAreaTbl, "pCgr_Area", 10, 4); 
pCgrPctFld = TableAddFieldFloat(unitAreaTbl, "pCgr_Pct", 10, 2);
pCgsAreaFld = TableAddFieldFloat(unitAreaTbl, "pCgs_Area", 10, 4); 
pCgsPctFld = TableAddFieldFloat(unitAreaTbl, "pCgs_Pct", 10, 2);
totUnitPctFld = TableAddFieldFloat(unitAreaTbl, "Total_Unit_Pct", 15, 2);

###################################################################
## loop through records in PolyToPoly table, and make corresponding
## records in new GraniteArea table
for i = 1 to basinInfoTbl.NumRecords
	{
	local array numeric elemList[1];		# array holding list of elements attached to current record
	local numeric numElements;				# number of elements attached to current record
	local array numeric recNumbers[1];	# array of record numbers (1) to attach
	local numeric j;							# counter for element loop

	local string sampID$ = BasinVectOut.poly.PolyToPoly[@i].REC_NO$;
	local numeric totalAreaKm = BasinVectOut.poly.PolyToPoly[@i].TotalArea_sq_km;

	# Add new record to GraniteArea table with values for first two fields
	# from those in current record in PolyToPoly table
	TableNewRecord(unitAreaTbl, sampID$, totalAreaKm);
	}

## loop through records in UnitArea table and sum total contributing
## area for each basin (as check) and areas of each geologic unit in 
## contributing area, then compute percentage for each unit

for i = 1 to unitAreaTbl.NumRecords
	{
	local numeric j, k, m, n;
	local string sampID$ = BasinVectOut.poly.UnitArea[@i].REC_NO$;
	local numeric recNum;	# record number in PolyToPoly table
	local numeric sumBasinArea;	# running sum of polygon areas for basin

	# numeric variables for summed area and percent of each rock unit in map area 
	local numeric cuArea, duArea, muArea, ouArea, qalArea, tgArea, tsArea, tvArea, pCgrArea, pCgsArea;
	local numeric cuPct, duPct, muPct, ouPct, qalPct, tgPct, tsPct, tvPct, pCgrPct, pCgsPct;
	local numeric totUnitPct;

	# array to hold element numbers of polygons PolyToPoly record is attached to
	local array numeric polyList[1];
	local numeric numAttachedPolys;		# number of polygons attached

	# array to hold record numbers of PERCENTAGE table records attached to each polygon
	local array numeric pctRecordList[1]; 
	local numeric numPctRecords;

	# call recursive function to get list of all upstream sample/basin IDs
	getUpstreamIDs( sampID$ );

	# loop through list of contributing basins to sum geologic unit areas
	for j = 1 to idList.GetNumItems()
		{
		# get record number for basin's record in the PolyToPoly table
		recNum = TableKeyFieldLookup(basinInfoTbl, "REC_NO", idList[j-1]);

		# get list of polygon elements attached to that record
		numAttachedPolys = TableGetRecordElementList(basinInfoTbl, recNum, polyList);

		for k = 1 to numAttachedPolys		# loop through polygons
			{
			sumBasinArea += BasinVectOut.poly[ polyList[k] ].POLYSTATS.Area;

			# get list of record numbers of PERCENTAGE table records attached to polygon
			numPctRecords = TableReadAttachment(percentTbl, polyList[k], pctRecordList);

			for m = 1 to numPctRecords		# loop through attached records to compute sums
				{
				n = pctRecordList[m];

				if ( BasinVectOut.poly.PERCENTAGE[@n].FORMATION$ == "Cu") then
					cuArea += BasinVectOut.poly.PERCENTAGE[@n].Area;

				else if ( BasinVectOut.poly.PERCENTAGE[@n].FORMATION$ == "Du") then
					duArea += BasinVectOut.poly.PERCENTAGE[@n].Area;

				else if ( BasinVectOut.poly.PERCENTAGE[@n].FORMATION$ == "Mu") then
					muArea += BasinVectOut.poly.PERCENTAGE[@n].Area;

				else if ( BasinVectOut.poly.PERCENTAGE[@n].FORMATION$ == "Ou") then
					ouArea += BasinVectOut.poly.PERCENTAGE[@n].Area;

				else if ( BasinVectOut.poly.PERCENTAGE[@n].FORMATION$ == "Qal") then
					qalArea += BasinVectOut.poly.PERCENTAGE[@n].Area;

				else if ( BasinVectOut.poly.PERCENTAGE[@n].FORMATION$ == "Tg") then
					tgArea += BasinVectOut.poly.PERCENTAGE[@n].Area;

				else if ( BasinVectOut.poly.PERCENTAGE[@n].FORMATION$ == "Ts") then
					tsArea += BasinVectOut.poly.PERCENTAGE[@n].Area;

				else if ( BasinVectOut.poly.PERCENTAGE[@n].FORMATION$ == "Tv") then
					tvArea += BasinVectOut.poly.PERCENTAGE[@n].Area;

				else if ( BasinVectOut.poly.PERCENTAGE[@n].FORMATION$ == "pCgr") then
					pCgrArea += BasinVectOut.poly.PERCENTAGE[@n].Area;

				else if ( BasinVectOut.poly.PERCENTAGE[@n].FORMATION$ == "pCgs") then
					pCgsArea += BasinVectOut.poly.PERCENTAGE[@n].Area;
				}
			}
		}

	# use summed basin area and summed unit area to calculate percentage
	cuPct = 100 * cuArea / sumBasinArea;
	duPct = 100 * duArea / sumBasinArea;
	muPct = 100 * muArea / sumBasinArea;
	ouPct = 100 * ouArea / sumBasinArea;
	qalPct = 100 * qalArea / sumBasinArea;
	tgPct = 100 * tgArea / sumBasinArea;
	tsPct = 100 * tsArea / sumBasinArea;
	tvPct = 100 * tvArea / sumBasinArea;
	pCgrPct = 100 * pCgrArea / sumBasinArea;
	pCgsPct = 100 * pCgsArea / sumBasinArea;

	# sum unit percentages for each basin as check
	totUnitPct = round( (cuPct + duPct + muPct + ouPct + qalPct + tgPct + tsPct + tvPct + pCgrPct + pCgsPct ) );

	# convert summed areas from square meters to square km
	sumBasinArea = sumBasinArea / 1000000;
	cuArea = cuArea / 1000000;
	duArea = duArea / 1000000;
	muArea = muArea / 1000000;
	ouArea = ouArea / 1000000;
	qalArea = qalArea / 1000000;
	tgArea = tgArea / 1000000;
	tsArea = tsArea / 1000000;
	tvArea = tvArea / 1000000;
	pCgrArea = pCgrArea / 1000000;
	pCgsArea = pCgsArea / 1000000;

	# write computed results to appropriate fields in UnitArea table
	TableWriteField(unitAreaTbl, i, "CheckArea", sumBasinArea);

	TableWriteField(unitAreaTbl, i, "Cu_Area", cuArea);
	TableWriteField(unitAreaTbl, i, "Cu_Pct", cuPct);

	TableWriteField(unitAreaTbl, i, "Du_Area", duArea);
	TableWriteField(unitAreaTbl, i, "Du_Pct", duPct);

	TableWriteField(unitAreaTbl, i, "Mu_Area", muArea);
	TableWriteField(unitAreaTbl, i, "Mu_Pct", muPct);

	TableWriteField(unitAreaTbl, i, "Ou_Area", ouArea);
	TableWriteField(unitAreaTbl, i, "Ou_Pct", ouPct);

	TableWriteField(unitAreaTbl, i, "Qal_Area", qalArea);
	TableWriteField(unitAreaTbl, i, "Qal_Pct", qalPct);

	TableWriteField(unitAreaTbl, i, "Tg_Area", tgArea);
	TableWriteField(unitAreaTbl, i, "Tg_Pct", tgPct);

	TableWriteField(unitAreaTbl, i, "Ts_Area", tsArea);
	TableWriteField(unitAreaTbl, i, "Ts_Pct", tsPct);

	TableWriteField(unitAreaTbl, i, "Tv_Area", tvArea);
	TableWriteField(unitAreaTbl, i, "Tv_Pct", tvPct);

	TableWriteField(unitAreaTbl, i, "pCgr_Area", pCgrArea);
	TableWriteField(unitAreaTbl, i, "pCgr_Pct", pCgrPct);

	TableWriteField(unitAreaTbl, i, "pCgs_Area", pCgsArea);
	TableWriteField(unitAreaTbl, i, "pCgs_Pct", pCgsPct);

	TableWriteField(unitAreaTbl, i, "Total_Unit_Pct", totUnitPct);

	# clear all sum variables for next pass
	sumBasinArea = 0;
	cuArea = 0;		duArea = 0;		muArea = 0;		ouArea = 0;		qalArea = 0;
	tgArea = 0;		tsArea = 0;		tvArea = 0;		pCgrArea = 0;	pCgsArea = 0;	totUnitPct = 0;

	# clear global stringlist of contributing basin IDs
	idList.Clear();
	}