# Estimation Walkthrough¶

[37]:

from shapely.geometry import Polygon
import numpy as np
%matplotlib inline

[38]:

import geopandas as gpd
from tobler import area_weighted

[39]:

from tobler.area_weighted import _area_tables_binning as area_tables
from tobler.area_weighted import area_interpolate


## Example: Two GeoDataFrames¶

[40]:

polys1 = gpd.GeoSeries([Polygon([(0,0), (10,0), (10,5), (0,5)]),
Polygon([(0,5), (0,10),  (10,10), (10,5)])])

polys2 = gpd.GeoSeries([Polygon([(0,0), (5,0), (5,7), (0,7)]),
Polygon([(5,0), (5,10),  (10,10), (10,0)]),
Polygon([(0,7), (0,10), (5,10), (5,7)  ])
])

df1 = gpd.GeoDataFrame({'geometry': polys1})
df2 = gpd.GeoDataFrame({'geometry': polys2})
df1['population'] = [ 500,  200]
df1['pci'] = [75, 100]
df1['income'] = df1['population'] * df1['pci']

df2['population'] = [ 500,  100, 200]
df2['pci'] = [75, 80, 100]
df2['income'] = df2['population'] * df2['pci']

ax = df1.plot(color='red', edgecolor='k')

[41]:

ax = df2.plot(color='green', alpha=0.5, edgecolor='k')

[42]:

res_union = gpd.overlay(df1, df2, how='union')
ax = res_union.plot(alpha=0.5, cmap='tab10')
df1.plot(ax=ax, facecolor='none', edgecolor='k');
df2.plot(ax=ax, facecolor='none', edgecolor='k');


## Area Table¶

[43]:

area_tables?

Signature: area_tables(source_df, target_df, spatial_index)
Docstring:
Construct area allocation and source-target correspondence tables using a spatial indexing approach
...

NOTE: this currently relies on Geopandas' spatial index machinery

Parameters
----------
source_df : geopandas.GeoDataFrame
GeoDataFrame containing input data and polygons
target_df : geopandas.GeoDataFramee
GeoDataFrame defining the output geometries
spatial_index : str
Spatial index to use to build the allocation of area from source to
target tables. It currently support the following values:
- "source": build the spatial index on source_df
- "target": build the spatial index on target_df
- "auto": attempts to guess the most efficient alternative.
Currently, this option uses the largest table to build the
index, and performs a bulk_query on the shorter table.

Returns
-------
tables : scipy.sparse.dok_matrix
File:      ~/projects/tobler/tobler/area_weighted/area_interpolate.py
Type:      function


[44]:

area_tables(df1, df2, 'auto')

[44]:

<2x3 sparse matrix of type '<class 'numpy.float32'>'
with 5 stored elements in Dictionary Of Keys format>

[45]:

area_tables(df2, df1, 'auto')

[45]:

<3x2 sparse matrix of type '<class 'numpy.float32'>'
with 5 stored elements in Dictionary Of Keys format>

[46]:

extensive_vars = ['population', 'income']
intensive_vars = ['pci']
estimates = area_interpolate(df1, df2, extensive_variables = extensive_vars,
intensive_variables = intensive_vars)

estimates

[46]:

population income pci geometry
0 289.999998 22749.999762 82.14286 POLYGON ((0.00000 0.00000, 5.00000 0.00000, 5....
1 350.000000 28750.000000 87.50000 POLYGON ((5.00000 0.00000, 5.00000 10.00000, 1...
2 59.999996 5999.999642 100.00000 POLYGON ((0.00000 7.00000, 0.00000 10.00000, 5...
[47]:

extensive_vars = ['population', 'income']
intensive_vars = ['pci']
estimates = area_interpolate(df2, df1, extensive_variables = extensive_vars,
intensive_variables = intensive_vars)
estimates

[47]:

population income pci geometry
0 407.142866 30785.714924 77.500000 POLYGON ((0.00000 0.00000, 10.00000 0.00000, 1...
1 392.857149 34714.286193 84.999997 POLYGON ((0.00000 5.00000, 0.00000 10.00000, 1...

## Non-exhuastive case¶

Here the first set of polygons have an envelope that does not coincide with that of the second dataframe.

[48]:

polys1 = gpd.GeoSeries([Polygon([(0,0), (12,0), (12,5), (0,5)]),
Polygon([(0,5), (0,10),  (10,10), (10,5)])])

polys2 = gpd.GeoSeries([Polygon([(0,0), (5,0), (5,7), (0,7)]),
Polygon([(5,0), (5,10),  (10,10), (10,0)]),
Polygon([(0,7), (0,10), (5,10), (5,7)  ])
])

df1 = gpd.GeoDataFrame({'geometry': polys1})
df2 = gpd.GeoDataFrame({'geometry': polys2})
df1['population'] = [ 500,  200]
df1['pci'] = [75, 100]
df1['income'] = df1['population'] * df1['pci']

df2['population'] = [ 500,  100, 200]
df2['pci'] = [75, 80, 100]
df2['income'] = df2['population'] * df2['pci']

ax = df1.plot(color='red', edgecolor='k')
df2.plot(ax=ax, color='green',edgecolor='k')

[48]:

<AxesSubplot:>

[49]:

extensive_vars = ['population']
intensive_vars = ['pci']
estimates = area_interpolate(df1, df2, extensive_variables = extensive_vars,
intensive_variables = intensive_vars)
estimates

[49]:

population pci geometry
0 289.999998 82.14286 POLYGON ((0.00000 0.00000, 5.00000 0.00000, 5....
1 350.000000 87.50000 POLYGON ((5.00000 0.00000, 5.00000 10.00000, 1...
2 59.999996 100.00000 POLYGON ((0.00000 7.00000, 0.00000 10.00000, 5...
[50]:

estimates.sum()

[50]:

population    699.999994
pci           269.642860
dtype: float64

[51]:

extensive_vars = ['population']
intensive_vars = ['pci']
estimates = area_interpolate(df1, df2, extensive_variables = extensive_vars,
intensive_variables = intensive_vars,
allocate_total=False)
estimates

[51]:

population pci geometry
0 248.333333 82.14286 POLYGON ((0.00000 0.00000, 5.00000 0.00000, 5....
1 308.333333 87.50000 POLYGON ((5.00000 0.00000, 5.00000 10.00000, 1...
2 60.000000 100.00000 POLYGON ((0.00000 7.00000, 0.00000 10.00000, 5...
[52]:

estimates.sum()

[52]:

population    616.666667
pci           269.642860
dtype: float64


When setting allocate_total=False the total population of a source zone is not completely allocated, but rather the proportion of total population is set to the area of intersection over the area of the source zone.

This will have no effect when the source df is df2 and the target df is df 1:

[53]:

extensive_vars = ['population']
estimates = area_interpolate(df2, df1, extensive_variables = extensive_vars)
estimates

[53]:

population geometry
0 407.142866 POLYGON ((0.00000 0.00000, 12.00000 0.00000, 1...
1 392.857149 POLYGON ((0.00000 5.00000, 0.00000 10.00000, 1...
[54]:

extensive_vars = ['population']
estimates = area_interpolate(df2, df1, extensive_variables = extensive_vars, allocate_total=False)
estimates

[54]:

population geometry
0 407.142857 POLYGON ((0.00000 0.00000, 12.00000 0.00000, 1...
1 392.857143 POLYGON ((0.00000 5.00000, 0.00000 10.00000, 1...