Main Page/Research/Smart Maps/dwelling units

From phurvitz
Jump to: navigation, search

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

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

The following source code can be used. Requirements:

  1. Existence of drive W for output. Use subst in Windows for this.
  2. Modify script in ArcToolbox to include arguments/parameters, matching the sys.argv[] lines.

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


<source lang="python">

  1. ---------------------------------------------------------------------------
  2. focal_du_per_ha.py
  3. Created on: Tue May 18 2010 02:48:36 PM
  4. (generated by ArcGIS/ModelBuilder)
  5. Usage: focal_du_per_ha <parcels> <Value_field> <Cellsize> <Neighborhood> <du_per_ha>
  6. ---------------------------------------------------------------------------
  1. Import system modules

import sys, string, os, arcgisscripting

  1. Create the Geoprocessor object

gp = arcgisscripting.create()

  1. Check out any necessary licenses

gp.CheckOutExtension("spatial")

  1. 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")

  1. output dir

outdir = "w:/" if os.path.exists(outdir) == 0: gp.AddMessage("Path does not exist.") exit()

  1. 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

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
  1. Local variables...

duperkm2 = outdir + "duperkm2.img" FocalSumR = outdir + "FocalSumR.img"

  1. Process: Polygon to Raster...

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

  1. Process: Focal Sum...

gp.AddMessage("FocalSum\n") gp.FocalStatistics_sa(duperkm2, FocalSumR, Neighborhood, "SUM", "DATA")

  1. Process: DU per hectare...

gp.AddMessage("Du per ha\n") gp.SingleOutputMapAlgebra_sa(FocalSumR + " / ( (Pi) * sqr ( " + Cellsize + " ) )", du_per_ha)

</source>