The basic function accepts a GeoPandas GeoDataFrame and the column name of the field for which you want to create a Dougenik, Chrisman, Niemayer cartogram.
Currently, a stopping rule has not been implemented, therefore a trial-and-error approach is required to arrive at a suitable number of iterations. Generally a good solution requires more than 20 iterations. The verbose option prints out the current iteration and error, a suitable number of iterations is generally arrived at when both mean error and the largest single error approach one. The precision of the error after 1 is effectively the percentage error - 1.16 = 16% error between desired and actual area.
As the output is a GeoPandas GeoDataFrame, plotting is straightforward using the built-in .plot() method.
NB The field used for creating the Dougenik, Chrisman, Niemayer cartogram should contain values which are positive and greater than zero. It is important that the input data be in a projected (Cartesian) coordinate system. The GeoDataFrame should be structured so that 1 row is 1 record.
The pseudocode for computing the cartogram is given in the reference cited, and is as follows:
For each polygon
Read and store PolygonValue (negative value illegal)
Sum PolygonValue into TotalValue
For each iteration (user controls when done)
For each polygon
Calculate area and centroid (using current boundaries)
Sum areas into TotalArea
For each polygon
Desired = (TotalArea * (PolygonValue/TotalValue))
Radius = SquareRoot(Area/pi)
Mass = SquareRoot(Desired/pi) - SquareRoot(Area/pi)
SizeError = Max(Area, Desired) / Min(Area, Desired)
ForceReductionFactor = 1 / (1 + Mean(SizeError))
For each boundary line; Read coordinate chain
For each coordinate pair
For each polygon centroid
Find angle, Distance from centroid to coordinate
If (Distance > Radius of polygon)
Fij = Mass * (Radius/Distance)
Else
Fij = Mass * (Distance$^{2}$/ Radius$^{2}$) * (4 - (3 *(Distance / Radius)))
Using Fij and angles, calculate vector sum
Multiply by ForceReductionFactor
Move coordinate accordingly
Write distorted line to output and plot result
Dougenik J., Chrisman N. and Niemeyer D. 1985. An algorithm to construct continuous area cartograms. Professional Geographer. 37(1): 75-81.
In [1]:
# Some required libraries.
import geopandas as gpd
import pandas as pd
import matplotlib.pyplot as plt
# import the non-contiguous cartogram function
from cartogrampy.dcn_cartogram import dcn_cartogram
%matplotlib inline
In [2]:
# Get some data for London
# Read in London Boroughs geojson and project to British National Grid.
ldn_boro = gpd.read_file('data/LDN_Boro.geojson').to_crs(epsg=27700)
# Get some population data for London Boroughs
ldn_pop = pd.read_excel('data/Pandas_Lon_Pop.xlsx')
# Merge the population data with the geospatial data
ldn_boro = ldn_boro.merge(ldn_pop, how='left', left_on='GSS_CODE', right_on='New Code')
In [3]:
# Compute DCN Cartogram - max error c. 4% after 10 iterations, after 16 < 1%.
ldn_dcn = dcn_cartogram(ldn_boro,2015,iterations=16,verbose=False)
In [4]:
# Map the outcome
fig, ax = plt.subplots(1,2,figsize=(16,8))
ldn_boro.plot(column = 2015,cmap = 'PuRd', scheme = 'Quantiles', k=4,legend=True, ax=ax[0])
ldn_dcn.plot(color='r',ax=ax[1])
ax[0].axis('equal')
ax[1].axis('equal')
Out[4]:
In [5]:
# Get some data for US States
usstate = gpd.read_file('data/US_State_2016_5m.geojson')
# Set crs
usstate.crs = {'init': u'epsg:4269'}
# Get continental US and project to NAD83 Contiguous US Albers.
usstate = usstate[~usstate['STUSPS'].isin(['AK', 'HI', 'AS', 'PR',
'GU', 'MP', 'VI'])].to_crs({'init': 'epsg:5070'})
# Read in state populations
state_pop = pd.read_excel('data/Pandas_US_Pop.xlsx')
# Merge population data
usstate = usstate.merge(state_pop, how='left', left_on='STUSPS', right_on='Code')
In [6]:
# Compute DCN Cartogram - max error c. 35% after 10 iterations, after 20 c.8%., after 30 c. 4%
us_dcn = dcn_cartogram(usstate,2016,iterations=30,verbose=True)
In [7]:
# Map the outcome
fig, ax = plt.subplots(1,2,figsize=(16,8))
usstate.plot(column = 2016,cmap = 'PuRd', scheme = 'Quantiles', k=4,legend=True, ax=ax[0])
us_dcn.plot(color='r',ax=ax[1])
ax[0].axis('equal')
ax[1].axis('equal')
Out[7]: