# 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$ = '<?xml version="1.0"?>
<!DOCTYPE root SYSTEM "smlforms.dtd">
<root>
<dialog id="mudlg" Title="Map Unit Areas" Buttons="">
<groupbox Name=" Select Input Vector Objects Created by SampleCatchments.sml: " ExtraBorder="4">
<pane Orientation="horizontal">
<pushbutton id="inCatchBtn" Name=" Catchment Poly Vector... " WidthGroup="1" OnPressed="OnCatchBtn()"/>
<edittext id="polyVecText" Width="40" ReadOnly="true"/>
</pane>
<pane Orientation="horizontal">
<pushbutton id="inSampBtn" Name=" Sample Point Vector... " WidthGroup="1" Enabled="false" OnPressed="OnSamplePtBtn ()"/>
<edittext id="ptVecText" Width="40" ReadOnly="true"/>
</pane>
<pane Orientation="horizontal">
<label id="sampTableLabel" Enabled="false">Sample Table Name </label>
<edittext id="sampTableText" Width="15" ReadOnly="true"/>
<label id="sampFieldLabel" Enabled="false">Unique ID Field </label>
<edittext id="sampFieldText" Width="15" ReadOnly="true"/>
</pane>
</groupbox>
<groupbox Name=" Select Map Unit Vector Object " ExtraBorder="4">
<pane Orientation="horizontal">
<pushbutton id="mapUnitBtn" Name=" Map Unit Vector... " WidthGroup="1" Enabled="false" OnPressed="OnMapUnitBtn()"/>
<edittext id="mapVecText" Width="40" ReadOnly="true"/>
</pane>
<pane Orientation="horizontal">
<label id="mapTableLabel" Enabled="false">Map Unit Table Name </label>
<edittext id="mapTableText" Width="15" ReadOnly="true"/>
<label id="mapFieldLabel" Enabled="false"> Map Unit ID Field </label>
<edittext id="mapFieldText" Width="15" ReadOnly="true"/>
</pane>
<pane Orientation="horizontal">
<togglebutton id="targetToggle" Name="Tabulate areas for single map unit " Enabled="false" OnChanged="OnTargetToggle()"/>
<label id="targetUnitLabel" Enabled="false"> Target map unit </label>
<combobox id="targetUnitCombo" Width="10" Enabled="false" OnSelection="OnMapUnitSelected()"/>
</pane>
</groupbox>
<pane Orientation="horizontal">
<pushbutton id="outProjFileBtn" Name=" Output Project File... " Enabled="false" OnPressed="OnOutputProjFileBtn()"/>
<edittext id="outProjFileText" Width="40" ReadOnly="true" />
</pane>
<pane Orientation="horizontal" HorizAlign="Right">
<edittext id="statusText" ReadOnly="true"/>
<pushbutton id="runBtn" Name=" Run " Enabled="false" OnPressed="onRun()"/>
<pushbutton Name=" Close " Enabled = "true" OnPressed="onClose()"/>
</pane>
</dialog>
</root>';
# 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();
Download