Created
May 18, 2013 05:54
-
-
Save danabauer/5603394 to your computer and use it in GitHub Desktop.
Dasymetric mapping (areal interpolation) script by Torrin Hultgren. Based on research by Jeremy Mennis. More information at http://astro.temple.edu/~jmennis/research/dasymetric/.
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
| # --------------------------------------------------------------------------- | |
| # 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 |
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment
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...