VERSION 1.0 CLASS
BEGIN
  MultiUse = -1  'True
END
Attribute VB_Name = "ThisDocument"
Attribute VB_GlobalNameSpace = False
Attribute VB_Creatable = False
Attribute VB_PredeclaredId = True
Attribute VB_Exposed = True
'Raster-Based Dasymetric Mapping Routine
'Developed at the University of Colorado in Boulder
'by Torrin Hultgren and Jeremy Mennis

'Apologies for the lack of built-in user interaction
'I don't have the time to code it right now, and since the future lies in python scripts in ModelBuilder,
'if I did have the time, I'd spend it on porting this script to python rather than coding something clumsy

Public Sub Dasymetric_Raster()
  On Error GoTo ERH
  Dim pMXDoc As esriArcMapUI.IMxDocument
  Set pMXDoc = ThisDocument
'These are the user defined variables, please change as necessary:
  
  Dim Percent As Variant
'This is how you set the selection criteria for representative source units
'It is the percent of the source unit that must be covered by an ancillary class (from 0 to 1)
  Percent = 0.9
  Dim popFieldName As Variant
'This is the name of the field in your population dataset that contains population count data
  popFieldName = "POP2000"
 
  Dim PresetClasses As New Collection
  Dim PresetDensities As New Collection
'This is how you preset certain classes - for example, a water class usually has 0 population
'Put your ancillary class code/value in the quotes - in the example, the water class has a raster value of 5
  PresetClasses.Add "5" 'Water
'Put your preset density in the first quotes, and the above class in the second quotes
  PresetDensities.Add "0", "5"
  PresetClasses.Add "3" 'Non-residential Developed
  PresetDensities.Add "0", "3"
'You may cut and paste the above lines as many times as necessary to include all classes
  
  Dim PathName As String
'Please type the pathname where the input data are stored
  PathName = "C:\users\mhessenflow\documents\rachels_dasymetric\output2"
  
  Dim pPopLayer As esriCarto.ILayer
'This is how you designate which layers in the table of contents are your input data
'Arc numbers the layers from top to bottom starting with zero.
'You can either put the vector population layer first and the raster ancillary layer second
'Or just put the appropriate numbers in the parentheses below.
  Set pPopLayer = pMXDoc.FocusMap.Layer(0)
  Dim pAncLayer As esriCarto.ILayer
  Set pAncLayer = pMXDoc.FocusMap.Layer(1)

'End of user defined variables
  
  'Time the script for reference purposes
  Dim StartTime As Variant
  Dim EndTime As Variant
  Dim Elapsed As Variant
  StartTime = Now
  
  Dim pErrorTrack As String
  pErrorTrack = "Initial Raster Procedures"
  Dim pConvOp As esriGeoAnalyst.IConversionOp
'  Set pConvOp = New esriGeoAnalyst.RasterConversionOp
        ' !!! MLH We set the pixel extent for the raster environment below,
        '         but we have already created this operator. Does that cause
        '         the problem with undefined pixel extents on the first run?
  Dim pPopFeatureLayer As esriCarto.IFeatureLayer
  Set pPopFeatureLayer = pPopLayer
  Dim pPopFeatureClass As esriGeoDatabase.IFeatureClass
  Set pPopFeatureClass = pPopFeatureLayer.FeatureClass
  Dim pPopGeoDataset As esriGeoDatabase.IGeoDataset
  Set pPopGeoDataset = pPopFeatureClass

  

  Dim pRasterLayer As esriCarto.IRasterLayer
  Set pRasterLayer = pAncLayer
  Dim pAncRaster As esriGeoDatabase.IRaster
  Set pAncRaster = pRasterLayer.Raster
  ' Get cell from pInputRaster and extent from pop dataset
  Dim pRasterProp As esriDataSourcesRaster.IRasterProps
  Dim pExtent As esriGeometry.IEnvelope
  Dim vCell As Double
  vCell = 5 'default value to avoid problems...
  Set pRasterProp = pAncRaster
  Set pExtent = pPopGeoDataset.Extent
  vCell = (pRasterProp.MeanCellSize.X + pRasterProp.MeanCellSize.Y) / 2
  
  ' Create a Rasterworkspace
  Dim pRasterWS As esriGeoDatabase.IWorkspace
  Set pRasterWS = SetRasterWorkspace(PathName)

  ' Create RasterAnalysis environment
  Dim pEnv As esriGeoAnalyst.IRasterAnalysisEnvironment
  Set pEnv = New esriGeoAnalyst.RasterAnalysis
  pEnv.SetCellSize esriGeoAnalyst.esriRasterEnvValue, vCell
  pEnv.SetExtent esriGeoAnalyst.esriRasterEnvValue, pExtent
  Set pEnv.OutWorkspace = pRasterWS
  ' Set to default so that it works for all the Ops in this session
  pEnv.SetAsNewDefaultEnvironment

  Set pConvOp = New esriGeoAnalyst.RasterConversionOp ' !!! MLH Maybe this line should go here.
  
'Convert the population enumeration districts to raster
  Dim pPopRasterDataset As esriGeoDatabase.IRasterDataset
  Set pPopRasterDataset = pConvOp.ToRasterDataset(pPopGeoDataset, "GRID", pRasterWS, "PopRaster")
  
'This section adds the rasterized population dataset to the map.
'It's not really a necessary step and could be omitted.
  Dim pPopRasterLayer As esriCarto.IRasterLayer
  Set pPopRasterLayer = New esriCarto.RasterLayer
  pPopRasterLayer.CreateFromDataset pPopRasterDataset
  pMXDoc.FocusMap.AddLayer pPopRasterLayer
'End of optional section __________________________________________________
    
  Dim ConversionTime As Variant
  ConversionTime = Now
    
'Combine the rasterized population districts with the landcover raster
  Dim pPopRBCollection As esriDataSourcesRaster.IRasterBandCollection
  Set pPopRBCollection = pPopRasterLayer.Raster
  Dim pAncRBCollection As esriDataSourcesRaster.IRasterBandCollection
  Set pAncRBCollection = pRasterLayer.Raster
  Dim pPopRasterBand As esriDataSourcesRaster.IRasterBand
  Set pPopRasterBand = pPopRBCollection.Item(0)
  Dim pAncRasterBand As esriDataSourcesRaster.IRasterBand
  Set pAncRasterBand = pAncRBCollection.Item(0)

' Create a raster bandcollection
  Dim pRBCollection As esriDataSourcesRaster.IRasterBandCollection
  Set pRBCollection = New esriDataSourcesRaster.Raster
  pRBCollection.AppendBand pPopRasterBand
  pRBCollection.AppendBand pAncRasterBand
   
  Dim pInputDataset As esriGeoDatabase.IGeoDataset
  Set pInputDataset = pRBCollection
' Declare the output raster object
  Dim pCombDataset As esriGeoDatabase.IGeoDataset
' Create the RasterLocalOp object
  Dim pLocalOp As esriSpatialAnalyst.ILocalOp
  Set pLocalOp = New esriSpatialAnalyst.RasterLocalOp
' Calls the method
  Set pCombDataset = pLocalOp.Combine(pInputDataset)

  Dim pCombRastLyr As esriCarto.IRasterLayer
  Set pCombRastLyr = New esriCarto.RasterLayer
  With pCombRastLyr
    .CreateFromRaster pCombDataset
    .Name = "Dasymetric_result"
  End With
  pMXDoc.FocusMap.AddLayer pCombRastLyr
  Dim pCombDispTable As esriCarto.IDisplayTable
  Set pCombDispTable = pCombRastLyr
  Dim pCombTable As esriGeoDatabase.ITable
  Set pCombTable = pCombDispTable.DisplayTable

  'Save the combined dataset
  Dim pDS As esriGeoDatabase.IDataset
  Dim pCombRast As esriGeoDatabase.IRaster
  Set pCombRast = pCombDataset
  Dim pCombRBColl As esriDataSourcesRaster.IRasterBandCollection
  Set pCombRBColl = pCombRast
  Set pDS = pCombRBColl.SaveAs("Dasy_Rast", pRasterWS, "GRID")

'Join the population table to the new raster output table
' ++ Create the MemoryRelationshipClass that defines what is to be joined
  Dim pMemRelClassFact As esriGeoDatabase.IMemoryRelationshipClassFactory
  Set pMemRelClassFact = New esriGeoDatabase.MemoryRelationshipClassFactory
  Dim pRelClass As esriGeoDatabase.IRelationshipClass
  Set pRelClass = pMemRelClassFact.Open("PopToRaster", pPopFeatureClass, "FID", pCombTable, "Popraster", _
                                        "forward", "backward", esriGeoDatabase.esriRelCardinalityOneToMany)
 ' ++ Perform the join
  Dim pRelQueryTableFact As esriGeoDatabase.IRelQueryTableFactory
  Dim pRelQueryTab As esriGeoDatabase.ITable
  Set pRelQueryTableFact = New esriGeoDatabase.RelQueryTableFactory
  Set pCombTable = pRelQueryTableFact.Open(pRelClass, False, Nothing, Nothing, "", True, False)

' Add it to the relates for the feature layer in Map
  Dim pRelClassCollEdit As esriCarto.IRelationshipClassCollectionEdit
  Set pRelClassCollEdit = pCombRastLyr
  pRelClassCollEdit.AddRelationshipClass pRelClass
  
  'Display the join
  Dim pDispRC As esriCarto.IDisplayRelationshipClass
  Set pDispRC = pCombRastLyr
  pDispRC.DisplayRelationshipClass pRelClass, esriGeoDatabase.esriLeftInnerJoin
 
'Export the resulting table to a dbf file
  Set pDS = Nothing
  Set pDS = pCombTable
  Dim pDSName As esriGeoDatabase.IDatasetName
  Set pDSName = pDS.FullName
  Dim pOutDSName As esriGeoDatabase.IDatasetName
  Set pOutDSName = New esriGeoDatabase.TableName
  pOutDSName.Name = "DasyTable"
  Dim pWkSpFactory As esriGeoDatabase.IWorkspaceFactory
  Set pWkSpFactory = New esriDataSourcesFile.ShapefileWorkspaceFactory
  Dim pWkSp As esriGeoDatabase.IWorkspace
  Set pWkSp = pWkSpFactory.OpenFromFile(PathName & "\", 0)
  Dim pWkSpDS As esriGeoDatabase.IDataset
  Set pWkSpDS = pWkSp
  Dim pWkSpName As esriGeoDatabase.IWorkspaceName
  Set pWkSpName = pWkSpDS.FullName
  Set pOutDSName.WorkspaceName = pWkSpName

  ' Export (Selection is ignored)
  Dim pExpOp As esriGeoDatabaseUI.IExportOperation
  Set pExpOp = New esriGeoDatabaseUI.ExportOperation
  pExpOp.ExportTable pDSName, Nothing, Nothing, pOutDSName, Application.hWnd
  
  ' add the table to map
  Dim pName As esriSystem.IName
  Dim pDasyTable As esriGeoDatabase.ITable
  Dim pStTab As esriCarto.IStandaloneTable
  Dim pStTabColl As esriCarto.IStandaloneTableCollection
  Set pName = pOutDSName
  Set pDasyTable = pName.Open
  Set pStTab = New esriCarto.StandaloneTable
  Set pStTab.Table = pDasyTable
  Set pStTabColl = pMXDoc.FocusMap
  pStTabColl.AddStandaloneTable pStTab
  pMXDoc.UpdateContents

  Dim IntersectTime As Variant
  IntersectTime = Now

'Create new fields
  CreateFields pDasyTable
  pErrorTrack = "Pre-Calculations"
'Collect unique tract ids
  Dim pPopIDs As New Collection
  Dim pQueryfilter As esriGeoDatabase.IQueryFilter
  Set pQueryfilter = Nothing
  Dim pPopIDField As Variant
  pPopIDField = "PopRaster"
  Set pPopIDs = GetUnique(pQueryfilter, pDasyTable, pPopIDField)

'Collect unique ancillary values
  Dim pAncIDs As New Collection
  Set pQueryfilter = Nothing
  Dim pAncIDField As Variant
  pAncIDField = VBA.Left(pAncLayer.Name, 10)
  Set pAncIDs = GetUnique(pQueryfilter, pDasyTable, pAncIDField)
  
    'Need to start an edit session to update the table
  Dim pEditor As esriEditor.IEditor
  Dim pID As New esriSystem.UID
  pID = "esriCore.Editor"
  Set pEditor = Application.FindExtensionByCLSID(pID)
  pEditor.StartEditing pWkSp
  
'Perform pre-selection calculations
  pEditor.StartOperation
  PreCalc pPopIDs, PresetClasses, PresetDensities, pPopIDField, pAncIDField, pDasyTable, popFieldName
  ' !!! MLH: Added popFieldName to argument list
  pEditor.StopOperation "Pre-Calculations"
  
  Dim PreCalcTime As Variant
  PreCalcTime = Now
  pErrorTrack = "Selection"
'Select representative tracts, probably by percent method comparing the count column with the populated area column
'Don't need to perform any special selection, since it can all be done with a simple queryfilter.
'For each ancillary class, find rows where Count > Percent * INHAB_CNT
'Select the rows and calculate density statistics on them.
'Calculate representative density by dividing total representative population by total representative area.  Save stats for end.
'Calculate densities for representative tracts, i.e. total pop/poparea for statsfile purposes?
  Dim pStats As esriGeoDatabase.IDataStatistics
  Dim pResults As esriSystem.IStatisticsResults
  Dim pCursor As esriGeoDatabase.ICursor
  Dim pCalc As esriGeoDatabaseUI.ICalculator
  
  Dim pPDArray As New Collection
  Dim pBegCount As New Collection
  Dim pBegMean As New Collection
  Dim pBegMin As New Collection
  Dim pBegMax As New Collection
  Dim pBegStdDev As New Collection
  Dim pTotAreaArray As New Collection
  Dim pTotPopArray As New Collection
  Dim pFlagClass As New Collection
  Dim pPopArray As New Collection
  Dim pCalcMethod As New Collection
  Dim pUnsampled As New Collection
  
  Dim pSampledArea As Double
  Dim pTotArea As Double
  Dim pTotPop As Double
  Dim pClass As Variant
  For Each pClass In pAncIDs
    Set pQueryfilter = New esriGeoDatabase.QueryFilter
    pQueryfilter.WhereClause = pAncIDField & " = " & pClass & " AND Count >= " & Percent & " * INHAB_CNT"
    'Now calculate the statistics for the selected tracts
    Set pCursor = pDasyTable.Search(pQueryfilter, False)
    Set pStats = New esriGeoDatabase.DataStatistics
    Set pStats.Cursor = pCursor
    pStats.Field = "INHAB_CNT" 'Want to be sure to only include populated area
    Set pResults = pStats.Statistics
    pTotArea = pResults.Sum
    pTotAreaArray.Add pTotArea, pClass
    pSampledArea = pSampledArea + pTotArea
    Set pCursor = Nothing
    Set pStats = Nothing
    Set pResults = Nothing
    
    Set pCursor = pDasyTable.Search(pQueryfilter, False)
    Set pStats = New esriGeoDatabase.DataStatistics
    Set pStats.Cursor = pCursor
    pStats.Field = popFieldName
    Set pResults = pStats.Statistics
    pTotPop = pResults.Sum
    pTotPopArray.Add pTotPop, pClass

    'Now collect the density statistics for analysis
    Set pCursor = pDasyTable.Search(pQueryfilter, False)
    Set pStats = New esriGeoDatabase.DataStatistics
    Set pStats.Cursor = pCursor
    pStats.Field = "INHAB_DENS"
    Set pResults = pStats.Statistics
    Dim BegCount As Double
    BegCount = pResults.Count
    pBegCount.Add BegCount, pClass
    Dim BegMean As Double
    BegMean = pResults.Mean
    pBegMean.Add BegMean, pClass
    Dim BegMin As Double
    BegMin = pResults.Minimum
    If BegCount = 0 Then BegMin = 0
    pBegMin.Add BegMin, pClass
    Dim BegMax As Double
    BegMax = pResults.Maximum
    If BegCount = 0 Then BegMax = 0
    pBegMax.Add BegMax, pClass
    Dim BegStdDev As Double
    BegStdDev = pResults.StandardDeviation
    pBegStdDev.Add BegStdDev, pClass

    Dim pNewDensity As Double
    If pDasyTable.RowCount(pQueryfilter) > 0 And pTotArea > 0 Then 'To avoid errors if there are no tract polygons selected
      pNewDensity = pTotPop / pTotArea
    Else: pNewDensity = 0
    End If
  
    'If this is one of our preset classes, flag it through the unsampled classes selection
    Dim pPresetC As Variant
    For Each pPresetC In PresetClasses
      If pClass = pPresetC Then
        pNewDensity = -1 'To exclude from our insufficient sample code
        pCalcMethod.Add "Preset", pClass 'Mark class as preset
        Exit For
      End If
    Next

    'Flag classes and regions that were insufficiently sampled.
    If BegCount < 2 And pNewDensity <> -1 Then
      pNewDensity = 0 'Set to zero so it doesn't contribute to the total
      pUnsampled.Add pClass 'Place in array for unsampled processing.
      pCalcMethod.Add "Areal Weight", pClass 'Mark class as having used the areal weight method
    ElseIf BegCount > 1 And pNewDensity <> -1 Then
      pCalcMethod.Add "Sampled", pClass
    End If
    If pNewDensity = -1 Then pNewDensity = PresetDensities(pClass) 'Remove flag, set to preset value.

    Dim pPDTotal As Double
    'put that landcover population density value into an array pPDArray,
    pPDArray.Add pNewDensity, pClass
    'sum all population density values into the total sum density pPDTotal
    pPDTotal = pPDTotal + pNewDensity
  
    'Calculate a population estimate for all dasymetric polys with this class by multiplying the class density by the poly area
    If pNewDensity <> 0 Then
      Set pCalc = New esriGeoDatabaseUI.Calculator
      Set pQueryfilter = New esriGeoDatabase.QueryFilter
      pQueryfilter.WhereClause = pAncIDField & " = " & pClass
      Set pCursor = pDasyTable.Update(pQueryfilter, False)
      With pCalc
        Set .Cursor = pCursor
            .Expression = "[Count] * " & pNewDensity
            .Field = "POPEST"
      End With
      pCalc.Calculate
      Set pCalc = Nothing
    End If
  
  Next 'next ancillary class

  Dim SelectionTime As Variant
  SelectionTime = Now
'Handle unsampled classes
  
  pErrorTrack = "initial population estimation"
  'Method for estimating unsampled ancillary class densities
  Dim pEstPop, pExtraPop As Double
  pEditor.StartOperation
  PopulationEst pPopIDs, pDasyTable, popFieldName, "Popraster", pErrorTrack
  pEditor.StopOperation "Population Estimation"
  
'Now calculate the population densities of the unsampled classes based on these estimates
  Dim pEstArea As Double
  For Each pClass In pUnsampled
    Set pCursor = Nothing
    Set pQueryfilter = New esriGeoDatabase.QueryFilter
    pQueryfilter.WhereClause = pAncIDField & " = " & pClass
    Set pCursor = pDasyTable.Search(pQueryfilter, False)
    Set pStats = New esriGeoDatabase.DataStatistics
    Set pStats.Cursor = pCursor
    pStats.Field = "POPEST"
    Set pResults = pStats.Statistics
    pEstPop = pResults.Sum
    Set pCursor = pDasyTable.Search(pQueryfilter, False)
    Set pStats = New esriGeoDatabase.DataStatistics
    Set pStats.Cursor = pCursor
    pStats.Field = "Count"
    Set pResults = pStats.Statistics
    pEstArea = pResults.Sum
    If pEstArea > 0 Then 'Because dividing by zero is bad, although this should never occur.
      pNewDensity = pEstPop / pEstArea
      'put the new landcover population density value into pPDArray,
      pPDArray.Remove pClass
      pPDArray.Add pNewDensity, pClass
      'add this new density into the total sum density pPDTotal
      pPDTotal = pPDTotal + pNewDensity
    End If
    'Put class density into a column
    pEditor.StartOperation
    Set pCalc = New esriGeoDatabaseUI.Calculator
    pNewDensity = pPDArray(pClass)
    'Add this value to the "CLSDENSITY" field for every row with that landcover,
    Set pQueryfilter = New esriGeoDatabase.QueryFilter
    pQueryfilter.WhereClause = pAncIDField & " = " & pClass
    Set pCursor = pDasyTable.Update(pQueryfilter, False)
    With pCalc
      Set .Cursor = pCursor
      .Expression = "[Count] * " & pNewDensity
      .Field = "POPEST"
    End With
    pCalc.Calculate
    Set pCalc = Nothing
    pEditor.StopOperation "Re-estimate Population"
  Next 'pClass in pUnsampled
  
  Dim SAMTime As Variant
  SAMTime = Now
  
  'Final Calculations are standardized, although there may be two options.
  pErrorTrack = "final calculations"
  Dim pRegQueryfilter As New esriGeoDatabase.QueryFilter
  pEditor.StartOperation
  FinalCalc pDasyTable, pRegQueryfilter, pPopIDs, "Popraster", popFieldName
  pEditor.StopOperation "Final Calculations"

  'All done with editing tables.
  pEditor.StopEditing (True)
  
  Dim FinalCalcTime As Variant
  FinalCalcTime = Now
  
  'Now generate a csv file with summary statistics for each ancillary class
  pErrorTrack = "stat file creation"
  Dim pStatCount, pMean, pMin, pMax, pStdDev As Variant
  Dim BegTotPop, BegTotArea, BegPopPerc, BegAreaPerc, PercSamArea As Variant
  Dim pCatArea, pCatPercArea, pCatPop, pCatPercPop As Variant
  Dim FinTotPop, FinTotArea As Variant
  Dim pSummary As String
  'Calculate the Original and final total population for comparison purposes.
  Set pCursor = pPopFeatureLayer.Search(Nothing, False)
  Set pStats = New esriGeoDatabase.DataStatistics
  Set pStats.Cursor = pCursor
  pStats.Field = popFieldName
  Set pResults = pStats.Statistics
  BegTotPop = pResults.Sum
  Set pCursor = pPopRasterBand.AttributeTable.Search(Nothing, False)
  Set pStats = New esriGeoDatabase.DataStatistics
  Set pStats.Cursor = pCursor
  pStats.Field = "Count"
  Set pResults = pStats.Statistics
  BegTotArea = pResults.Sum
  Set pCursor = pDasyTable.Search(Nothing, False)
  Set pStats = New esriGeoDatabase.DataStatistics
  Set pStats.Cursor = pCursor
  pStats.Field = "NEWPOP"
  Set pResults = pStats.Statistics
  FinTotPop = pResults.Sum
  Set pCursor = pDasyTable.Search(Nothing, False)
  Set pStats = New esriGeoDatabase.DataStatistics
  Set pStats.Cursor = pCursor
  pStats.Field = "TRCTCOUNT"
  Set pResults = pStats.Statistics
  FinTotArea = pResults.Sum
  
  Dim output As Object
  Dim pSuffix As Variant
  Dim outtext As Object
  Dim pAncArray As New Collection
  Set output = CreateObject("Scripting.FileSystemObject")
  Set outtext = output.CreateTextFile(PathName & "\Dasymetric_Stats" & pSuffix & ".csv", True)
  pSummary = "Ancillary Class,Rep Sample Count,Rep % Area,% Sampled Area,Rep % Pop,Rep Mean,Rep Min,Rep Max,Rep StdDev,Est Density," & _
             "Sample Count,% Total Area,% Total Pop,Mean Density,Min Density,Max Density,StdDev,Calc Method"
  outtext.writeline (pSummary)
  Set pQueryfilter = New esriGeoDatabase.QueryFilter
  Set pAncArray = GetUnique(pQueryfilter, pDasyTable, pAncIDField)
  For Each pClass In pAncArray
    Set pQueryfilter = New esriGeoDatabase.QueryFilter
    pQueryfilter.WhereClause = pAncIDField & "= " & pClass
    Set pCursor = pDasyTable.Search(pQueryfilter, False)
    Set pStats = New esriGeoDatabase.DataStatistics
    Set pStats.Cursor = pCursor
    pStats.Field = "NEWDENSITY"
    Set pResults = pStats.Statistics
    pStatCount = pResults.Count
    pMean = pResults.Mean
    pMin = pResults.Minimum
    pMax = pResults.Maximum
    pStdDev = pResults.StandardDeviation
    Set pCursor = pDasyTable.Search(pQueryfilter, False)
    Set pStats = New esriGeoDatabase.DataStatistics
    Set pStats.Cursor = pCursor
    pStats.Field = "TRCTCOUNT"
    Set pResults = pStats.Statistics
    pCatArea = pResults.Sum
    pCatPercArea = pCatArea / FinTotArea
    Set pCursor = pDasyTable.Search(pQueryfilter, False)
    Set pStats = New esriGeoDatabase.DataStatistics
    Set pStats.Cursor = pCursor
    pStats.Field = "NEWPOP"
    Set pResults = pStats.Statistics
    pCatPop = pResults.Sum
    pCatPercPop = pCatPop / FinTotPop
    BegPopPerc = pTotPopArray(pClass) / BegTotPop
    BegAreaPerc = pTotAreaArray(pClass) / BegTotArea
    If pSampledArea <> 0 Then
      PercSamArea = pTotAreaArray(pClass) / pSampledArea
    Else
      PercSamArea = 0
    End If
    pSummary = pClass & ", " & pBegCount(pClass) & ", " & BegAreaPerc & ", " & PercSamArea & ", " & BegPopPerc & ", " & pBegMean(pClass) & ", " & pBegMin(pClass) & ", " & pBegMax(pClass) & ", " & pBegStdDev(pClass) & ", " & pPDArray(pClass) & ", " & _
               pStatCount & ", " & pCatPercArea & ", " & pCatPercPop & ", " & pMean & ", " & pMin & ", " & pMax & ", " & pStdDev & ", " & pCalcMethod(pClass)
    outtext.writeline (pSummary)
  Next 'pClass in pAncArray
  
  Dim SelectType As String
  Dim Prompt As String
  outtext.writeline (SelectType)
    EndTime = Now
    Elapsed = DateDiff("s", StartTime, EndTime)
    Prompt = "The program took " & Elapsed & " seconds to run."
  outtext.writeline (Prompt)
    Elapsed = DateDiff("s", StartTime, ConversionTime)
    Prompt = "The program took " & Elapsed & " seconds to perform the vector to raster conversion."
  outtext.writeline (Prompt)
    Elapsed = DateDiff("s", ConversionTime, IntersectTime)
    Prompt = "The program took " & Elapsed & " seconds to perform the intersection."
  outtext.writeline (Prompt)
    Elapsed = DateDiff("s", IntersectTime, PreCalcTime)
    Prompt = "The program then took " & Elapsed & " seconds to perform the precalculations."
  outtext.writeline (Prompt)
    Elapsed = DateDiff("s", PreCalcTime, SelectionTime)
    Prompt = "The program then took " & Elapsed & " seconds to perform all selections and class calculations."
  outtext.writeline (Prompt)
    Elapsed = DateDiff("s", SelectionTime, SAMTime)
    Prompt = "The program then took " & Elapsed & " seconds to perform the smart areal weighting."
  outtext.writeline (Prompt)
    Elapsed = DateDiff("s", SAMTime, FinalCalcTime)
    Prompt = "The program then took " & Elapsed & " seconds to perform the final calculations."
  outtext.writeline (Prompt)
    Elapsed = DateDiff("s", FinalCalcTime, EndTime)
    Prompt = "The program then took " & Elapsed & " seconds to generate the summary table."
  outtext.writeline (Prompt)
  outtext.Close
  Set pPDArray = Nothing
  Set pAncArray = Nothing
  Set pBegCount = Nothing
  Set pBegMean = Nothing
  Set pBegMin = Nothing
  Set pBegMax = Nothing
  Set pBegStdDev = Nothing
  Set pTotAreaArray = Nothing
  Set pTotPopArray = Nothing
  Set pCalcMethod = Nothing
  Set pUnsampled = Nothing
  Set pCursor = Nothing
  Set pMXDoc = Nothing
  Set pCalc = Nothing
  Set pQueryfilter = Nothing
  Set pEditor = Nothing

Set pRasterWS = Nothing
Set pEnv = Nothing
Set pConvOp = Nothing
Set pRBCollection = Nothing
Set pLocalOp = Nothing
pRelClassCollEdit.RemoveAllRelationshipClasses
Exit Sub

'Handle regions...



'Saving the code below, not sure when or whether we'll use it at some point.
'Remove the old join, join the new final table back to the raster table


'  Join the ancillary table - this might have to be done with standalone tables
'  Dim pAncTable As ITable
'  Set pAncTable = pRasterLayer
'  Set pMemRelClassFact = New MemoryRelationshipClassFactory
'  Dim pDasyRelClass As IRelationshipClass
'  Set pDasyRelClass = pMemRelClassFact.Open("DasyTable", pAncTable, "Value", pCombTable, "LC3RASTERd", _
'                                            "forward", "backward", esriRelCardinalityOneToMany)
'  Dim pDasyTable As ITable
'  Set pRelQueryTableFact = New RelQueryTableFactory
'  Set pDasyTable = pRelQueryTableFact.Open(pDasyRelClass, True, Nothing, Nothing, "", True, True)
'' Add it to the relates for the feature layer in Map
'  pRelClassCollEdit.AddRelationshipClass pDasyRelClass
'  'Display the join
'  pDispRC.DisplayRelationshipClass pDasyRelClass, esriLeftInnerJoin


Exit Sub
ERH:
    MsgBox Err.Description & ", Error Number " & Err.Number
    MsgBox pErrorTrack
    Set pRasterWS = Nothing
    Set pEnv = Nothing
    Set pConvOp = Nothing
    Set pRBCollection = Nothing
    Set pLocalOp = Nothing
    pRelClassCollEdit.RemoveAllRelationshipClasses
    If pEditor.EditState = esriEditor.esriStateEditing Then pEditor.StopEditing True
End Sub

Private Function SetRasterWorkspace(ByVal PathName As String) As esriGeoDatabase.IWorkspace
    ' Given a pathname, returns the raster workspace object for that path
    On Error GoTo ERH
    Dim pWSF As esriGeoDatabase.IWorkspaceFactory
    Set pWSF = New esriDataSourcesRaster.RasterWorkspaceFactory
    Dim pWS As esriGeoDatabase.IWorkspace
    Set pWS = pWSF.OpenFromFile(PathName, 0)
    Set SetRasterWorkspace = pWS
    Exit Function
ERH:
    Set SetRasterWorkspace = Nothing
End Function

'Create all of the new fields in the output table, and input table if necessary
Private Sub CreateFields(pDasyTable As esriGeoDatabase.ITable)
    Dim pField As esriGeoDatabase.IField
    Dim pFieldEdit As esriGeoDatabase.IFieldEdit
    Dim pName As Variant
    Dim pFieldNames As New Collection
    Set pField = New esriGeoDatabase.Field
    Set pFieldEdit = pField
    pFieldNames.Add "TRCTCOUNT"
'    If pDasyTable.FindField("POPAREA") = -1 Then pFieldNames.Add "POPAREA" 'Add only if it didn't exist already
    pFieldNames.Add "INHAB"
    pFieldNames.Add "INHAB_CNT"
    pFieldNames.Add "INHAB_DENS"
    pFieldNames.Add "POPEST"
    pFieldNames.Add "TOTALFRACT"
    pFieldNames.Add "NEWPOP"
    pFieldNames.Add "NEWDENSITY"
    For Each pName In pFieldNames
      With pFieldEdit
          .Editable = True
          .Name = pName
          .Type = esriGeoDatabase.esriFieldTypeDouble
      End With
      pDasyTable.AddField pField
    Next 'pName in pFieldNames
End Sub

'Get the unique values in the for the Layer, Field, and QueryFilter and put them in an array
Private Function GetUnique(pQueryfilter As esriGeoDatabase.IQueryFilter, pTable As esriGeoDatabase.ITable, pFieldName As Variant) As Collection
  Dim pCursor As esriGeoDatabase.ICursor
  Dim pData As esriGeoDatabase.IDataStatistics
  Dim pEnumvarP As esriSystem.IEnumVariantSimple
  Set GetUnique = New Collection
  Set pCursor = pTable.Search(pQueryfilter, False)
  Set pData = New esriGeoDatabase.DataStatistics
  pData.Field = pFieldName
  Set pData.Cursor = pCursor
  Set pEnumvarP = pData.UniqueValues
  Dim pPop As Variant
  pPop = pEnumvarP.Next
  Dim i As Long
  For i = 1 To pData.UniqueValueCount
    GetUnique.Add pPop
    pPop = pEnumvarP.Next
  Next i
  Set pCursor = Nothing
  Set pData = Nothing
End Function

'Perform a number of calculations prior to the selection process
Private Sub PreCalc(pPopArray, PresetClasses As Collection, PresetDensities As Collection, PopKeyName, AncCatName As Variant, _
                    pDasyTable As esriGeoDatabase.ITable, popFieldName As Variant)
' !!! MLH: Added popFieldName as an argument

'Populate binary inhabited/noninhabited column.
'Make sure we have a total pixel count for each tract (carryover from rasterized tracts?)
'Calculate the inhabited pixel count for each tract.
'All of the above calculations can probably be considered precalculations, and done in one full pass through the tracts.
  
  Dim pCursor As esriGeoDatabase.ICursor
  Dim pCalc As esriGeoDatabaseUI.ICalculator
  Dim pSearchFilter As esriGeoDatabase.IQueryFilter
  Dim pUpdateFilter As esriGeoDatabase.IQueryFilter
  Dim pStats As esriGeoDatabase.IDataStatistics
  Dim pResults As esriSystem.IStatisticsResults
  Dim pPresetC As Variant
  Dim pClassNum As Variant
  Dim pPop As Variant
  'Populate column designating inhabited ancillary classes from presets.
  Set pCalc = New esriGeoDatabaseUI.Calculator
  Set pCursor = pDasyTable.Update(Nothing, False)
  With pCalc
     Set .Cursor = pCursor
         .Expression = 1
         .Field = "INHAB"
  End With
  pCalc.Calculate
  Set pCalc = Nothing
  Dim i As Integer
  i = 1
  For Each pPresetC In PresetClasses
    Set pSearchFilter = New esriGeoDatabase.QueryFilter
    pSearchFilter.WhereClause = AncCatName & " = " & pPresetC
    If pDasyTable.RowCount(pSearchFilter) = 0 Then
      PresetClasses.Remove (i)  'Program tends to crash if a preset class is not found in the table.  Remove from array to prevent this.
    Else 'If the class does exist...
      If PresetDensities(pPresetC) = 0 Then 'Only if the class has been preset to zero
        Set pUpdateFilter = New esriGeoDatabase.QueryFilter
        pUpdateFilter.WhereClause = AncCatName & " = " & pPresetC
        Set pCursor = pDasyTable.Update(pUpdateFilter, False)
        Set pCalc = New esriGeoDatabaseUI.Calculator
        With pCalc
          Set .Cursor = pCursor
              .Expression = 0
              .Field = "INHAB"
        End With
        pCalc.Calculate
        Set pCalc = Nothing
      End If
    End If
    i = i + 1
  Next 'pPresetC In PresetClasses
  For Each pPop In pPopArray
  'First populate the TRCTCOUNT field by summing the COUNT field
    Set pSearchFilter = New esriGeoDatabase.QueryFilter
    pSearchFilter.WhereClause = PopKeyName & " = " & pPop
    Set pCursor = pDasyTable.Search(pSearchFilter, False)
    Set pStats = New esriGeoDatabase.DataStatistics
    Set pStats.Cursor = pCursor
    pStats.Field = "Count"
    Set pResults = pStats.Statistics
    Set pCalc = New esriGeoDatabaseUI.Calculator
    Set pUpdateFilter = New esriGeoDatabase.QueryFilter
    pUpdateFilter.WhereClause = PopKeyName & " = " & pPop
    Set pCursor = pDasyTable.Update(pUpdateFilter, False)
    With pCalc
      Set .Cursor = pCursor
          .Expression = pResults.Sum
          .Field = "TRCTCOUNT"
    End With
    pCalc.Calculate
  'Clear variables before moving on
      Set pSearchFilter = Nothing
      Set pCursor = Nothing
      Set pStats = Nothing
      Set pResults = Nothing
      Set pCalc = Nothing
      Set pUpdateFilter = Nothing
  'Calculate the total populated pixel count of each population unit and fill the appropriate columns
    Set pSearchFilter = New esriGeoDatabase.QueryFilter
    pSearchFilter.WhereClause = PopKeyName & " = " & pPop & " AND INHAB = 1"
    Set pCursor = pDasyTable.Search(pSearchFilter, False)
    Set pStats = New esriGeoDatabase.DataStatistics
    Set pStats.Cursor = pCursor
    pStats.Field = "Count"
    Set pResults = pStats.Statistics
    Set pCalc = New esriGeoDatabaseUI.Calculator
    Set pUpdateFilter = New esriGeoDatabase.QueryFilter
    pUpdateFilter.WhereClause = PopKeyName & " = " & pPop & " AND INHAB = 1"
    If pDasyTable.RowCount(pUpdateFilter) <> 0 Then 'Not happy if population polys ended up clipped.
      Set pCursor = pDasyTable.Update(pUpdateFilter, False)
      With pCalc
      Set .Cursor = pCursor
          .Expression = pResults.Sum
          .Field = "INHAB_CNT"
      End With
      pCalc.Calculate
    End If
    Set pCursor = Nothing
    Set pStats = Nothing
    Set pResults = Nothing
    Set pCalc = Nothing
    Set pSearchFilter = Nothing
    Set pUpdateFilter = Nothing
  Next 'pPop in pPopArray
'Calculate the population density of the inhabited area (POPFIELDNAME / INHABCNT) in people/pixel, purely for comparison
  Set pCalc = New esriGeoDatabaseUI.Calculator
  Set pCursor = pDasyTable.Update(Nothing, False)
  With pCalc
  Set .Cursor = pCursor
      .Expression = "[" + popFieldName + "] / [TRCTCOUNT]" '!!! MLH: Replaced hard-coded TOTPOP field name with popFieldName
      .Field = "INHAB_DENS"
  End With
  pCalc.Calculate
  'Clear variables again before exiting sub
  Set pCursor = Nothing
  Set pCalc = Nothing
End Sub

'Subroutine for estimating population of unsampled classes based on existing data
Private Sub PopulationEst(pPopArray, pTable As esriGeoDatabase.ITable, popFieldName, PopKeyName As Variant, pErrorTrack As Variant)
  Dim pCursor As esriGeoDatabase.ICursor
  Dim pCalc As esriGeoDatabaseUI.ICalculator
  Dim pStats As esriGeoDatabase.IDataStatistics
  Dim pResults As esriSystem.IStatisticsResults
  Dim pQueryfilter As esriGeoDatabase.IQueryFilter
  Dim pEstPop, pTotalPop, pExtraPop, pEstArea As Double
  Dim pIndex As Long
  Dim pPop As Variant
  Set pQueryfilter = New esriGeoDatabase.QueryFilter
  pQueryfilter.WhereClause = "POPEST = 0 AND INHAB = 1"
  If pTable.RowCount(pQueryfilter) <> 0 Then 'Skip the entire operation if there are no empty rows.
  For Each pPop In pPopArray
    Set pQueryfilter = New esriGeoDatabase.QueryFilter
    pQueryfilter.WhereClause = PopKeyName & " = " & pPop & " and POPEST = 0 AND INHAB = 1"
    If pTable.RowCount(pQueryfilter) <> 0 Then 'Skip this operation if there are no empty rows in this source polygon.
      Set pQueryfilter = New esriGeoDatabase.QueryFilter
      pQueryfilter.WhereClause = PopKeyName & " = " & pPop
      pQueryfilter.SubFields = "POPEST"
      Set pCursor = pTable.Search(pQueryfilter, False)
      Set pStats = New esriGeoDatabase.DataStatistics
      Set pStats.Cursor = pCursor
      pStats.Field = "POPEST"
      Set pResults = pStats.Statistics
      pEstPop = pResults.Sum
      Set pQueryfilter = New esriGeoDatabase.QueryFilter
      Set pCursor = Nothing
      Set pQueryfilter = New esriGeoDatabase.QueryFilter
      pQueryfilter.WhereClause = PopKeyName & " = " & pPop
      pQueryfilter.SubFields = popFieldName
      Set pCursor = pTable.Search(pQueryfilter, False)
      pIndex = pCursor.FindField(popFieldName)
      pTotalPop = pCursor.NextRow.Value(pIndex)
      If pEstPop < pTotalPop Then 'Otherwise there's no remaining population to distribute and we can leave POPEST=0
        Set pQueryfilter = New esriGeoDatabase.QueryFilter
        pQueryfilter.WhereClause = PopKeyName & " = " & pPop & " and POPEST <> 0"
        Set pCursor = pTable.Search(pQueryfilter, False)
        Set pStats = New esriGeoDatabase.DataStatistics
        Set pStats.Cursor = pCursor
        pStats.Field = "Count"
        Set pResults = pStats.Statistics
        pEstArea = pResults.Sum
        pExtraPop = pTotalPop - pEstPop
        Set pCursor = Nothing
        Set pQueryfilter = Nothing
        Set pCalc = Nothing
        Set pCalc = New esriGeoDatabaseUI.Calculator
        Set pQueryfilter = New esriGeoDatabase.QueryFilter
        pQueryfilter.WhereClause = PopKeyName & " = " & pPop & " and POPEST = 0"
        Set pCursor = pTable.Update(pQueryfilter, False)
        'For this calculation we want to distribute the extra population amongst the unsampled inhabited dasymetric polygons
        'using areal weighting = [NEWAREA] / pExtraArea.  pExtraArea = [POPAREA] - pEstArea
        With pCalc
        Set .Cursor = pCursor
            .PreExpression = "Dim pExtraArea" & vbNewLine & _
                             "Dim EstPop" & vbNewLine & _
                             "pExtraArea = [INHAB_CNT] - " & pEstArea & vbNewLine & _
                             "If pExtraArea > 0 and [INHAB] = 1 Then" & vbNewLine & _
                             "EstPop = ([Count] / pExtraArea) * " & pExtraPop & vbNewLine & _
                             "Else EstPop = 0" & vbNewLine & _
                             "End If"
            .Expression = "EstPop"
            .Field = "POPEST"
        End With
        pErrorTrack = "The program fails on Tract " & pPop & " with an estimated population of " & pEstPop & _
                      " an actual population of " & pTotalPop & " and an estimated area of " & pEstArea
        pCalc.Calculate
      End If
    End If
  Next 'pPop in pPopArray
  End If
  Set pCalc = Nothing
  Set pCursor = Nothing
  Set pQueryfilter = Nothing
  Set pStats = Nothing
  Set pResults = Nothing
End Sub

'Final calculations
Private Sub FinalCalc(pDasyTable As esriGeoDatabase.ITable, pRegQueryfilter As esriGeoDatabase.IQueryFilter, pPopArray As Collection, _
                      PopKeyName, popFieldName As Variant)
  Dim pCalc As esriGeoDatabaseUI.ICalculator
  Dim pCursor As esriGeoDatabase.ICursor
  Dim pStats As esriGeoDatabase.IDataStatistics
  Dim pResults As esriSystem.IStatisticsResults
  Dim pQueryfilter As esriGeoDatabase.IQueryFilter
  Dim pPop As Variant
  Dim pPDFARsum As Double
  'Find the sum of all the PDF_AR values for each Population District and store as pPDFARsum
  'Calculate the total fraction as [PDF_AR] / pPDFARsum
  For Each pPop In pPopArray
    Set pQueryfilter = New esriGeoDatabase.QueryFilter
    pQueryfilter.WhereClause = PopKeyName & " = " & pPop
    Set pCursor = pDasyTable.Search(pQueryfilter, False)
    Set pStats = New esriGeoDatabase.DataStatistics
    Set pStats.Cursor = pCursor
    pStats.Field = "POPEST"
    Set pResults = pStats.Statistics
    pPDFARsum = pResults.Sum
    If pPDFARsum > 0 Then 'Because dividing by zero is bad.
      Set pCalc = New esriGeoDatabaseUI.Calculator
      Set pCursor = pDasyTable.Update(pQueryfilter, False)
      With pCalc
        Set .Cursor = pCursor
            .Expression = "[POPEST] / " & pPDFARsum
            .Field = "TOTALFRACT"
      End With
      pCalc.Calculate
      Set pCalc = Nothing
    End If
  Next

  'Calculate the actual population of the new, intersected polygon by multiplying
  'the total fraction times the population of that tract
  Set pCalc = New esriGeoDatabaseUI.Calculator
  Set pCursor = pDasyTable.Update(pRegQueryfilter, False)
   With pCalc
      Set .Cursor = pCursor
      .Expression = "[TOTALFRACT] * [" & popFieldName & "]"
      .Field = "NEWPOP"
   End With
  pCalc.Calculate
  Set pCalc = Nothing
  
  'Lastly calculate a density column for statistical purposes.
  Set pCalc = New esriGeoDatabaseUI.Calculator
  Set pCursor = pDasyTable.Update(pRegQueryfilter, False)
   With pCalc
      Set .Cursor = pCursor
      .Expression = "[NEWPOP] / [Count]"
      .Field = "NEWDENSITY"
   End With
  pCalc.Calculate
  Set pCalc = Nothing
  Set pCursor = Nothing
End Sub

Private Sub UIButtonControl1_Click()
    Dasymetric_Raster
End Sub