Area Weighted Join vs. Standard Join

Students from an advanced suitability analysis course this summer needed to create a report that specified the percent of the join layer that intersected the target layer. For example, they needed to calculate the percent of each use from a landuse shapefile that intersect each zip code in Texas. I wrote a quick script for the class that generated the report they needed, but the implications are astounding to me.

Here is what I mean by astounding. As a test, I calculated the median household income and the total population within a 1-mile radius around each dance club in Arlington, TX comparing the following methods: (1) Total population using ArcMap's standard spatial join tool, (2) Total population using an area-weighted summation, (3) Average median household income using standard spatial join, (4) Average median household income using an area-weighted average.

The following table displays the results.



Field A displays the dance club's name. Fields B, C, and D above display the difference between using the standard ArcMap spatial join tool and a weighted-average spatial join tool when calculating average household income. Fields E, F, and G display the differences when calculating the total population.

The differences in both cases are quite high. I am thinking of myself and all of the students who I have seen naively rely on the standard spatial join tool for these types of calculations. Wow...

Why is there such a difference?

An area-weighted spatial join between two polygons comes in two flavors, depending on whether ti is calculating an average or a sum.

If it is calculating an average, the formula is [area-percent] * [value] + [area-percent] * [value]... The most important consideration is the percentages of the join features that are within each target feature. For example, in a particular zip code, there might be 3 block groups. Let's further suppose that block group 1 comprises 50%, block group 2 comprises 35%, and block group 3 comprises 15%.

If it is calculating a sum, the most important consideration is the percent of the join feature that actually intersects the target feature. The formula is ( [% area intersects target feature] * [value] + [% area intersects target feature] * [value] ) / number of intersecting join features. This is why you will see a much larger error when using the standard spatial join tool for summations than for averages. If 2% of a block group intersects a zip code, the standard tool will include the entire population of the block group instead of only 2%.

Is This a Perfect Solution?

No. This assumes a perfectly even distribution within each join feature. It is, however, a huge improvement.

Where Can I Get the Script?

Download it here. Extract the compressed archive and you will see three Python scripts and an ArcGIS toolbox. Open ArcMap or ArcCatalog, ensure ArcToolbox is visible, and add the Spatial Join Tools.tbx (single-click).

Caveat: These scripts are first drafts and have not been tested on any systems other than the ArcINFO Desktop 9.1 & 9.2 systems here at UT Arlington. There is no documentation. Also, the scripts run on the slow side. Eventually these will be optimized, but at the current time they are presented as is.

Description of the three tools:

1. Average Area Weighted: Use this tool to calculate an area-weighted average spatial join between two polygons.
2. Sum Area Weighted Join: Use this tool to calculate an area-weighted summation spatial join between two polygons.
3. Percent Area Report: Use this tool to generate a report that specifies the percent of the join layer that intersected the target layer.

ArcMap2GMap for ArcGIS 9.0, 9.1, and 9.2

As I just posted earlier, I have devised band-aid solutions that will allow all of my scripts to operate in ArcGIS 9.whatever, but that I must create separate scripts for each version.


Just created ArcMap2GMap scripts for ArcGIS 9.2 and 9.0, in addition to 9.1. Click here to download from ArcScripts.

Native Support for Geoprocessor in ArcGIS 9.2...Doh!

OK, so finally re-installed ArcGIS 9.2 after my rash initial installation went awry. Now experiencing first-hand that all of the Python scripts I developed in for ArcGIS 9.0 & 9.1 do not work in ArcGIS 9.2.

Why? As the ESRI documentation states here, "At ArcGIS 9.2, there is native Python support for geoprocessing scripting."

What does this mean in practical terms? Replace the COM connection code at the tippety-top of your scripts.


Replace:
import win32com.client
gp = win32com.client.Dispatch("esriGeoprocessing.GpDispatch.1")

With:
import arcgisscripting
gp = arcgisscripting.create()

Now, how to make a single script that will work with either ArcGIS 9.1 or 9.2? I really do not know. The only method that I am aware of to detect the version number is to pull that info from the registry. For me, that is way, way not worth it. Until a nice solution comes along, I will create different scripts for different versions.

While I am discussing Python scripting differences, there is still an unresolved issue (for me) using the searchcursor to access the geometry object for points between 9.0 and 9.1/9.2. Sent a query a while back to the ArcView-L list, but no one was able to help. Which is a shame because the last time I posted a scripting query (about constructing multipart polygons and inner circles) I received a fantastic answer within a couple of days from Nathan Warmer (ESRI).

Anyway, here is a snippet that highlights the difference:

rows = gp.SearchCursor(inputFC)
row = rows.Next()
# For each row
while row:
feat = row.shape
LUArray = feat.GetPart()
# Following line is required for 9.1 & 9.2
pnt = LUArray
# Following 2 lines must be removed for version 9.1 & 9.2 They are essential for 9.0.
# LUArray.Reset()
# pnt = LUArray.Next()
row = rows.Next()
...

I do see that it is no longer necessary to store point features in object arrays, but I have been unable to devise a solution that would work seamlessly across all ArcGIS versions.

Anyway, just thought I would post this difference here while I had it in my mind.

Python Script: Shape to Text to Shape

Except for ironing out a few remaining buglets and cleaning up the code, I have just about completed a script that will convert a shapefile of any geometry into a text file, and then will convert that text file back into a shapefile. (I previously wrote about the development of this script here and here.)

Script configured for ArcGIS 9.1: Download
This will not work for ArcGIS 9.0 and has not yet been tested for 9.2. For a narrative of my ArcGIS 9.2 woes, go here.

Purpose

This script has two potential uses: Edit an existing shapefile or Create a new shapefile in text format.

Description

There are actually two scripts: shp2exch.py and exch2shp.py that are launched independent of each other.

I initially developed the script so that faculty in our Earth & Environmental Science department can disassemble shapefiles, apply tectonic rotation formulas to the resulting text files, and then reassemble the shapefiles.

• shp2txt converts points, polylines (including multipart), and polygons (including multipart and inner circles) to a text file using an exchange format. The text file can be edited by hand or of course another application can be developed to edit these text files.
• txt2shp converts an existing text file in the exchange format to a shapefile.

Instructions

Download and unzip the contents. Open ArcMap or ArcCatalog 9.1. View the ArcToolbox pane. Right-click on ArcToolbox and select Add Toolbox. Browse to the directory where you unzipped the scripts, and click once on exchanger, and click Open.

Now you can expand the exchanger toolbox to launch the two scripts.

There are also two executable files that launch forms external of ArcGIS. However, ArcGIS still needs to be installed on the system and everything runs slower. You can give it a shot, but I advise using the toolbox within ArcGIS.

The Exchange Format

The text file uses a CSV comma delimited (comma-separated values) file structure with a extension. I used a .csv extension so that the file can be launched directly into Excel by double-clicking.

The text file uses an exchange format that follows the following structure:

• First line contains all of the attribute (field) names in the shapefile, plus xLatitude, xLongitude, and geometry
• Second line contains the field type of the attributes
• The third line through the last line contain the X and Y value for each vertex.
  
• In addition, each new feature, part, or inner circle
  is prefaced with a line beginning with NEW and the attributes for that feature or part.

Here is an example of a text file representing a U.S. states polygon shapefile. Note the multipart polygon of Hawaii.
Here is an example of a text file representing freeway lines in North Central Texas.

Known Issues

1. Coordinate system is preserved because the XY units are preserved in the text file. However, I have been having trouble preserving the complete projection information, so the resulting shapefile projections are not defined.
   
2. Currently only works with shapefiles, and not geodatabase feature classes. This is due to poor planning as the script currently relies on the FID field, and does not look for an OBJECTID field.
3. Probably does not work with ArcGIS 9.2. After last week's troubles, I do not expect to install and fix any 9.2 bugs within the next couple of weeks.
4. Code is currently a bit sloppily written. Needs to be cleaned up a bit.

ArcMap2GMap v.2

Just completed a major revision of ArcMap2GMap, the Python script that converts shapefiles from ArcMap to the HTML/JavaScripts necessary to overlay the layers on top of a Google Map. To download the file, see the ArcScripts entry. See here for the initial release of ArcMap2GMap last year.

An example of a webpage generated by ArcMap2GMap can be viewed here.

The scripts, VBA, HTML, and JavaScripts were written by myself and Kaushal Gala, my GRA who has since gone on to bigger and brighter things.

Here is a brief description (see the readme file for more details):

ArcMap2GMap converts shapefiles to Google Map overlays, including all HTML and JavaScript files necessary to run. The generated HTML file includes the following functionality: multiple layers visible/invisible, geocoding, proximities (top 10 closest visible points), and driving directions.

Converts point shapefiles, polyline shapefiles, and polygon shapefiles. Handles polygons with multipart and/or donut holes (inner circles). Options include the ability to simplify polylines (and polygons) to minimize the number of vertices Google Maps needs to render.

To Run:

Extract the contents directly to the C:\ drive. This will create the following directory: C:\ArcMap2GMap. Then launch the C:\ArcMap2GMap\ArcMap2GMap.mxd ArcMap project file.

Requires ArcGIS Desktop 9.1 Not tested yet on 9.2.

Developed by the University of Texas at Arlington Library: GIS Program (http://library.uta.edu/gis). Please direct all comments to gislib@uta.edu.