Skip to main content
3 of 14
deleted 23 characters in body

How to fix this Heatmap script from failing at Kernel Density Raster creation using arcpy (ArcMap 10.5)?

I have a heatmap script that will eventually print 10 PDF maps. At the moment it is failing at Kernel Density Raster creation. Here is the error:

Reading .CSV file and plotting XY ...

Saving xy points in WGS 1984 to a layer for online use later if needed ... > Projecting xy points to TX FIPS ...

Making rasters ...

Traceback (most recent call last): File "\coacd.org\dfs\SWS\Code Enforcement\ACD Admin Operations\GIS\Operations\Council District Analysis Reporting\Auto_HeatMap_WIP\autoheatmap.py", line 55, in outKDens = KernelDensity(selection_subset, "None", cell_size=100) File "C:\Program Files (x86)\ArcGIS\Desktop10.5\ArcPy\arcpy\sa\Functions.py", line 395, in KernelDensity method) File "C:\Program Files (x86)\ArcGIS\Desktop10.5\ArcPy\arcpy\sa\Utils.py", line 53, in swapper result = wrapper(*args, **kwargs) File "C:\Program Files (x86)\ArcGIS\Desktop10.5\ArcPy\arcpy\sa\Functions.py", line 386, in Wrapper method) File "C:\Program Files (x86)\ArcGIS\Desktop10.5\ArcPy\arcpy\geoprocessing_base.py", line 510, in return lambda *args: val(*gp_fixargs(args, True)) arcgisscripting.ExecuteError: ERROR 999999: Error executing function. An invalid SQL statement was used. An invalid SQL statement was used. An invalid SQL statement was used. ERROR 010029: Unable to create the raster G:\Code Enforcement\ACD Admin Operations\GIS\Operations\Council District Analysis Reporting\Auto_HeatMap_WIP\shps\autoLayers.gdb\KernelD_csv_1. ERROR 010067: Error in executing grid expression. Failed to execute (KernelDensity).

... Script Finished at 9:45:20.58

Here is the script:

    import xlrd
    import os
    import arcpy
    from arcpy import env
    import arcpy.sa 
    from arcpy.sa import *
    
    workspace = r"\\coacd.org\dfs\SWS\Code Enforcement\ACD Admin Operations\GIS\Operations\Council District Analysis Reporting\Auto_HeatMap_WIP\shps"
    
    autoLayerGDB = os.path.join(workspace, "autoLayers.gdb")
    arcpy.env.workspace = autoLayerGDB 
    arcpy.env.overwriteOutput = True
    mxd_path = workspace 
    SOURCE_DATA = r"\\coacd.org\dfs\SWS\Code Enforcement\ACD Admin Operations\GIS\Operations\Council District Analysis Reporting\Auto_HeatMap_WIP\Source_Data\SOURCE_DATA.csv"
    print("\tReading .CSV file and plotting XY ...\n")
    spatial_ref = arcpy.SpatialReference("WGS 1984")
    output_layer = "sample_xy_points"
    
    saved_Layer = workspace +"\\sample_xy_points.lyr"
    in_x_field = "Longitude"
    in_y_field = "Latitude"
    arcpy.MakeXYEventLayer_management(SOURCE_DATA, in_x_field, in_y_field, output_layer, spatial_ref,)
    
    print("\tSaving xy points in WGS 1984 to a layer for online use later if needed ...\n")
    #Copy output to make this permanent and assign an object ID field
    arcpy.SaveToLayerFile_management(output_layer, saved_Layer)
    arcpy.CopyFeatures_management(output_layer, "coded_xy_points")
    #project layer into Texas FIPS     
    print("\tProjecting xy points to TX FIPS ...\n")
    repr_spatial_ref = arcpy.SpatialReference("NAD 1983 StatePlane Texas Central FIPS 4203 (US Feet)")
    arcpy.Project_management(output_layer, "coded_xy_points_FIPS4203", repr_spatial_ref)
    #project layer into Texas FIPS 
    print("\tMaking rasters ...\n")
    arcpy.CheckOutExtension('Spatial')
    #use 'CURRENT' if running from arcmap, when published use MXD on disk
    mxd = arcpy.mapping.MapDocument(workspace+"\\working_environment.mxd")
    df = arcpy.mapping.ListDataFrames(mxd)[0] # address location data frame and turn on select layers
    df.zoomToSelectedFeatures()
    
    #working_gdb = workspace+"\\autoLayers.gdb"
    
    arcpy.env.extent = df.extent #set extent to current one
    arcpy.env.overwriteOutput = True #overwrite existing files
    
    my_point_features = autoLayerGDB +"\\output_layer"
    scratch_output = autoLayerGDB +"\\scratch"
        
    
    #Select each district via script
    # # Within selected features, further select only those cities which have a population > 10,000   
    for i in range(1, 11):
        my_clause = 'COUNCILDISTRICT = %d' % i
        my_filename = "kden_d%d" % i
        selection_subset = arcpy.MakeFeatureLayer_management(output_layer, scratch_output, where_clause=my_clause)
        outKDens = KernelDensity(selection_subset, "None", cell_size=100)
        outKDens.save(autoLayerGDB+"\\"+my_filename)
        print("\tRaster created for District %d" % i)
    
    ## now making citywide district raster
    my_clause = 'COUNCILDISTRICT IS NOT NULL'
    my_filename = "kden_dAll"
    selection_subset = arcpy.MakeFeatureLayer_management(output_layer, scratch_output, where_clause=my_clause)
    outKDens = KernelDensity(selection_subset, "None", cell_size=50)
    outKDens.save(autoLayerGDB+"\\"+my_filename)
    print("\tRaster created for Citywide View")
    
    mxd = arcpy.mapping.MapDocument(r"\\coacd.org\dfs\SWS\Code Enforcement\ACD Admin Operations\GIS\Operations\Council District Analysis Reporting\Auto_HeatMap_WIP\shps\pdf_template.mxd")
    df = arcpy.mapping.ListDataFrames(mxd)[0] # address location data frame and turn on select layers
    
    # visibleMapLayers = ["District Outline","Lakes","case_density","District Fill","Basemap"]
    visibleMapLayers = ["District_Outlines_SRC","Lakes","case_density","District_Fills_SRC","Basemap"]
    
    mxd.activeView = "PAGE_LAYOUT" #pan to location and refresh layout
    arcpy.RefreshActiveView() #update view
    
    scratch_output = "\\scratch"
    
    outline_lyr = arcpy.mapping.ListLayers(mxd, "District_Outlines_SRC", df)[0]
    fill_lyr = arcpy.mapping.ListLayers(mxd, "District_Fills_SRC", df)[0]
    
    for i in range(1, 11):
        print("\tstarting with council district %d" % i)
        my_clause = 'COUNCILDISTRICT = %d' % i
        fill_clause = 'COUNCILDISTRICT <> %d' % i
        title_text = "Council District %d" % i
        raster_name = "kden_d%d" % i
        pdf_name = "District_%d_%s.pdf" % ( i , _time_tag )
        #-----------REFRESH TEMPLATE------------------
        for lyr in arcpy.mapping.ListLayers(mxd, "", df): #set visibility and remove older layers
            lyr.visible = False
            if lyr.name in visibleMapLayers:
                lyr.visible = True
            if lyr.name == "District_Fills_SRC":  
                lyr.definitionQuery = fill_clause
            if lyr.name == "District_Outlines_SRC":  
                lyr.definitionQuery = my_clause
            if lyr.name == "case_density":
                lyr.replaceDataSource(autoLayerGDB, "FILEGDB_WORKSPACE", raster_name)