Spatial Interpolation Toolbox

This is the home to Spatial Interpolation Toolbox, a graphical user interface (GUI) for interpolating geographic vector data.

Table of Contents

• Introduction
• Installation and Setup
  • Conda
  • Launching Spatial Interpolation Toolbox
• Description and Example Usage
  • Areal Weighting
    • Description
    • Usage
  • Binary Method
    • Description
    • Usage
  • Limiting Variable Method
    • Description
    • Usage
  • N-Class Method
    • Description
    • Usage
  • Parcel Method
    • Description
  • Cadastral-Based Expert Dasymetric System
    • Description
• Troubleshooting
• Sources
• Special Thanks
• License

Introduction

Spatial Interpolation Toolbox is a Python-based GUI that is able to interpolate spatial data in vector format. Spatial Interpolation Toolbox currently implements six different forms of spatial interpolation that are based on existing bodies of research.

Many forms of spatial interpolation are quite involved, and not easily accessible to data analysts. The goal of the Spatial Interpolation Toolbox is to offer a simplified and free way to interpolate spatial data.

Installation and Setup

Spatial Interpolation Toolbox depends on the following Python Packages:

• geopandas
• pandas
• pyqt

Conda

These packages depend on several low-level libraries for geospatial analysis, and can be challenging to install. For this reason, I recommend installing all of the dependencies using the conda package manager. Conda can be obtained by installing the Anaconda Distribution (A free Python distribution tailored for data science).

Once you have conda installed, launch Anaconda Prompt. Then, download this repo and change directory to the folder where you've downloaded the repo.

git clone https://github.com/mikeRobWard/spatial-interpolation-toolbox

cd spatial-interpolation-toolbox

Once there, create a new conda environment using the requirements.txt file, then activate the environment

conda create --name spatial_toolbox --file requirements.txt

conda activate spatial_toolbox

Launching Spatial Interpolation Toolbox

Now that you have the repo downloaded, conda installed, and your new environment running, you can launch Spatial Interpolation Toolbox from Anaconda Prompt using this command:

python sp_int_toolbox.py

When you're finished, you can either close the GUI or press Ctrl+C in the Anaconda Prompt terminal.

Description and Example Usage

Areal Weighting

Description

The areal weighting method interpolates data into target polygons by using the ratio of intersected area to source area. It accepts two shapefiles – a source and target, a list of columns to be interpolated, and an optional suffix for the new name of the interpolated column in the target shapefile. In the function, the source polygons are reindexed for summing later, and the areas of the source polygons are calculated. The source and target are intersected, and the area of the intersected polygons is calculated. Each intersected area is divided by the source area that encapsulates it for its areal weight. The function then iterates through the selected columns and multiplies each value by the areal weight. The target shapefile is returned with the interpolated columns.

Usage

For this example, we will be using open data from Philadelphia. The first shapefile is crash data aggregated by Traffic Analysis zone (TAZ). The second shapefile is Census Block Groups.

To begin, lets take a look at our two shapefiles in your preferred GIS viewer:

In this example, we want to interpolate the number of crashes from TAZ in the source layer, to Census Block group in our target layer. We can see from the crash-data attributes that the field for aggregated crashes is named Count_

Let's input the fields into Spatial Interpolation Toolbox:

This will output a new shapefile to the directory that you chose. Lets open the new shapefile and compare it to our source shapefile, the crash data by TAZ:

Binary Method

Description

This method accepts two shapefiles - a source shapefile which should contain the values that will be interpolated - and an ancillary shapefile containing a column with categorical geographic data such as land use types. The function also takes an input called exclusion field which allows the user to pass in the name of the column that contains the categorical data. Another input, exclusion value(s), allows the user to pass in a list of values that appear within the exclusion field column. The values that are passed to exclusion field will be dropped from the ancillary shapefile before it is spatially intersected with the source shapefile so that the geography of the exclusionary values is not included in the intersected shapefile. This effectively turns the ancillary shapefile into a mask, which masks out the geography of all values that were passed to exclusion value. Returns a shapefile that has the masked areas clipped and the interpolation values disaggregated based on areal weight to the non-clipped zones.

Usage

In this example, we will use the Philadelphia crash data again, but this time we will use a land use shapefile as an ancillary data source. Let's take a look at our data:

This method will use the land use shapefile to mask out certain land use types from the crash data shapefile. Car crashes definitely don't happen on water, and there may be other land use types you'd want to mask out. For this example, let's assume that we want to interpolate the car crash data to just residential land use. Here's what our inputs will look like:

The field containing land use types is named C_DIG1, which contains a numbered value corresponding to the land use type. Residential corresponds to 1, so we will exclude all other values. Land use shapefiles can contain hundreds of thousands of polygons, so be patient if you are processing a large area. The output of this interpolation should look similar to this:

Limiting Variable Method

Description

The limiting variable method interpolates data into disaggregated target polygons by setting thresholds to area-class categories. It accepts two shapefiles – a source and ancillary (landuse most common), the area-class column in the ancillary shapefile, a dictionary for specifying thresholds to each area-class, a list of columns to be interpolated, an optional source identifier, and an optional suffix for the new name of the interpolated column. Source polygons are reindexed and the area of each is calculated, an intersection is performed, and intersected areas are calculated. The values of the dictionary are placed in a new threshold field and their keys are matched with the specified area-class category. After the areal weight is found, a copy of the dictionary is made with values of none or 0 removed (these correspond to the class with no threshold). Starting with the most restrictive, the specified columns are multiplied by their areal weight and clipped at the specified threshold per square unit. The area that has been used is decremented from the source area and areal weight is recalculated. The most restrictive class is then removed from the dictionary, and this process repeats until all the classes have been removed from the dictionary. Finally, the remaining data is interpolated into the class with no restriction. The target shapefile is returned with interpolated columns.

Note: The thresholds for limiting variable are values per square unit. Your square units will vary depending on your projection.

Usage

For this example, we can continue to use the Philadelphia crash data and Philadelphia land use data. Our starting data will look like this:

And our inputs to Spatial Interpolation Toolbox will look something like this:

And the output of the limiting variable function will look like this:

N-Class Method

Description

The n-class method interpolates data into disaggregated target polygons by assigning weights to area-class categories. It accepts two shapefiles – a source and ancillary (landuse most common), the area-class column in the ancillary shapefile, a dictionary for specifying percentages to each area-class, a list of columns to be interpolated, an optional source identifier, and an optional suffix for the new name of the interpolated column. Like areal weighting, source polygons are reindexed and the area of each is calculated.The dictionary values are placed into a new percentage field and their keys are matched with the specified area-class column. After intersecting, the areal weight for each new polygon is calculated and multiplied by its corresponding user-defined percentage. Each of those products is then divided by the sum of all the products per source polygon. That fraction is called class_weight and is multiplied by column values for interpolation. The target shapefile is returned with interpolated columns.

Usage

For testing the n-class method, we can continue using the Philadelphia crash data and Philadelphia land use data. Our starting data will look like this again:

The inputs for n-class method are very similar to the limiting variable method, but instead of passing in thresholds based on square units, we pass in percentages as a decimal for our thresholds. The percentages should add up to 100%, regardless of how many classes you are splitting between. For this example, we'll assign 75% to residential, 20% to commercial, and 5% to industrial:

The output should look something like this:

Parcel Method

Description

The parcel based method disaggregates population from a large geography to the tax lot level by using residential area and number of residential units as proxies for population distribution. It accepts two shapefiles, a zone shapefile with geography and population, and a parcel shapefile which contain geography, total units per parcel, residential units per parcel, building area per parcel, and residential area per parcel. This method returns a shapefile at the tax lot level that has two calculated columns of disaggregated population, one based on residential area and one based on residential units.

Cadastral-Based Expert Dasymetric System

Description

The CEDS method works in conjunction with the parcel based method to determine whether adjusted residential area or number of residential units are a more accurate determinant when disaggregating population. The CEDS method accepts three shapefiles, two zone shapefiles that must nest with each other and contain geometry and population, and a parcel shapefile that contains geometry, total units per parcel, residential units per parcel, building area per parcel, and residential area per parcel. The parcel based method is called twice inside the CEDS method, once using the larger zone shapefile as an input, and once using the smaller nested zone shapefile as an input to the parcel method. The populations at the tax lot level that were derived from the large zone are then reaggregated back up to the small zone level. The absolute value of the difference between the large zone based populations and small zone estimated population are then calculated. Finally, for each parcel, if the absolute difference between the large zone based population and the small zone estimated population based on residential units is less than or equal to the absolute difference between the large zone population and small zone estimated population based on adjusted residential area, then the population estimate from the small zone based on residential units is determined to be the more accurate disaggregation. Otherwise, the population estimate from the small zone based on adjusted residential area is determined by the CEDS method to be the more accurate measure of disaggregation. This method returns one shapefile at the tax lot level with the parcel based method calculations, plus an additional column that contains the selected outcome of the CEDS method.

Troubleshooting

Common issues:

• Projection - Spatial Interpolation Toolbox will not reproject shapefiles for you. You must ensure that your shapefiles are in the correct projections before using this tool.
• Field type - The fields that you interpolate must be stored numerically (int / float). If you pass in a field that's stored as a string it may fail, even if it actually only contains numerical values.
• Error Handling - In addition to error pop-ups, Spatial Interpolation Toolbox will print exception messages to Anaconda prompt's terminal window. If the program fails during interpolation and the cause isn't clear, check the terminal for additional insight.

If you encounter a problem or bug that isn't answered here, please open an issue 👍

Sources

• Mennis, Jeremy (2009). Dasymetric Mapping for Estimating Population in Small Areas. Geography Compass, pp 727-745

• Hultgren, T. and J. Mennis (2006). Intelligent Dasymetric Mapping and Its Application to Areal Interpolation. Cartography and Geographic Information Science 33, pp 179-194

• Maantay, Juliana Astrud, Andrew R. Maroko, and Christopher Herrmann. "Mapping population distribution in the urban environment: The cadastral-based expert dasymetric system (CEDS)." Cartography and Geographic Information Science 34.2 (2007): 77-102.

• Cory L. Eicher & Cynthia A. Brewer (2001) Dasymetric Mapping and ArealInterpolation: Implementation and Evaluation, Cartography and Geographic Information Science,28:2, 125-138, DOI: 10.1559/152304001782173727

Special Thanks

• John Fitzgibbons - John contributed the code behind the Areal Weighting, Limiting Variable, and N-Class methods
• Lee Hachadoorian - Lee helped formulate the idea for this project, and provided great feedback along the way

License

Spatial Interpolation Toolbox
Copyright (C) 2021 Michael Ward

This program is free software: you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.

This program is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
GNU General Public License for more details.

You should have received a copy of the GNU General Public License
along with this program.  If not, see 
   .