This is the home to Spatial Interpolation Toolbox, a graphical user interface (GUI) for interpolating geographic vector data.
- Introduction
- Installation and Setup
- Description and Example Usage
- Troubleshooting
- Sources
- Special Thanks
- License
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.
Advanced forms of spatial interpolation are not easily accessible to all data analysts. The goal of the Spatial Interpolation Toolbox is to offer a convenient and free way to interpolate spatial data.
Spatial Interpolation Toolbox depends on the following Python Packages:
geopandaspandaspyqt
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
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 using the toolbox, you can either close the GUI or enter Ctrl+C in the Anaconda Prompt terminal.
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.
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:
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.
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. The output of this interpolation should look similar to this:
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.
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:
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.
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:
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.
For the parcel method, we will use tax lot data from NYC's MapPLUTO, and population at the census block group level from TIGER/Line. Our data will look like this to begin:
The data that we will be interpolating is population, which is currently aggregated in census block groups. Using the parcel method, the population can be disaggregated into individual parcels. Our inputs should look like this:
Note: This method can take a long time to execute when calculating an area this large with hundreds of thousands of parcel polygons
The parcel method will interpolate population into two new fields which are calculated from different inputs. One of the new fields is named ara_derived (derived from adjusted residential area), and the other field is named ru_derived (derived from number of residential units). Below are the results of the parcel method, one map for each interpolation type:
There are subtle differences in the outcomes of interpolation based on whether adjusted residential area or number of residential units are used as proxies for approximating population density. For a more in-depth explanation of the parcel based method, I recommend reading "Mapping population distribution in the urban environment: The cadastral-based expert dasymetric system" by Maantay et al. The next method, Cadastral-Based Expert Dasymetric System, can be used to determine which proxy is a more accurate determinant of population.
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.
Like the parcel method, we'll be using census block groups containing population, and parcel data. In addition, we are also using a larger census zone shapefile (which also contains population) that the smaller census zone nests inside, in this case census tracts.
The CEDS method works perfectly with census data, but theoretically will work with any two geographies that nest without intersecting
For our inputs, the field that we are interpolating (population) needs to have the same field name in both source shapefiles (tracts and block groups). Other than that condition, the inputs for CEDS are almost identical to the parcel method.
The mapped output of these inputs should look similar to this (the field expert_sys is mapped here):
The dataframe that results from the CEDS method contains both the ru_derived and ara_derived interpolations for population, as well as a new field named expert_sys. As seen in the dataframe below, CEDS determines whether to use ru_derived or ara_derived to measure population, on a census block group basis. In Block Group 3 of GEOID 360610271003 CEDS chooses the ru_derived population, then chooses the ara_derived population for block group 1 of GEOID 360610277001.
- Projection - Spatial Interpolation Toolbox will not reproject shapefiles. You must ensure that your shapefiles are in the correct projections before using this tool.
- Error Handling - In addition to error pop-ups, Spatial Interpolation Toolbox will print status messages and 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 addressed here or in an already open issue, please open a new issue 👍
-
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
- 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
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 <https://www.gnu.org/licenses/>.