### 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();
}
Download