MapUnitAreas.sml

# MapUnitAreas.sml

## standalone script to tabulate areas of different map units occurring in contributing areas of catchments, using
## data produced by SampleCatchments.sml.

##### Inputs:
#	1) relocated sample point vector and sample catchment polygon vector objects created by SampleCatchments.sml
# 2) an overlapping geologic, soil, or other map with polygons attributed with a map unit identifier.

##### Outputs:
#	1) Catchment vector with UnitAreas table tabulating map unit areas and
#		area percentages for each catchment (summed over all contributing catchments);
#		records in this table are directly attached to the catchment polygons
#	2) Sample point vector with copy of the UnitAreas table with records directly
#		attached to the sample points
#	3) Intersect vector: logical intersect (VectorAND) of the catchment polygons and map
#		unit polygons, used to tabulate the map unit areas for each catchment.

# The UnitAreas tables include 
# -- CatchmentArea field that has area (sq Km) of all contributing catchment polygons 
# 		(computed for each catchment polygon by SampleCatchments.sml and read from PolyToPoly table)
# -- TotalArea_sqkm field that shows summed area of all map unit subpolygons found in the Intersect vector
#			for the current catchment polygon and its contributing catchment polygons; should equal CatchmentArea value
# -- TotalUnitPct field = 100 * TotalArea_sqkm / CatchmentArea; should equal 100% (check for correct processing)

##### Custom interface window allows user to:
#	-- choose input catchment and sample point vectors
# -- choose the table and field from the sample point vector that has the unique sample ID
# -- choose the table and field from the map unit vector that has the map unit identifiers
# -- optionally choose a particular target map unit to compile areas, otherwise they are
#			compiled for all map units.
# -- choose Project File for output vectors, which are automatically named: catchment and sample point
# 		vectors get input name + "areas"; Intersect vector named "Intersect".

### Script checks:  #######
# Sample point vector:
#	1) has point elements
# 2) has a table named PointToPoint created by SampleCatchements.sml
# 3) the sample  ID field selected by user is either string or numeric

# Catchment polygon vector:
# 1) has polygons
# 2) has a table named PolyToPoly created by SampleCatchments.sml
# 3) has a table with the same name as the sample ID table selected by user from the sample point vector

# Map unit vector 
# 1) has polygons
# 2) the map unit ID field selected by user is either string or numeric
# 3) is reprojected if necessary to match the CRS of the catchment vector to ensure proper result from intersect operation.

###  March 23, 2009
###  Randy Smith
###  MicroImages, Inc.

### Revised 2 July 2009: changed dialog from modal to modeless, enabling multiple runs with same data (must
			# choose new output Project File). 

###########  Global Variable Declarations  ###################

class RVC_VECTOR InSampPtVec, InCatchPolyVec, MapVec;
class RVC_VECTOR OutSampPtVec, OutCatchPolyVec, IntersectVec;
class SR_COORDREFSYS catchCRS, mapCRS;

class RVC_OBJITEM inSampPtObjItem, inCatchPolyObjItem, MapUnitObjItem;
string sampTblName$, sampFieldName$, mapUnitTblName$, mapUnitFieldName$;
class FILEPATH inputFilepath;
class STRING targetUnitStr = "";
class STRING sampIdFieldTypeStr, mapUnitIdFieldTypeStr;
class STRING outProjFilename$;

class DATABASE inCatchPolyDB;			# polygon database for input catchment vector
class DBTABLEINFO mapTbl, catchUnitAreaTbl, sampUnitAreaTbl;
class STRINGLIST mapUnitList;		# stringlist with single entry for each map unit

class STRINGLIST catchmentIDlist;		# stringlist with IDs for current catchment polygon and its upstream catchments 
string message$;

### Variables for the script dialog and controls
numeric err;		# for error handling
class STRING xml$;
class XMLDOC dlgdoc;
class XMLNODE dlgnode;

class GUI_DLG dlgwin;
class GUI_CTRL_PUSHBUTTON runBtn;
class GUI_CTRL_COMBOBOX targetUnitCombo;
class GUI_CTRL_EDIT_STRING ptVecText, polyVecText, mapVecText;
class GUI_CTRL_EDIT_STRING outProjFileText, targetUnitText, statusText;
class GUI_CTRL_LABEL targetUnitLabel;


####################################################
###########  User-defined procedures and functions  #############

####################################################
### function to return stringlist with single entry for each map unit identifier
func class STRINGLIST getMapUnitList (class DBTABLEINFO dbInfo, class STRING fieldName$)
	{
	local numeric i, j;							# counters
	local numeric unitInList;			# flag to indicate whether unit has been found in list
	local class STRING unit$;
	local class STRINGLIST list;

#	print("getMapUnitList function called.");
#	printf("Number of records in table = %d\n", dbInfo.NumRecords);
#	printf("mapUnitIdFieldTypeStr = %s\n", mapUnitIdFieldTypeStr);
#	printf("mapUnitFieldName$ = %s\n", mapUnitFieldName$);

	for i = 1 to dbInfo.NumRecords
		{
		# read map unit identifier field in record 
		if (mapUnitIdFieldTypeStr == "string")
			{
			unit$ = TableReadFieldStr(dbInfo, fieldName$, i);
			}
		else
			unit$ = NumToStr(TableReadFieldNum(dbInfo, fieldName$, i) );

		if (unit$ <> "")	# if string is not empty
			{
			# if stringlist is empty, add the unit identifier to the list
			if (list.GetNumItems() == 0) then
				list.AddToEnd(unit$);

			# otherwise loop through the stringlist to see if unit identifier is already in list
			else
				{
				unitInList = 0;
				for j = 1 to list.GetNumItems()
					{
					if (list[j - 1] == unit$) then unitInList = 1;
					}

				if (unitInList == 0) then
					list.AddToEnd(unit$);
				}		# end else
			}  # end if (unit$ <> "")
		}  # end for i = 1 to info.NumRecords

		list.Sort();		# sort list alphabetically

	return list;
	}  # end getMapUnitList()


####################################################
### function to return stringlist with IDs for current catchment and its upstream catchments
proc getUpstreamIDs
	(class STRING IDstr, 						# id string
	class STRING sampFld$,				# sample field name
	class DBTABLEINFO sampTbl, 	# sample table in the output point vector
	class DBTABLEINFO ptpTbl )		# PointToPoint table in the output point vector
	{
	local numeric j;		# counter
	local string fieldbase$ = "UpSample";		# base string for numbered UpSample fields
	local class STRINGLIST tempList;

	# add current sample ID to global ID List;
	catchmentIDlist.AddToEnd(IDstr);

	# get record number in the point sample table for the current sample ID;
	# this table is directly attached to the point elements
	local numeric recNum = TableKeyFieldLookup(sampTbl, sampFld$, IDstr);

	# get the element number of the sample point this record is attached to
	local array numeric ptNums[1];			# array of element numbers; should be just one
	TableGetRecordElementList(sampTbl, recNum, ptNums, "point");  

	# get the record in the PointToPoint table that is attached to this point
	local array numeric recNums[1];		# array of record numbers for attached records; should be just 1
	TableReadAttachment(ptpTbl, ptNums[1], recNums, "point");

	# get number of upstream samples immediately adjacent to the current basin from PointToPoint table
	local numeric numAdjUp = TableReadFieldNum(ptpTbl, "NumUpSamples", recNums[1]);
#	printf("Number of upstream samples = %d\n", numAdjUp);

	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
		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);
			local class STRING tempIDstr;

			# read ID from the current upstream field
			if (sampIdFieldTypeStr == "string") then 
				tempIDstr = TableReadFieldStr(ptpTbl, field$, recNums[1] );
			else
				tempIDstr = NumToStr( TableReadFieldNum(ptpTbl, field$, recNums[1] ) );

			tempList.AddToEnd(tempIDstr);		# add ID as string to temporary stringlist
			}

		# loop through temporary list of UpSample fields to get their adjacent upstream basin IDs
		# and call this function to check for basins further upstream
		for j =1 to numAdjUp
			{ 
			getUpstreamIDs( tempList[j-1], sampFld$, sampTbl, ptpTbl);
			}

		}	# end processing upstream samples / polygons

	} # end getUpstreamIDs( )



##########  User-defined procedures for the dialog window  #########

### procedure called by Catchment Poly Vector pushbutton
proc OnCatchBtn ()
	{
	local string message$;
	local numeric failed;
	local class RVC_GEOREFERENCE catchGeoref;

	err = DlgGetObject("Select USDBASIN vector object with catchment polygons:", "Vector", inCatchPolyObjItem, "ExistingOnly"); 
	if (err == 0)
		{
		InCatchPolyVec.Open(inCatchPolyObjItem, "Read");		# open the catchment vector to check database tables

		# check that selected vector has polygons
		if (InCatchPolyVec.$Info.NumPolys == 0)
			{
			failed = 1;
			message$ = "Selected vector does not have polygon elements.\n";
			}
		else
			{
			# check that catchment vector has a table named PolyToPoly from SampleCatchments script
			inCatchPolyDB =OpenVectorPolyDatabase(InCatchPolyVec); 

			if (!TableExists(inCatchPolyDB, "PolyToPoly") )
				{
				failed = 1;
				message$ = "Selected catchment vector does not have the required PolyToPoly table from SampleCatchments.sml.\n";
				}
			else
				{
				# get georeference and coordinate reference system from the catchment vector
				InCatchPolyVec.GetDefaultGeoref(catchGeoref);
				catchCRS = catchGeoref.GetCoordRefSys();

				polyVecText.SetValueStr(inCatchPolyObjItem.GetDisplayString() );
				dlgwin.GetCtrlByID("inSampBtn").SetEnabled(1);			# enable the Sample Point Vector button
				}
			}

		if (failed > 0)
			{
			message$ += "Please select a valid catchment vector object.";
			PopupMessage(message$);
			polyVecText.SetValueStr("");
			return;
			}
		}
	}

### procedure called by Sample Point Vector pushbutton
proc OnSamplePtBtn ()
	{
	local numeric failed;

	err = DlgGetObject("Select Result vector object with relocated sample points", "Vector", inSampPtObjItem, "ExistingOnly");
	if (err == 0)
		{
		InSampPtVec.Open(inSampPtObjItem, "Read");		# open the vector object and show info in dialog

		# check that selected vector has point elements.
		if (InSampPtVec.$Info.NumPoints == 0)		# no point elements; show error message
			{
			failed = 1;
			message$ = "The selected vector object has no point elements.\n";
			}
		else						# selected vector has point elements
			{
			local class DATABASE ptDB;		# get the point database
			ptDB = OpenVectorPointDatabase(InSampPtVec);

			# check that selected vector has PointToPoint table from SampleCatchments.sml
			if (!TableExists(ptDB, "PointToPoint") )
				{
				failed = 1;
				message$ = "The selected vector does not have the required PointToPoint table from SampleCatchments.sml.\n";
				}
			else
				{
				# set name of selected vector in dialog
				ptVecText.SetValueStr(inSampPtObjItem.GetDisplayString() );

				inputFilepath = InSampPtVec.$Info.Filename;

				# prompt to select table.field with sample IDs
				PopupMessage("Please select sample point table and field containing the unique sample ID.");
				PopupSelectTableField(ptDB, sampTblName$, sampFieldName$);

				if (sampTblName$ == "" or sampFieldName$ == "")
					{
					message$ = "The table and field for the unique sample ID have not been specified.\n";
					message$ += "Please reselect the Result vector and select the sample ID table and field.";
					PopupMessage(message$);
					ptVecText.SetValueStr("");
					return;
					}

				# check that the previously selected catchment polygon vector has a table with the same name
				if (!TableExists(inCatchPolyDB, sampTblName$) )
					{
					message$ = sprintf("The previously-selected catchment vector object does not have a table named %s.\n", sampTblName$);
					message$ += "Please select a different table for the sample id or select a catchment vector that includes this table.";
					PopupMessage(message$);
					}
				else		# catchment vector has a sample table with the same name
					{
					# get table info and field info to get field type for the sample ID field (numeric or string);
					local class DBTABLEINFO sampTblInfo;
					sampTblInfo = DatabaseGetTableInfo(ptDB, sampTblName$);
					local class DBFIELDINFO sampFieldInfo;
					sampFieldInfo =FieldGetInfoByName(sampTblInfo, sampFieldName$);
					sampIdFieldTypeStr =sampFieldInfo.Type;

					# check that selected ID field is string or numeric
					if (sampIdFieldTypeStr == "string" or sampIdFieldTypeStr == "numeric")
						{
						# show selected table and field names in dialog
						dlgwin.GetCtrlByID("sampTableLabel").SetEnabled(1);
						dlgwin.GetCtrlByID("sampTableText").SetValueStr(sampTblName$);
						dlgwin.GetCtrlByID("sampFieldLabel").SetEnabled(1);
						dlgwin.GetCtrlByID("sampFieldText").SetValueStr(sampFieldName$);

						dlgwin.GetCtrlByID("mapUnitBtn").SetEnabled(1);		# enable the Map Unit Vector button
						}
					else
						{
						message$ = "Sample ID field must be string or numeric.\n";
						message$ += "Please reselect the Result vector and select the sample ID table and field.";
						PopupMessage(message$);
						ptVecText.SetValueStr("");
						return;
						}
					} # end else polygon sample table exists
				} # end else point vector has PointToPoint table
			}

		if (failed > 0)
			{
			message$ += "Please select a valid sample point result vector object.";
			PopupMessage(message$);
			return;
			}
		}
	}

### procedure called by Map Unit Vector pushbutton
proc OnMapUnitBtn ()
	{
	local class RVC_GEOREFERENCE mapGeoref;
	local string message$;

	err = DlgGetObject("Select vector object with map unit polygons:", "Vector", MapUnitObjItem, "ExistingOnly");
	if (err == 0)
		{
		MapVec.Open(MapUnitObjItem, "Read");	# open the vector object and show info in dialog
		mapVecText.SetValueStr(MapUnitObjItem.GetDisplayString() );

		local class DATABASE mapUnitDB;	# prompt to select table.field with map unit IDs
		mapUnitDB = OpenVectorPolyDatabase(MapVec);
		PopupMessage("Please select the table and field containing the unique map unit ID.");
		PopupSelectTableField(mapUnitDB, 	mapUnitTblName$, mapUnitFieldName$);

		# check for Table Field dialog being canceled.
		if (mapUnitTblName$ == "" or mapUnitFieldName$ == "")
			{
			message$ = "The table and field for the map unit identifier have not been specified.\n";
			message$ += "Please reselect the Map Unit vector and select the map unit ID table and field.";
			PopupMessage(message$);
			mapVecText.SetValueStr("");
			return;
			}

		# get table info and field info to get field type for the map unit ID field (numeric or string);
		mapTbl = DatabaseGetTableInfo(mapUnitDB, mapUnitTblName$);
		local class DBFIELDINFO mapFieldInfo;
		mapFieldInfo =FieldGetInfoByName(mapTbl, mapUnitFieldName$);
		mapUnitIdFieldTypeStr =mapFieldInfo.Type;

		# check that selected ID field is string or numeric
		if (mapUnitIdFieldTypeStr == "string" or mapUnitIdFieldTypeStr == "numeric")
			{
			MapVec.GetDefaultGeoref(mapGeoref);
			mapCRS = mapGeoref.GetCoordRefSys();

			# show selected table and field names in dialog
			dlgwin.GetCtrlByID("mapTableLabel").SetEnabled(1);
			dlgwin.GetCtrlByID("mapTableText").SetValueStr(mapUnitTblName$);
			dlgwin.GetCtrlByID("mapFieldLabel").SetEnabled(1);
			dlgwin.GetCtrlByID("mapFieldText").SetValueStr(mapUnitFieldName$);

			# get list of map unit identifiers
			mapUnitList = getMapUnitList(	mapTbl, mapUnitFieldName$);

			dlgwin.GetCtrlByID("targetToggle").SetEnabled(1);		# enable toggle for single map unit
			dlgwin.GetCtrlByID("outProjFileBtn").SetEnabled(1);		# enable Output Project File button
			}
		else
			{
			message$ = "Map Unit ID field must be string or numeric.\n";
			message$ += "Please reselect the Map Unit vector and select the map unit ID table and field.";
			PopupMessage(message$);
			mapVecText.SetValueStr("");
			return;
			}
		}
	}

### procedure called by toggle button to tabulate areas for single map unit
proc OnTargetToggle ()
	{
	local numeric i;		# counter

	# get handles for map unit label and text field
	targetUnitLabel = dlgwin.GetCtrlByID("targetUnitLabel");
	targetUnitText = dlgwin.GetCtrlByID("targetUnitText");

	if (targetUnitStr == "")			# toggle was off, no target map unit string had been set
		{
		targetUnitLabel.SetEnabled(1);

		# get handle for the map unit combobox, populate the combobox, and enable it
		targetUnitCombo = dlgwin.GetCtrlByID("targetUnitCombo");

		for i = 1 to mapUnitList.GetNumItems()
			{
			targetUnitCombo.AddItem(mapUnitList[i-1], mapUnitList[i-1]);
			}

		targetUnitCombo.SetSelectedItemID( mapUnitList[0] );
		targetUnitCombo.SetEnabled(1);
		targetUnitStr = targetUnitCombo.GetSelectedItemID();
		}
	else													# toggle was on, target map unit string had been set, 
		{
		targetUnitStr = "";			# set blank target map unit string
		targetUnitLabel.SetEnabled(0);
		targetUnitText.SetEnabled(0);
		}
	}

### procedure called when target map unit is selected
proc OnMapUnitSelected ()
	{
	targetUnitStr = targetUnitCombo.GetSelectedItemID();
	}

### procedure called when Output Project File button is pressed
proc OnOutputProjFileBtn ()
	{
	local string message$;
	local string startDir$ = inputFilepath.GetPath();
	outProjFilename$ = GetOutputFileName(startDir$, "Select Project File to save output objects to:", "rvc");

	# check if user pressed "Cancel" in the selection dialog
	if (outProjFilename$ == startDir$)
		{
		message$ = "No output Project File name selected.\n";
		message$ += "Please select an output Project File.";
		PopupMessage(message$);
		return;
		}

	outProjFileText.SetValueStr(outProjFilename$);
	runBtn.SetEnabled(1);

	# Delete the file if it exists.  The user would have received a "Overwrite this file?" warning.
	if (fexists(outProjFilename$) )
		DeleteFile(outProjFilename$);

	CreateProjectFile(outProjFilename$, "Output objects from MapUnitAreas SML script");
	statusText.SetValueStr("");
	}


### procedure called when Close button is pressed
proc onClose ()
	{
	dlgwin.Close(0);
	Exit();
	}

#######################################
### procedure called when Run button is pressed
proc onRun ()
	{
	local numeric i, j, k;			# counters
	local class STRING areaFldStr, pctFldStr, idStr, mapUnitStr;
	local numeric id;
	local numeric catchmentArea, unitArea, sumUnitAreas;
	local numeric recNum, numAttached;
	local array numeric recordArray[1];
	local array numeric elemArray[1];

	clear();			# clear console

	############  CATCHMENT POLYGON VECTOR  ########
	# Make copy of catchment polygon vector in output Project File
	print("Making copy of catchment polygon vector.");
	statusText.SetValueStr("Making copy of catchment polygon vector.");
	CreateVector(OutCatchPolyVec,  outProjFilename$, "USDBASINareas", "Catchment polygons with tabulated map unit areas", "Polygonal");
	VectorCopyElements(InCatchPolyVec, OutCatchPolyVec);
	CloseVector(InCatchPolyVec);

	# open output catchment vector's polygon database
	class DATABASE catchPolyDB;
	catchPolyDB = OpenVectorPolyDatabase(OutCatchPolyVec);

	# get the PolyToPoly table info and info for the sample ID field 
	class DBTABLEINFO polyToPoly;
	polyToPoly = DatabaseGetTableInfo(catchPolyDB, sampTblName$);
	class DBFIELDINFO sampID_PtoP;
	sampID_PtoP = FieldGetInfoByName(polyToPoly, sampFieldName$); 

	# get the sample table info and info for the sample ID field to set as the key
	# for foreign key field in new unit area table
	class DBTABLEINFO outCatchSampTbl;
	outCatchSampTbl = DatabaseGetTableInfo(catchPolyDB, sampTblName$);
	class DBFIELDINFO sampID_ocSampTbl;

	##############  SAMPLE POINT VECTOR  ################
	# Make copy of relocated sample point vector in output Project File
	print("Making copy of sample point vector.");
	statusText.SetValueStr("Making copy of sample point vector.");
	CreateVector(OutSampPtVec, outProjFilename$, "SampPTareas", "Relocated sample points with tabulated map unit areas.", "Planar"); 
	VectorCopyElements(InSampPtVec, OutSampPtVec);
	CloseVector(InSampPtVec);

	# open the point database of the sample point vector and get info
	# for the PointToPoint table and sample table
	local class DATABASE sampPtDB;
	sampPtDB = OpenVectorPointDatabase(OutSampPtVec);
	local class DBTABLEINFO sampPtToPtTbl, sampPtTbl;
	sampPtToPtTbl = DatabaseGetTableInfo(sampPtDB, "PointToPoint");
	sampPtTbl = DatabaseGetTableInfo(sampPtDB, sampTblName$);


	#############   MAP UNIT VECTOR   ###################
	# compare CRS of catchment and map unit vectors and reproject
	# map unit vector to match if necessary.
	class RVC_VECTOR MapUnitVec;
	CreateTempVector(MapUnitVec, "Polygonal");

	if (mapCRS.Name <> catchCRS.Name)
		{
		print("Reprojecting the map unit vector to match the catchment vector.");
		statusText.SetValueStr("Reprojecting the map unit vector to match the catchment vector.");
		local class TRANS2D_MAPGEN trans;
		trans.InputCoordRefSys = mapCRS;
		trans.OutputCoordRefSys = catchCRS;

		VectorWarp(MapVec, MapUnitVec, trans, 10);
		}
	else
		VectorCopyElements(MapVec, MapUnitVec);

	#############  INTERSECT VECTOR  ###############
	# Do Logical Intersect ("AND") with map unit vector (source) and catchment polygon vector (operator)
	# Result is vector with basin polygons subdivided by map unit, with original basin and map unit records
	#    attached to the subpolygons.
	print("Doing logical Intersect of map unit vector and catchment polygon vector.");
	statusText.SetValueStr("Doing logical Intersect of map unit vector and catchment polygon vector.");
	CreateVector(IntersectVec, outProjFilename$, "Intersect", "Intersect of catchment polygons and map unit polygons.");
	IntersectVec = VectorAND(OutCatchPolyVec, MapUnitVec);

	# update the standard attributes for the polygons in the Intersect vector (to compute areas)
	VectorToolkitInit(IntersectVec);
	VectorUpdateStdAttributes(IntersectVec);

	# open the Intersect vector's polygon database and get handles for its sample, map unit,
	# and POLYSTATS tables
	class DATABASE intersectPolyDB;
	intersectPolyDB = OpenVectorPolyDatabase(IntersectVec); 
	class DBTABLEINFO intersectSampTbl, intersectMapUnitTbl, intersectPolystatsTbl;
	intersectSampTbl = DatabaseGetTableInfo(intersectPolyDB, sampTblName$);
	intersectMapUnitTbl = DatabaseGetTableInfo(intersectPolyDB, mapUnitTblName$);
	intersectPolystatsTbl = DatabaseGetTableInfo(intersectPolyDB, "POLYSTATS");


	############  CATCHMENT UNIT AREA TABLE  ################
	# make new unit area polygon table for output catchment vector
	print("Making unit area polygon table for the output catchment vector.");
	statusText.SetValueStr("Making unit area polygon table for the output catchment vector.");
	class DBTABLEINFO  catchUnitAreaTbl;
	catchUnitAreaTbl = TableCreate(catchPolyDB, "UnitAreas", "Map unit areas and cumulative areas for catchments");
	catchUnitAreaTbl.OneToOneAttachment = 1;

	# create field for sample ID 
	class DBFIELDINFO sampIDcatchFld;

	if (sampIdFieldTypeStr == "string")
		sampIDcatchFld = TableAddFieldString(catchUnitAreaTbl, sampFieldName$, 30);
	else
		sampIDcatchFld = TableAddFieldInteger(catchUnitAreaTbl, sampFieldName$, 30);

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

	# add fields for area and percent area for each map unit
	class DBFIELDINFO areaFld, pctFld;

	if (targetUnitStr == "")			# no target map unit selected; add fields for all map units
		{
		message$ = sprintf("Making fields for %d map units\n", mapUnitList.GetNumItems());
		print(message$);
		statusText.SetValueStr(message$);
		for i = 1 to mapUnitList.GetNumItems()
			{
			areaFldStr = mapUnitList[i-1] + "_Area"; 
			areaFld = TableAddFieldFloat(catchUnitAreaTbl, areaFldStr, 10, 4);
			areaFld.UnitType = "Area";
			areaFld.Units = "square kilometers";

			pctFldStr = mapUnitList[i-1] + "_Pct";
			pctFld = TableAddFieldFloat(catchUnitAreaTbl, pctFldStr, 10, 1);
			}
		}
	else					# user selected a target map unit; make fields for only this unit
		{
		print("Making fields for the target map unit.");
		statusText.SetValueStr("Making fields for the target map unit.");	
		areaFldStr = targetUnitStr + "_Area";
		areaFld = TableAddFieldFloat(catchUnitAreaTbl, areaFldStr, 10, 4);
		areaFld.UnitType = "Area";
		areaFld.Units = "square kilometers";

		pctFldStr = targetUnitStr + "_Pct";
		pctFld = TableAddFieldFloat(catchUnitAreaTbl, pctFldStr, 10, 1);
		}

	# add fields for area and percent area for all map units as another check
	class DBFIELDINFO totAreaFld, totPctFld;
	totAreaFld = TableAddFieldFloat(catchUnitAreaTbl, "TotalArea_sqkm", 10, 4);
	totAreaFld.UnitType ="Area";
	totAreaFld.Units = "square kilometers";

	totPctFld = TableAddFieldFloat(catchUnitAreaTbl, "TotalUnitPct", 10, 1);


	#############  ADD RECORDS WITH IDS  ######################
	# add a record to the unit area polygon table for each polygon in output catchment vector
	# and attach to the polygons
	print("Adding records to unit area polygon table.");
	statusText.SetValueStr("Adding records to unit area polygon table.");

	for i = 1 to polyToPoly.NumRecords
		{
		# get the polygon attached to the current PolyToPoly record
		numAttached = TableGetRecordElementList(polyToPoly, i, elemArray);

		# get the value from the TotalArea_sq_km field in the PolyToPoly table;
			# this is the total area of all contributing catchment polygons
		catchmentArea = TableReadFieldNum(polyToPoly, "TotalArea_sq_km", i);

		# get sample ID for the current record in the PolyToPoly table and
		# write new record to unit area polygon table with ID and total area values
		if (sampIdFieldTypeStr	=="string")
			{
			idStr = TableReadFieldStr(polyToPoly, sampFieldName$, i);
			recNum = TableNewRecord(catchUnitAreaTbl, idStr, catchmentArea);
			}
		else
			{
			id = TableReadFieldNum(polyToPoly,  sampFieldName$, i);
			recNum = TableNewRecord(catchUnitAreaTbl, id, catchmentArea);
			}

		recordArray[1] = recNum;
		TableWriteAttachment(catchUnitAreaTbl, elemArray[1], recordArray, 1, "polygon");
		}

	############  COMPUTE UNIT AREAS  ####################### 
	local numeric areaSqKm, unitPct, totalAreaPct;

	print("Computing map unit areas.");
	statusText.SetValueStr("Computing map unit areas.");

	if (targetUnitStr == "")		# if no target unit, computing areas for all map units
		numeric unitAreas[];				# set up hash to compile areas for map units (unit IDs used as keys to hash)

	else
		local numeric sumTargetAreas;		# simple numeric variable to accumulate areas for target unit

	# loop through records in unit area polygon table and sum total contributing
	# area for each basin (as check) and areas of each map unit in contributing
	# area, then compute percentage for each unit
	for i = 1 to catchUnitAreaTbl.NumRecords
		{
		# get sample ID for current record in the unit area polygon table as a string
		if (sampIdFieldTypeStr	=="string")
			idStr = TableReadFieldStr(catchUnitAreaTbl, sampFieldName$, i); 
		else
			idStr = NumToStr( TableReadFieldNum(catchUnitAreaTbl, sampFieldName$, i) );

		# call recursive function to get stringlist of all upstream sample/catchment IDs
		getUpstreamIDs(idStr, sampFieldName$, sampPtTbl, sampPtToPtTbl);

		#######################
		# loop through stringlist of contributing catchment IDs for the current catchment
		for j = 1 to catchmentIDlist.GetNumItems()
			{
			# get the record number for the current ID in the Intersect Vector sample table;
				# this record is attached to all unit (sub)polygons for the current catchment
			recNum = TableKeyFieldLookup(intersectSampTbl, sampFieldName$, catchmentIDlist[j-1]);

			# get element numbers of all (sub)polygons this sample number is attached to
			numAttached = TableGetRecordElementList(intersectSampTbl, recNum, elemArray);

			######################
			# loop through the intersect subpolygons to compile polygon areas and unit areas
			for k = 1 to numAttached
				{
				# get number of attached record in the Intersect Vector POLYSTATS table
				TableReadAttachment(intersectPolystatsTbl, 	elemArray[k], recordArray);

				# get the area from the POLYSTATS table and add to the running sum of subpolygon areas
				unitArea = TableReadFieldNum(intersectPolystatsTbl, "Area", recordArray[1] );
				sumUnitAreas +=unitArea;

				ResizeArrayClear(recordArray, 1);		# clear the record array

				# get number of attached record in the IntersectVector map unit table
				TableReadAttachment(intersectMapUnitTbl, elemArray[k], recordArray);

				# read map unit identifier field in record 
				if (mapUnitIdFieldTypeStr == "string")
					mapUnitStr = TableReadFieldStr(intersectMapUnitTbl, mapUnitFieldName$, recordArray[1]);
				else
					mapUnitStr = NumToStr(TableReadFieldNum(intersectMapUnitTbl, mapUnitFieldName$, recordArray[1]) );

				ResizeArrayClear(recordArray, 1);

				if (targetUnitStr == "")			# no target unit specified, add area for the current unit to its index in the unit hash
					unitAreas[mapUnitStr] += unitArea;

				else			# record area only if polygon matches target unit
					if (mapUnitStr == targetUnitStr)  then  sumTargetAreas += unitArea;

				} # end of loop through intersect subpolygons

				# clear the array of intersect polygon element numbers for the next pass
				ResizeArrayClear(elemArray, 0);

			}	# end of loop through stringlist of contributing catchment IDs

			catchmentIDlist.Clear();		# clear the stringlist of catchment IDs for the next pass

			# convert summed unit areas from square meters to square km
			sumUnitAreas = sumUnitAreas / 1000000;

		######################################
		# write computed results for the current catchment to appropriate fields in the 
		# output catchment vector's unit area table

		if 	(targetUnitStr == "")		# write area and percentage values for each map unit to the record
			{
			# loop through map unit stringlist to get summed unit areas for the current catchment
			for k = 1 to mapUnitList.GetNumItems()
				{
				# get summed area for the current map unit from the area hash and convert to square km
				areaSqKm = unitAreas[ mapUnitList[k-1] ] / 1000000;

				# write area if value is not null
				if (!IsNull(areaSqKm) )
					{
					unitPct = 100 * areaSqKm / sumUnitAreas;
					TableWriteField(catchUnitAreaTbl, i, sprintf("%s_Area",  	mapUnitList[k-1] ), areaSqKm);
					TableWriteField(catchUnitAreaTbl, i, sprintf("%s_Pct", mapUnitList[k-1] ), unitPct);
					}
				}
			} # end writing info for all map units

		else			# write area and percentage values for just the target map unit
			{
			areaSqKm = sumTargetAreas / 1000000;

			# write area if value is not null
			if ( !IsNull(areaSqKm) )
				{
				unitPct = 100 * areaSqKm / sumUnitAreas;

				TableWriteField(catchUnitAreaTbl, i, sprintf("%s_Area", targetUnitStr), areaSqKm);
				TableWriteField(catchUnitAreaTbl, i, sprintf("%s_Pct", targetUnitStr), unitPct);
				}
			}

		# compute total area percentage and write total unit area and percentage to current record
		totalAreaPct = 100 * sumUnitAreas / TableReadFieldNum(catchUnitAreaTbl, "CatchmentArea", i);

		TableWriteField(catchUnitAreaTbl, i, "TotalArea_sqkm", sumUnitAreas);
		TableWriteField(catchUnitAreaTbl, i, "TotalUnitPct", totalAreaPct);

		if (targetUnitStr == "")
			unitAreas.Clear();				# clear the unit area hash
		else
			sumTargetAreas = 0;			# reset the target area sum

		}	# end of computing unit areas

	print("Copying unit area table to the output sample point vector.");
	statusText.SetValueStr("Copying unit area table to the output sample point vector.");


	##################################################
	#######  UNIT AREA TABLE FOR OUTPUT SAMPLE POINT VECTOR#####

	# copy unit area table to the output sample point vector's point database
	TableCopy(catchPolyDB, catchUnitAreaTbl, sampPtDB);

	# get info for the new point table and for its sample ID field
	local class DBTABLEINFO ptUnitAreaTbl;
	ptUnitAreaTbl = DatabaseGetTableInfo(sampPtDB, "UnitAreas");
	local class DBFIELDINFO sampIDptUnitAreaTbl;
	sampIDptUnitAreaTbl = FieldGetInfoByName(ptUnitAreaTbl, sampFieldName$);

	# loop through records in point unit area table to attach records to points
	for i = 1 to ptUnitAreaTbl.NumRecords
		{ 
		# get the sample ID for the current record in the point unit area table and
		# find the corresponding record number in the sample point table
		if (sampIdFieldTypeStr	=="string")
			{
			idStr = TableReadFieldStr(ptUnitAreaTbl, sampFieldName$, i);
			recNum = TableKeyFieldLookup(sampPtTbl, sampFieldName$, idStr);
			}
		else
			{
			id = TableReadFieldNum(ptUnitAreaTbl,  sampFieldName$, i);
			recNum = TableKeyFieldLookup(sampPtTbl, sampFieldName$, id);
			}

		# get the point element number attached to this record in the sample point table
		numAttached = TableGetRecordElementList(sampPtTbl, recNum, elemArray);

		recordArray[1] = i;

		# attach the current record in the point unit area table to the corresponding point element	 
		TableWriteAttachment(ptUnitAreaTbl, elemArray[1], recordArray, 1);

		}

	print("Done.");
	statusText.SetValueStr("Done.");

	# disable Run button and clear output Project File; must select new output if want to run again
	runBtn.SetEnabled(0);
	outProjFileText.SetValueStr("");

	}		# end OnRun()


####################################################
###################   Main Program  ######################

# dialog specification
xml$ = '


	

		
			
				
				
			
			
				
				
			
			
				Sample Table Name 
				
				Unique ID Field 
				
			
		
		
			
				
				
			
			
				Map Unit Table Name 
				
				  Map Unit ID Field 
				
			
			
				
				  Target map unit 
				
			
		
		
			
			
		
		
			
			
			
		
	
';

# parse XML text for the dialog into memory; 
# return an error code (number < 0 ) if there are syntax errors
err = dlgdoc.Parse(xml$);

if (err < 0) 
	{
	PopupError(err); 	# Popup an error dialog. "Details" button shows syntax errors.
	Exit();
	}

# get the dialog element from the parsed XML document and
# show error message if the dialog element can't be found
dlgnode = dlgdoc.GetElementByID("mudlg");

if (dlgnode == 0) 
	{
	PopupMessage("Could not find dialog node in XML document");
	Exit();
	}

# Set the XML dialog element as the source for the GUI_DLG class instance
# we are using for the dialog window.
dlgwin.SetXMLNode(dlgnode);

dlgwin.CreateModeless();

## get handles for dialog controls
ptVecText = dlgwin.GetCtrlByID("ptVecText");
polyVecText = dlgwin.GetCtrlByID("polyVecText");
mapVecText = dlgwin.GetCtrlByID("mapVecText");
outProjFileText = dlgwin.GetCtrlByID("outProjFileText");

statusText = dlgwin.GetCtrlByID("statusText");
runBtn = dlgwin.GetCtrlByID("runBtn");

dlgwin.Open();
WaitForExit();