View source for Main Page/Research/Smart Maps/dwelling units

Create a smart map raster that represents the count or density of dwelling units 

# add a field to the table for dwelling units per km<sup>2</sup>.
#: greater precision than per ft<sup>2</sup>
# calculate the field, e.g.
#: <tt>du_per_km2 = [num_units]*10763910.42/[shape_area]</tt>
# convert to raster using the field
#: make a note of cell size, e.g., 10 m cell size = 32.808 ft
# focalsum of the raster
#: note the radius, e.g., 10 cells = 100 m
#: gives sum of du/km<sup>2</sup> for the circle
# divide by area of the circle to give density
#: note unit conversions

The following source code can be used. Requirements:
# Existence of drive W for output. Use <tt>subst</tt> in Windows for this.
# Modify script in ArcToolbox to include arguments/parameters, matching the <tt>sys.argv[]</tt> lines.

Note: this assumes you want the density to be in units/ha.



<source lang="python">
# ---------------------------------------------------------------------------
# focal_du_per_ha.py
# Created on: Tue May 18 2010 02:48:36 PM
#   (generated by ArcGIS/ModelBuilder)
# Usage: focal_du_per_ha <parcels> <Value_field> <Cellsize> <Neighborhood> <du_per_ha> 
# ---------------------------------------------------------------------------

# Import system modules
import sys, string, os, arcgisscripting

# Create the Geoprocessor object
gp = arcgisscripting.create()

# Check out any necessary licenses
gp.CheckOutExtension("spatial")

# Load required toolboxes...
gp.AddToolbox("C:/Program Files (x86)/ArcGIS/ArcToolbox/Toolboxes/Spatial Analyst Tools.tbx")
gp.AddToolbox("C:/Program Files (x86)/ArcGIS/ArcToolbox/Toolboxes/Conversion Tools.tbx")

# output dir
outdir = "w:/"
if os.path.exists(outdir) == 0:
	gp.AddMessage("Path does not exist.")
	exit()
	
# Script arguments...
parcels = sys.argv[1]
if parcels == '#':
 parcels = "PIN_ADDRESS" # provide a default value if unspecified

Value_field = sys.argv[2]
if Value_field == '#':
 Value_field = "du_per_km2" # provide a default value if unspecified

Cellsize = sys.argv[3]
if Cellsize == '#':
 Cellsize = "32.808" # provide a default value if unspecified
StrCellAreaMeters = str((float(Cellsize) / 3.2808)**2)

Neighborhood = sys.argv[4]
if Neighborhood == '#':
 Neighborhood = "Circle 10 CELL" # provide a default value if unspecified

du_per_ha = sys.argv[5]
if du_per_ha == '#':
 du_per_ha = outdir + "du_per_ha.img" # provide a default value if unspecified

# Local variables...
duperkm2 = outdir + "duperkm2.img"
FocalSumR = outdir + "FocalSumR.img"

# Process: Polygon to Raster...
gp.AddMessage("Polygon to raster\n")
gp.PolygonToRaster_conversion(parcels, Value_field, duperkm2, "CELL_CENTER", "NONE", Cellsize)

# Process: Focal Sum...
gp.AddMessage("FocalSum\n")
gp.FocalStatistics_sa(duperkm2, FocalSumR, Neighborhood, "SUM", "DATA")

# Process: DU per hectare...
gp.AddMessage("Du per ha\n")
gp.SingleOutputMapAlgebra_sa(FocalSumR + " /  ( (Pi) * sqr( " + StrCellAreaMeters + " ))", du_per_ha)
</source>