Dasymetric mapping (areal interpolation) script by Torrin Hultgren. Based on research by Jeremy Mennis. More information at http://astro.temple.edu/~jmennis/research/dasymetric/.

vectdasy.py

# ---------------------------------------------------------------------------
# VectDasy.py
# By Torrin Hultgren and Jeremy Mennis
# For use with ArcGIS 9.2 - 2007
# Usage: VectDasy <PopFeatureClass> <PopFieldName> <AreaFieldName>
# <PopKeyName> <AncFeatureClass> <AncCatName> <PresetTable> <PresetField>
# <OutFeatureClass> <SelectMethod> <Percent> <SampleMin> 
# ---------------------------------------------------------------------------

import string, arcgisscripting
gp = arcgisscripting.create()

# Load required toolboxes...
gp.AddToolbox("C:/Program Files/ArcGIS/ArcToolBox/Toolboxes/Data Management Tools.tbx")
gp.AddToolbox("C:/Program Files/ArcGIS/ArcToolbox/Toolboxes/Analysis Tools.tbx")

# Script arguments...
PopFeatureClass = gp.GetParameterAsText(0) # the name of the population feature class - this requires the full path
PopFieldName = gp.GetParameterAsText(1) # the field containing population counts
AreaFieldName = gp.GetParameterAsText(2) # the area field in the population layer 
PopKeyName = gp.GetParameterAsText(3) # the field in the population layer that uniquely identifies all of the population polygons
AncFeatureClass = gp.GetParameterAsText(4) # the name of the ancillary layer - this requires the full path
AncCatName = gp.GetParameterAsText(5) # the name of the field in the Ancillary layer containing category information
OutWorkspace = gp.GetParameterAsText(6) # Output path
OutFeatureClass = gp.GetParameterAsText(7) # Output dataset name
PresetTable = gp.GetParameterAsText(11) # Table with ancillary categorys and preset values
PresetField = gp.GetParameterAsText(12) # Field in preset table with preset values - field with class values assumed to be the same as AncCatName
SelectMethod = gp.GetParameterAsText(9) # HAVE_THEIR_CENTER_IN=centroid method, CONTAINED_BY=wholly contained method, and PERCENT_AREA=percent area method
Percent = gp.GetParameterAsText(10) # optional parameter - percent value for percent area method - default = 0.95
SampleMin = gp.GetParameterAsText(8) # Minimum number of source units to ensure a representative sample - default = 2

## To run the script as a standalone, rather than from ArcToolbox,
## comment out the lines above, uncomment the following lines and set the arguments manually

##PopFeatureClass = "D:\\Dasymetric\\SampleData\\tracts_Clip.shp"
##PopFieldName = "TOTPOP"
##AreaFieldName = "Area"
##PopKeyName = "Unique"
##AncFeatureClass = "D:\\Dasymetric\\SampleData\\lc_2k.shp"
##AncCatName = "LC"
##PresetTable = "#"
##PresetField = "#"
##OutWorkspace = "D:/Dasymetric/SampleData/"
##OutFeatureClass = "test_out_centroid"
##SelectMethod = "HAVE_THEIR_CENTER_IN"
##Percent = "#"
##SampleMin = "5"

try:
    gp.Workspace = OutWorkspace
    if gp.Describe(OutWorkspace).WorkspaceType == 'FileSystem':
        TableSuffix = '.dbf'
    else:
        TableSuffix = ''
    # Make layers from the feature classes for selection purposes
    def GetLayerName(FeatureClass):
        return (FeatureClass[(FeatureClass.rfind('\\')+1):]).split(".")[0]
    PopLayerName = GetLayerName(PopFeatureClass)
    PopLayer = gp.MakeFeatureLayer(PopFeatureClass,PopLayerName + "_lyr")
    AncLayerName = GetLayerName(AncFeatureClass)
    AncLayer = gp.MakeFeatureLayer(AncFeatureClass,AncLayerName + "_lyr")
    def NameCheck(Name):
        j,OkName = 1,""
        while not OkName:
            tList = gp.ListTables(Name)
            fList = gp.ListFeatureClasses(Name)
            if tList.Next() or fList.Next():
                Name = Name[:-2] + "_" + str(j)
            else:
                OkName = Name
            j = j + 1
        return OkName
    # AncDissolveName = NameCheck("AncDissolve")
    # To ensure that no population is lost, make sure that only source units contained by the extent of the
    # - ancillary dataset are used in the procedure
    # Dissolve_management (in_features, out_feature_class, dissolve_field, statistics_fields, multi_part)
    # AncDissolve = gp.Dissolve_management(AncFeatureClass, AncDissolveName)
    # AncDLayer = gp.MakeFeatureLayer(AncDissolve, AncDissolveName + "_lyr")
    # SelectLayerByLocation_management (in_layer, overlap_type, select_features, search_distance, selection_type)
    # PopLayer = gp.SelectLayerByLocation(PopLayer, "CONTAINED_BY", AncDLayer, 0, "NEW_SELECTION")
    # gp.Delete_management(AncDLayer)
    # gp.Delete_management(AncDissolve)
    
    # Process: Intersect...
    # Intersect_analysis (in_features, out_feature_class, join_attributes, cluster_tolerance, output_type)
    print "Performing Intersection..."
    gp.AddMessage("Performing Intersection...")
    OutLayerName = NameCheck(GetLayerName(OutFeatureClass))
    OutFeatureClass = gp.Intersect_analysis(PopLayer + ";" + AncLayer, OutLayerName, "ALL", "", "INPUT")
    # Create a population working table so that no changes to the source data are necessary
    print "Creating a working table for calculations..."
    gp.AddMessage("Creating a working table for calculations...")
    PopWorkingTableName = NameCheck("PopWorkingTable")
    # Frequency_analysis (in_table, out_table, frequency_fields, summary_fields)
    PopWorkingTable = gp.Frequency_analysis(PopLayer, PopWorkingTableName + TableSuffix, PopKeyName + ";" + PopFieldName + ";" + AreaFieldName)
    # Look up Field Properties
    FieldTypeLookup = {"Integer":"LONG","String":"TEXT","Single":"FLOAT","Double":"DOUBLE","SmallInteger":"SHORT","Date":"DATE"}
    Fields = gp.ListFields(AncFeatureClass, AncCatName)
    Field = Fields.Next()
    AncCatFType = Field.Type
    FieldProps = [FieldTypeLookup[AncCatFType],Field.Length]
    print "Adding fields to working table and calculating density..."
    gp.AddMessage("Adding fields to working table and calculating density...")
    # AddField_management (in_table, field_name, field_type, field_precision, field_scale, field_length, field_alias, field_is_nullable, field_is_required, field_domain)
    for Field in ["DENSITY","POPAREA","PDENSITY"]:
        gp.AddField_management(PopWorkingTable, Field, "DOUBLE")
    gp.AddField_management(PopWorkingTable, "REPCAT", FieldProps[0], "", "", FieldProps[1])
    # MakeTableView_management (in_table, out_view)
    PopWTableView = gp.MakeTableView_management(PopWorkingTable, "PopWTableView")
    gp.CalculateField_management(PopWTableView, "DENSITY", "[" + PopFieldName + "] / [" + AreaFieldName + "]")
    # Get appropriate wrappers for fieldnames for queries - " " or [  ]
    Wrp = {}
    pFC, aFC, oFC = PopFeatureClass, AncFeatureClass, OutFeatureClass
    for fc in [pFC, aFC, oFC]:
        description = gp.Describe(fc)
        if description.DataType == 'FeatureClass':
            Wrp[fc], fL = ['[',']'], 255
        else:
            Wrp[fc], fL = ['"','"'], 10 # If the output feature class is a shapefile, fieldnames will be truncated.
    
    # Process: Add Field...
    # Create a list of the fields to add
    print "Adding fields to the output table..."
    gp.AddMessage("Adding fields to the output table...")
    for Field in ["POPAREA","POPEST","REMAREA","DENSEST","TOTALFRACT","NEWPOP","NEWDENSITY"]:
        gp.AddField_management(OutFeatureClass, Field, "DOUBLE")

    # Process: gather unique IDs
    # Gather unique ancillary IDs into a table
    print "Collecting unique ancillary categories and building indices..."
    gp.AddMessage("Collecting unique ancillary categories and building indices..")
    AncCategoryTableName = NameCheck("AncCategoryTable")
    AncCatTable = gp.Frequency_analysis(OutFeatureClass, AncCategoryTableName + TableSuffix, AncCatName[:fL])
    # Gather the IDs from the table to a python list object - GetValues is a function defined below
    def GetValues(Table, Field):
        InList = []
        if Table and Table <> "#":
            FieldList = gp.ListFields(Table, Field)
            FieldObj = FieldList.Next()
            Rows = gp.SearchCursor(Table)
            Row = Rows.Next()
            if FieldObj.Type == "String":
                while Row:
                    InList.append(Row.GetValue(Field))
                    Row = Rows.Next()
            else:
                while Row:
                    InList.append(str(int(Row.GetValue(Field))))
                    Row = Rows.Next()
            Rows, Row, FieldList, FieldObj = None, None, None, None
        return InList
    InAncCatList = GetValues(AncCatTable, AncCatName)
    OutAncCatList = InAncCatList
    if AncCatFType == "String":
        InAncCatList = ["'" + AncCat + "'" for AncCat in InAncCatList]
        OutAncCatList = ['"' + AncCat + '"' for AncCat in OutAncCatList]
    gp.Delete_management(AncCatTable)
    PopFeatureClass, PopLayerName, AncLayerName, AncDissolveName, FieldTypeLookup = None, None, None, None, None
    Fields, description, fc, Field, AncCategoryTableName = None, None, None, None, None

    # Create an index on both source unit ID and ancillary ID to speed processing
    # AddIndex_management (in_table, fields, index_name, unique, ascending)
    gp.AddIndex_management(OutFeatureClass, AncCatName[:fL], AncCatName + "_atx")
    gp.AddIndex_management(OutFeatureClass, PopKeyName[:fL], PopKeyName + "_atx") 
    
    # Calculate area of all the intersected polygons
    # - note that the units are determined by the linear unit of the dataset's spatial reference
    # Calculate areas function adapted from ESRI Calculate Areas Python script
    if Wrp[oFC][0] == '"':
        print "Calculating area of intersected polygons..."
        gp.AddMessage("Calculating area of intersected polygons...")
        gp.AddField_management(OutFeatureClass, "NEWAREA", "DOUBLE")
        pFields = gp.ListFields(OutFeatureClass)
        pField = pFields.Next()
        while pField:
            sName = (pField.Name).upper()
            sType = pField.Type
            if sType == "Geometry":
                sShapeField = sName
                break
            pField = pFields.Next()
        pFields = None
        pRows = gp.UpdateCursor(OutFeatureClass)
        pRow = pRows.Next()
        while pRow <> None:
            pGeometry = pRow.GetValue(sShapeField)
            sArea = pGeometry.Area 
            pRow.SetValue("NEWAREA",sArea)
            pRows.UpdateRow(pRow)
            pRow = pRows.Next()
        pRows = None
        NewArea = "NEWAREA"
        del pField, sName, sType, sShapeField, pRow, pGeometry, sArea
    elif Wrp[oFC][0] == '[':
        NewArea = "Shape_Area"
    
    # Create dictionary object of ancillary classes and their preset densities, as well as list of classes preset to zero
    UninhabList, UnsampledList = [], []
    InPresetCatList = GetValues(PresetTable, AncCatName)
    OutPresetCatList = InPresetCatList
    if AncCatFType == "String":
        InPresetCatList = ["'" + AncCat + "'" for AncCat in InPresetCatList]
        OutPresetCatList = ['"' + AncCat + '"' for AncCat in OutPresetCatList] 
    PresetValList = GetValues(PresetTable, PresetField)
    i = 0
    for PresetCat in InPresetCatList:
        if float(PresetValList[i]) == 0:
            UninhabList.append(PresetCat)
        i = i + 1
    i = None

    # Make layer from the feature class for selection purposes
    OutLayer = gp.MakeFeatureLayer(OutFeatureClass,OutLayerName + "_lyr")
    
    # Create a summary table of the "populated" area of each population source unit
    print "Creating summary table with the populated area of each source unit..."
    gp.AddMessage("Creating summary table with the populated area of each source unit...")
    # SelectLayerByAttribute_management (in_layer_or_view, selection_type, where_clause)
    if UninhabList:
        OutLayer = gp.SelectLayerByAttribute_management(OutLayer, "NEW_SELECTION", Wrp[oFC][0] + AncCatName[:fL] + Wrp[oFC][1] + " NOT IN (" + ", ".join(UninhabList) + ")")
    else:
        OutLayer = gp.SelectLayerByAttribute_management(OutLayer, "CLEAR_SELECTION")
    InhabAreaTableName = NameCheck("InhabAreaTable")
    # Statistics_analysis (in_table, out_table, statistics_fields, case_field)
    InhabAreaTable = gp.Frequency_analysis(OutLayer, InhabAreaTableName + TableSuffix, PopKeyName[:fL], NewArea)
    # AddJoin_management (in_layer_or_view, in_field, join_table, join_field, join_type)
    gp.AddJoin_management(OutLayer, PopKeyName[:fL], InhabAreaTable, PopKeyName[:fL], "KEEP_COMMON")
    # CalculateField_management (in_table, field, expression)
    gp.CalculateField_management(OutLayer, OutLayerName + ".POPAREA", "[" + InhabAreaTableName + "." + NewArea + "]")
    # RemoveJoin_management (in_layer_or_view, join_name)
    gp.RemoveJoin_management(OutLayer, InhabAreaTableName)
    # Do the same for the POPAREA field in the PopWorkingTable
    gp.AddJoin_management(PopWTableView, PopKeyName[:fL], InhabAreaTable, PopKeyName[:fL], "KEEP_COMMON")
    gp.CalculateField_management(PopWTableView, PopWorkingTableName + ".POPAREA", "[" + InhabAreaTableName + "." + NewArea + "]")
    gp.RemoveJoin_management(PopWTableView, InhabAreaTableName)
    gp.Delete_management(InhabAreaTable)
    del InhabAreaTableName

    # The default value for a number field is zero in a shapefile, but is null in a geodatabase
    # These nulls cause problems in summary operations, so define a function to set them to zero.
    def RemoveNulls(Table,Field,Wrp,OutFeatureClass): 
        Table = gp.SelectLayerByAttribute_management(Table, "NEW_SELECTION", Wrp[oFC][0] + Field + Wrp[oFC][1] + " IS NULL")
        if gp.GetCount_management(Table):
            gp.CalculateField_management(Table, Field, "0")
    RemoveNulls(OutLayer,"POPAREA",Wrp,OutFeatureClass)
    RemoveNulls(PopWTableView,"POPAREA",Wrp,OutFeatureClass)
    # Select all non-zero population units (dividing by zero to get density is bad)
    POPAREAWhereClause = Wrp[oFC][0] + "POPAREA" + Wrp[oFC][1] + " > 0"
    PopWTableView = gp.SelectLayerByAttribute_management(PopWTableView, "NEW_SELECTION", POPAREAWhereClause)
    gp.CalculateField_management(PopWTableView, "PDENSITY", "[" + PopFieldName[:fL] + "] / [POPAREA]")
    
    # Begin selection process...
    print "Selecting representative source units..."
    gp.AddMessage("Selecting representative source units...")
    # The goal is to calculate "representative" population density values for each ancillary class
    # - by selecting all population units that "represent" that ancillary class and summing population and area.
    # - There are three methods for selecting "representative" units: Centroid, Wholly Contained, and Percent Area
    if SelectMethod == "PERCENT_AREA":
        # Percent area method
        # ***We are assuming the same selection set here - only inhabited classes
        PercentAreaTableName = NameCheck("PercentAreaTable")
        Outlayer = gp.SelectLayerByAttribute_management(OutLayer, "NEW_SELECTION", POPAREAWhereClause)
        # Frequency_analysis (in_table, out_table, frequency_fields, summary_fields)
        PercentAreaTable = gp.Frequency_analysis(OutLayer, PercentAreaTableName + TableSuffix, PopKeyName[:fL] + ";" + AncCatName[:fL] + ";POPAREA", NewArea)
        gp.AddField_management(PercentAreaTable, "Percent", "DOUBLE")
        # CalculateField_management (in_table, field, expression)
        gp.CalculateField_management(PercentAreaTable, "Percent", "[" + NewArea + "] / [POPAREA]")
        for InAncCat, OutAncCat in zip(InAncCatList, OutAncCatList):
            WhereClause = Wrp[oFC][0] + "Percent" + Wrp[oFC][1] + " >= " + Percent
            WhereClause = WhereClause + " AND " + Wrp[oFC][0] + AncCatName[:fL] + Wrp[oFC][1] + " = " + InAncCat
            PopSelSetName = NameCheck("PopSelSet")
            # TableSelect_analysis (in_table, out_table, where_clause)
            PopSelSet = gp.TableSelect_analysis(PercentAreaTable, PopSelSetName, WhereClause)
            # GetCount_management (in_rows) 
            Count = gp.GetCount_management(PopSelSet)
            print InAncCat + " has " + str(Count) + " representative source units."
            gp.AddMessage(InAncCat + " has " + str(Count) + " representative source units.")
            # If the selection set is not empty...
            if Count >= long(SampleMin):
                # AddJoin_management (in_layer_or_view, in_field, join_table, join_field, join_type)
                gp.AddJoin_management(PopWTableView, PopKeyName[:fL], PopSelSet, PopKeyName[:fL], "KEEP_COMMON")
                # Designate these source units as representative of the current ancillary category by putting the category value in the REPCAT field
                gp.CalculateField_management(PopWTableView, PopWorkingTableName + ".REPCAT", OutAncCat)
                # RemoveJoin_management (in_layer_or_view, join_name)
                gp.RemoveJoin_management(PopWTableView, PopSelSetName)                
            # Flag this class if it was insufficiently sampled and not preset
            elif InAncCat not in InPresetCatList:
                UnsampledList.append(InAncCat)
            # End loop for this ancillary class...
            gp.Delete_management(PopSelSet)
        gp.Delete_management(PercentAreaTable)
 
    else:  # For Centroid or Wholly Contained Methods...
        for InAncCat, OutAncCat in zip(InAncCatList, OutAncCatList):
            WhereClause = Wrp[aFC][0] + AncCatName[:fL] + Wrp[aFC][1] + " = " + InAncCat
            AncLayer = gp.SelectLayerByAttribute_management(AncLayer, "NEW_SELECTION", WhereClause)
            # Select ancillary class polygons
            if SelectMethod == "HAVE_THEIR_CENTER_IN": # Centroid Selection Method
                PopLayer = gp.SelectLayerByLocation_management(PopLayer, SelectMethod, AncLayer, 0, "NEW_SELECTION") 
            elif SelectMethod == "CONTAINED_BY": # Wholly Contained Method
                PopLayer = gp.SelectLayerByLocation_management(PopLayer, SelectMethod, AncLayer, 0, "NEW_SELECTION")
            # GetCount_management (in_rows) 
            Count = gp.GetCount_management(PopLayer)
            print InAncCat + " has " + str(Count) + " representative source units."
            gp.AddMessage(InAncCat + " has " + str(Count) + " representative source units.")
            # If the selection set is not empty...
            if Count >= long(SampleMin):
                PopSelSetName = NameCheck("PopSelSet")
                # Create a new, temporary table containing the IDs of the selected source units
                PopSelSet = gp.Frequency_analysis(PopLayer, PopSelSetName, PopKeyName)
                # Join this table to the working table
                gp.AddJoin_management(PopWTableView, PopKeyName[:fL], PopSelSet, PopKeyName[:fL], "KEEP_COMMON")
                # Designate these source units as representative of the current ancillary category by putting the category value in the REPCAT field
                gp.CalculateField_management(PopWTableView, "REPCAT", OutAncCat)
                # Remove the join and clean up the temporary table
                gp.RemoveJoin_management(PopWTableView, PopSelSetName)
                gp.Delete_management(PopSelSet)
                PopSelSetName = None
            # Flag this class if it was insufficiently sampled and not preset
            elif InAncCat not in InPresetCatList:
                UnsampledList.append(InAncCat)
        # End loop for this ancillary class...
    # End divide over selection method...
    AncFeatureClass, SelectMethod, Percent, SampleMin, PopLayer, InAncCatList = None, None, None, None, None, None
    AncLayer, OutAncCatList, PopWorkingTableName, UninhabList, POPAREAWhereClause = None, None, None, None, None
    
    # For each ancillary class (listed in the REPCAT field) calculate sum of population and area and statistics
    # - (count, mean, min, max, stddev) of densities further analysis
    print "Calculating statistics for selected classes..."
    gp.AddMessage("Calculating statistics for selected classes...")
    # Make sure there are representative classes.
    WhereClause = Wrp[oFC][0] + "REPCAT" + Wrp[oFC][1] + " IS NOT NULL"
    if AncCatFType == "String":
        WhereClause = WhereClause + " AND " + Wrp[oFC][0] + "REPCAT" + Wrp[oFC][1] + " <> ''"
    aPopWTableView = gp.MakeTableView_management(PopWorkingTable, "aPopWTableView", WhereClause)
    AncDensTableName = NameCheck("SamplingSummaryTable")
    if gp.GetCount_management(aPopWTableView):
        # Statistics_analysis (in_table, out_table, statistics_fields, case_field)
        AncDensTable = gp.Statistics_analysis(aPopWTableView, AncDensTableName + TableSuffix, PopFieldName[:fL] + " SUM; " + AreaFieldName[:fL] + " SUM; DENSITY MEAN; DENSITY MIN; DENSITY MAX; DENSITY STD; POPAREA SUM; PDENSITY MEAN; PDENSITY MIN; PDENSITY MAX; PDENSITY STD;" , "REPCAT")
        gp.AddField_management(AncDensTable, "CLASSDENS", "DOUBLE")
        # Calculate an initial population estimate for each polygon in this class by multiplying the representative class densities by the polygon areas
        gp.CalculateField_management(AncDensTable, "CLASSDENS", "[" + ("SUM_" + PopFieldName)[:fL] + "] / [" + "SUM_POPAREA"[:fL] + "]")
        # Add a field that designates these classes as "Sampled"
        gp.AddField_management(AncDensTable, "METHOD", "TEXT", "", "", "7")
        gp.CalculateField_management(AncDensTable, "METHOD", '"Sampled"')
        gp.AddField_management(AncDensTable, "PreDens", "DOUBLE")
        # For all sampled classes that are not preset, calculate a population estimate for every intersected polygon by joining the AncDensTable and multiplying the class density by the polygon area.
        print "Calculating first population estimate for sampled and preset classes..."
        gp.AddMessage("Calculating first population estimate for sampled and preset classes...")
        gp.AddJoin_management(OutLayer, AncCatName[:fL], AncDensTable, "REPCAT", "KEEP_COMMON")
        if PresetTable and PresetTable <> "#":
            WhereClause = Wrp[oFC][0] + OutLayerName + "." + AncCatName[:fL] + Wrp[oFC][1] + " NOT IN (" + ", ".join(InPresetCatList) + ")"
            OutLayer = gp.SelectLayerByAttribute_management(OutLayer, "NEW_SELECTION", WhereClause)
        else:
            OutLayer = gp.SelectLayerByAttribute_management(OutLayer, "CLEAR_SELECTION")
        gp.CalculateField_management(OutLayer, OutLayerName + ".POPEST", "[" + OutLayerName + "." + NewArea + "] * [" + AncDensTableName + ".CLASSDENS]")
        gp.RemoveJoin_management(OutLayer, AncDensTableName)
    else:
        # CreateTable_management (out_path, out_name, template, config_keyword)
        AncDensTable = CreateTable_management(OutWorkspace, AncDensTableName + TableSuffix)
        gp.AddField_management(AncDensTable, "REPCAT", FieldProps[0], "", "", FieldProps[1])
        gp.AddField_management(AncDensTable, "CLASSDENS", "DOUBLE")
        gp.AddField_management(AncDensTable, "METHOD", "TEXT", "", "", "7")
        gp.AddField_management(AncDensTable, "PreDens", "DOUBLE")

    if PresetTable and PresetTable <> "#":
        print "Adding preset values to the summary table..."
        gp.AddMessage("Adding preset values to the summary table...")
        # Now, for the preset classes, calculate a population estimate for every intersected polygon by joining the Preset Table and multiplying the preset density by the polygon area
        OutLayer = gp.SelectLayerByAttribute_management(OutLayer, "CLEAR_SELECTION")
        gp.AddJoin_management(OutLayer, AncCatName[:fL], PresetTable, AncCatName, "KEEP_COMMON")
        PresetTableName = GetLayerName(PresetTable)
        gp.CalculateField_management(OutLayer, "POPEST", "[" + OutLayerName + "." + NewArea + "] * [" + PresetTableName + "." + PresetField + "]")
        gp.RemoveJoin_management(OutLayer,PresetTableName)
        # Add these preset values to the AncDensTable for comparison purposes.
        i = 0
        for InPresetCat in InPresetCatList:
            AncDensTableView = gp.MakeTableView_management(AncDensTable, "AncDensTableView", Wrp[oFC][0] + "REPCAT" + Wrp[oFC][1] + " = " + InPresetCat)
            Count = gp.GetCount_management(AncDensTableView)
            if Count > 0:
                gp.CalculateField_management(AncDensTableView, "PreDens", PresetValList[i])
                gp.CalculateField_management(AncDensTableView, "METHOD", "'Preset'")
            else:
                rows = gp.InsertCursor(AncDensTable)
                row = rows.NewRow()
                row.PreDens = PresetValList[i]
                row.METHOD = "Preset"
                row.REPCAT = OutPresetCatList[i]
                rows.InsertRow(row)
                row, rows = None, None
            i = i + 1
        OutPresetCatList, InPresetCatList, i = None, None, None
    RemoveNulls(OutLayer,"POPEST",Wrp,OutFeatureClass)
    GetLayerName, AreaFieldName, PresetTable, PresetField = None, None, None, None
    OutWorkspace, FieldProps, PresetValList, AncDensTableName = None, None, None, None
    
    # Intelligent areal weighting for unsampled classes
    # - for every source population unit sum the initial population estimates and compare
    # - the result to the actual count for the unit. Distribute any residual population 
    # - among the remaining unsampled inhabited dasymetric polygons
    print "Performing intelligent areal weighting for unsampled classes..."
    gp.AddMessage("Performing intelligent areal weighting for unsampled classes...")
    if UnsampledList:
        OutLayer = gp.SelectLayerByAttribute_management(OutLayer, "NEW_SELECTION", Wrp[oFC][0] + AncCatName[:fL] + Wrp[oFC][1] + " IN (" + ", ".join(UnsampledList) + ")")
        gp.CalculateField_management(OutLayer, "REMAREA", "[" + NewArea + "]")
        RemoveNulls(OutLayer,"REMAREA",Wrp,OutFeatureClass)
        RemainderTableName = NameCheck("RemainderTable")
        OutLayer = gp.SelectLayerByAttribute_management(OutLayer, "CLEAR_SELECTION") # To clear previous selection set
        RemainderTable = gp.Frequency_analysis(OutLayer, RemainderTableName + TableSuffix, PopKeyName[:fL] + ";" + PopFieldName[:fL], "POPEST;REMAREA")
        gp.AddField_management(RemainderTable, "POPDIFF", "DOUBLE")
        gp.CalculateField_management(RemainderTable, "POPDIFF", "[" + PopFieldName[:fL] + "] - [POPEST]")
        gp.AddJoin_management(OutLayer, PopKeyName[:fL], RemainderTable, PopKeyName[:fL], "KEEP_COMMON")
        WhereClause = Wrp[oFC][0] + RemainderTableName + ".POPDIFF" + Wrp[oFC][1] + " > 0"
        WhereClause = WhereClause + " AND " + Wrp[oFC][0] + RemainderTableName + ".REMAREA" + Wrp[oFC][1] + " <> 0"
        gp.SelectLayerByAttribute_management(OutLayer, "NEW_SELECTION", WhereClause)
        gp.CalculateField_management(OutLayer, OutLayerName + ".POPEST", "[" + RemainderTableName + ".POPDIFF] * [" + OutLayerName + ".REMAREA] / [" + RemainderTableName + ".REMAREA]")
        gp.RemoveJoin_management(OutLayer, RemainderTableName)
        gp.Delete_management(RemainderTable)

        # Calculate population density values for these unsampled classes
        # - for every unsampled ancillary class, sum total area and total population estimated using intelligent areal weighting.  Calculate class representative density.
        WhereClause = Wrp[oFC][0] + AncCatName[:fL] + Wrp[oFC][1] + " IN (" + ", ".join(UnsampledList) + ")"
        WhereClause = WhereClause + " AND " + Wrp[oFC][0] + PopFieldName[:fL] + Wrp[oFC][1] + " <> 0"
        OutLayer = gp.SelectLayerByAttribute_management(OutLayer, "NEW_SELECTION", WhereClause)
        AncDensTable2Name = NameCheck("AncDensTable2")
        AncDensTable2 = gp.Frequency_analysis(OutLayer, AncDensTable2Name + TableSuffix, AncCatName[:fL], "POPEST;POPAREA")
        gp.AddField_management(AncDensTable2, "CLASSDENS", "DOUBLE")
        WhereClause = Wrp[oFC][0] + "POPAREA" + Wrp[oFC][1] + " > 0"
        AncDensTable2View = gp.MakeTableView_management(AncDensTable2, "AncDensTable2View", WhereClause)
        gp.CalculateField_management(AncDensTable2View, "CLASSDENS", "[POPEST] / [POPAREA]")
        gp.AddJoin_management(OutLayer, AncCatName[:fL], AncDensTable2, AncCatName[:fL], "KEEP_COMMON")
        # - Again recalculate population estimate field (POPEST) using new representative density for final stats analysis
        gp.CalculateField_management(OutLayer, OutLayerName + ".POPEST", "[" + OutLayerName + "." + NewArea + "] * [" + AncDensTable2Name + ".CLASSDENS]")
        gp.RemoveJoin_management(OutLayer, AncDensTable2Name)
        def GetNumVals(Table, Field):
            InList = []
            if Table:
                Rows = gp.SearchCursor(Table)
                Row = Rows.Next()
                while Row:
                    InList.append(Row.GetValue(Field))
                    Row = Rows.Next()
                Rows, Row = None, None
            return InList
        # - Lastly, add these IAW values to the AncDensTable
        IAWValList = GetNumVals(AncDensTable2, "CLASSDENS")
        UnsampledList = GetNumVals(AncDensTable2, AncCatName)
        for i in range(0, len(IAWValList)):
            rows = gp.InsertCursor(AncDensTable)
            row = rows.NewRow()
            row.CLASSDENS = IAWValList[i]
            row.METHOD = "IAW"
            row.REPCAT = UnsampledList[i]
            rows.InsertRow(row)
            del row, rows
        gp.Delete_management(AncDensTable2)
        IAWValList, UnsampledList = None, None

    # Perform final calculations to ensure pycnophylactic integrity
    print "Performing final calculations to ensure pycnophylactic integrity..."
    gp.AddMessage("Performing final calculations to ensure pycnophylactic integrity...")
    # For each population source unit, sum the population estimates,
    # - which do not necessarily sum to the actual population of the source,
    # - and use the ratio of the estimates to the estimated total to distribute the actual total.
    PopEstSumTableName = NameCheck("PopEstSumTable")
    OutLayer = gp.SelectLayerByAttribute_management(OutLayer, "CLEAR_SELECTION") # To clear previous selection set
    PopEstSumTable = gp.Frequency_analysis(OutLayer, PopEstSumTableName + TableSuffix, PopKeyName[:fL], "POPEST")
    gp.AddJoin_management(OutLayer, PopKeyName[:fL], PopEstSumTable, PopKeyName[:fL], "KEEP_COMMON")
    # The ratio of each dasymetric unit's population estimate to this sum is called the total fraction
    WhereClause = Wrp[oFC][0] + PopEstSumTableName + ".POPEST" + Wrp[oFC][1] + " <> 0"
    OutLayer = gp.SelectLayerByAttribute_management(OutLayer, "NEW_SELECTION", WhereClause)
    # [TOTALFRACT] = [POPEST] / PopEstSum
    gp.CalculateField_management(OutLayer, "TOTALFRACT", "[" + OutLayerName + ".POPEST] / [" + PopEstSumTableName + ".POPEST]")
    gp.RemoveJoin_management(OutLayer, PopEstSumTableName)
    gp.Delete_management(PopEstSumTable)
    # The total fraction times the actual population is the true dasymetric estimate
    # [NEWPOP] = [TOTALFRACT] * [source unit pop] = [source unit pop] * [POPEST] / PopEstSum 
    gp.CalculateField_management(OutLayer, "NEWPOP", "[" + PopFieldName[:fL] + "] * [TOTALFRACT]")
    # Calculate a final density value for statistical purposes
    WhereClause = Wrp[oFC][0] + NewArea + Wrp[oFC][1] + " <> 0"
    OutLayer = gp.SelectLayerByAttribute_management(OutLayer, "NEW_SELECTION", WhereClause)
    # [NEWDENSITY] = [NEWPOP] / NewArea
    gp.CalculateField_management(OutLayer, "NEWDENSITY", "[NEWPOP] / [" + NewArea + "]")
    
    # Lastly create an official output statistics table
    print "Creating a final summary table"
    gp.AddMessage("Creating a final summary table")
    FinalSummaryTableName = NameCheck("FinalSummaryTable")
    gp.Statistics_analysis(OutFeatureClass, FinalSummaryTableName + TableSuffix, "NEWPOP SUM; NEWDENSITY MEAN; NEWDENSITY MIN; NEWDENSITY MAX; NEWDENSITY STD", AncCatName[:fL])
    FinalSummaryTableName, NewArea, OutFeatureClass, WhereClause, PopFieldName = None, None, None, None, None
    PopEstSumTableName, RemoveNulls, PresetTableName, OutLayerName, PopKeyName = None, None, None, None, None
    RemainderTableName, AncCatName, PopFieldName, NameCheck, GetValues = None, None, None, None, None
    Wrp, fL, pFC, aFC, oFC = None, None, None, None, None 
    del gp
except:  # If an error occurred print the message to the screen
    print gp.GetMessages(2)
    gp.AddMessage(gp.GetMessages(2))
    del gp

Hi @danabauer might you be willing to chat with me about this next week? I'm at the Census Bureau and we're starting a new OSS project wherein something to this end would be really handy...