<SYSTEM_TASK:>
Given the following text description, write Python code to implement the functionality described below step by step
<END_TASK>
<USER_TASK:>
Description:
 Spatial Analysis
"Everything is related to everything else, but near things are more related than distant things." -Waldo Tobler
If the strength of a relationship between entities increases with their proximity, then spatial analysis/modeling is essential for understanding the relationship's process and pattern. Today we focus on exploratory spatial data analysis (ESDA) to discover patterns in spatial data.
Overview of today's topics
Step1:  1. Data prep
1.1. California census tract geometries
Step2:  1.2. California tract-level census variables
Step3:  1.3. Merge the data
Step4:  2. Initial exploration
Let's do some quick mapping and analysis of distributions and correlations for a couple variables of interest.
Step5:  Looks like we have some missing values. We'll spatially interpolate them later.
Visually, it appears that these two variables may be negatively correlated? In general, where one is high, the other is low.
Step6:  Every 1 person/km^2 increase in density is associated with a m change in median household income.
Step7:  3. Spatial weights matrix
Spatial analysis depends on spatial relationships. A spatial weights matrix defines the spatial relationships among our units of analysis (tracts, in this case). It tells how they're spatially connected to one another. These weights can take on many different forms. Pick the right form for your theoretical needs, including
Step8:  3.2. Contiguity-based weights
Step9:  How many neighbors does the average tract have?
Step10:  Plot a census tract of interest, along with its neighbors
Step11:  3.3. Distance-based weights
Step12:  3.4. Distance-based weights
Step13:  3.5. Standardizing weights
A spatial weights matrix with raw values (e.g., binary 1s and 0s for neighbor/not) is not always the best for analysis. Some sort of standardization is useful. Typically, we want to apply a row-based transformation so every row of the matrix sums up to 1. We'll see some examples of why this matters in practice shortly.
Step14:  PySAL supports the following transformations
Step15:  4. Spatial interpolation
Interpolation lets you estimate unobserved values based on observed values. With spatial data, you can perform spatial interpolation by filling in missing data points based on nearby values. This assumes positive spatial autocorrelation exists
Step16:  4. Spatial lag
Spatial lag tells us how values locate near other (dis)similar values. While spatial interpolation filled in unobserved (missing) values using nearby values, spatial lag lets us compare observed values to their nearby values.
Here we calculate the spatial lag of a variable. If the spatial weights matrix is row-standardized (important), then the spatial lag is the average value of an observation's neighbors, however "neighbor" is defined in the matrix.
Step17:  5. Spatial autocorrelation
Spatial autocorrelation is a central question in ESDA. Statistical models typically assume that the observations are independent of each other. This assumption is violated when a variable's value at one location is correlated with its value at nearby locations.
Such spatial autocorrelation is common in the real world due to proximity-based spillover effects. For example, a home's value may be a function of its own characteristics and accessibility, but it is also a function of nearby homes' values. In other words, homes near one another tend to have similar home values.
positive spatial autocorrelation
Step18:  5.1. Moran's I
Moran's I measures global spatial autocorrelation
Step19:  If we generated a large number of maps with the same values but randomly allocated over space, and calculated Moran's I for each of these maps, only 1/1000 of them would display a larger absolute value than the one we computed from the real-world data set. Thus there is a 1/1000 chance of getting the observed value of Moran's I if the spatial distribution of our variable is random. We can conclude that the variable's distribution is statistically significantly postively spatially autocorrelated.
Step20:  5.2. Moran plots
A Moran plot scatter plots the spatially-lagged values (y-axis) vs the original variable's values (x-axis).
Step21:  Notice the 95% confidence interval shading and the positive slope. Given the p-value of Moran's I that we calculated earlier, we can conclude that the slope of the line is statistically-significantly different from zero.
More useful, however, is a standardized Moran plot. Moran's I is the slope of the line in the standardized Moran plot, which makes this all a bit easier to conceptualize.
Step22:  5.3. LISAs
Local Indicators of Spatial Autocorrelation
Step23:  Now map the tracts, colored according to their LISA quadrants, to identify clusters
<ASSISTANT_TASK:>
Python Code:
import geopandas as gpd
import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
import pysal as ps
import seaborn as sns
from scipy.stats import stats
np.random.seed(0)
Explanation: Spatial Analysis
"Everything is related to everything else, but near things are more related than distant things." -Waldo Tobler
If the strength of a relationship between entities increases with their proximity, then spatial analysis/modeling is essential for understanding the relationship's process and pattern. Today we focus on exploratory spatial data analysis (ESDA) to discover patterns in spatial data.
Overview of today's topics:
  - Tobler's first law of geography
  - spatial weights
  - spatial interpolation
  - spatial lag
  - spatial autocorrelation
  - hot spot mapping
Today we will conduct an exploratory spatial analysis of LA county household income.
End of explanation
# load CA tracts, display shape
tracts_ca = gpd.read_file('../../data/tl_2017_06_tract/')
tracts_ca = tracts_ca.set_index('GEOID')
tracts_ca.shape
# what variables are present?
tracts_ca.columns
# inspect the first 5 rows
tracts_ca.head()
# retain LA county only (and drop channel island tracts)
tracts_ca = tracts_ca[tracts_ca['COUNTYFP']=='037'].drop(index=['06037599100', '06037599000'])
tracts_ca.shape
# project spatial geometries to a meter-based projection for SoCal
crs = '+proj=utm +zone=11 +ellps=WGS84 +datum=WGS84 +units=m +no_defs'
tracts_ca = tracts_ca.to_crs(crs)
Explanation: 1. Data prep
1.1. California census tract geometries
End of explanation
# load CA tract-level census variables
df_census = pd.read_csv('../../data/census_tracts_data_ca.csv', dtype={'GEOID10':str}).set_index('GEOID10')
df_census.shape
df_census.columns
df_census.head()
Explanation: 1.2. California tract-level census variables
End of explanation
# merge tract geometries with census variables
tracts = tracts_ca.merge(df_census, left_index=True, right_index=True, how='left')
tracts.shape
# calculate pop density in persons per sq km
# turn any infinities into nulls
tracts['pop_density'] = tracts['total_pop'] / (tracts['ALAND'] / 1e6)
tracts = tracts.replace([np.inf, -np.inf], np.nan)
tracts.columns
Explanation: 1.3. Merge the data
End of explanation
# descriptive stats
tracts['med_household_income'].describe()
# descriptive stats
tracts['pop_density'].describe()
# inspect these variables' statistical distributions
fig, axes = plt.subplots(1, 2, figsize=(12, 2))
ax1 = sns.boxplot(ax=axes[0], x=tracts['med_household_income'])
ax2 = sns.boxplot(ax=axes[1], x=tracts['pop_density'])
# map a couple variables to inspect their spatial distributions
cols = ['pop_density', 'med_household_income']
for col in cols:
    ax = tracts.dropna(subset=[col]).plot(column=col,
                                          scheme='NaturalBreaks',
                                          cmap='plasma',
                                          figsize=(4, 4),
                                          legend=True,
                                          legend_kwds={'bbox_to_anchor': (1.7, 1)})
    ax.set_title(col)
    _ = ax.axis('off')
Explanation: 2. Initial exploration
Let's do some quick mapping and analysis of distributions and correlations for a couple variables of interest.
End of explanation
# calculate correlation coefficient and p-value
subset = tracts.dropna(subset=['pop_density', 'med_household_income'])
r, p = stats.pearsonr(x=subset['pop_density'],
                      y=subset['med_household_income'])
print('r={:.4f}, p={:.4f}'.format(r, p))
# quick and dirty scatter plot with matplotlib
fig, ax = plt.subplots()
sc = ax.scatter(x=subset['pop_density'],
                y=subset['med_household_income'],
                s=1)
# estimate a simple linear regression model with scipy
# what if you log transform your variables first?
m, b, r, p, se = stats.linregress(x=subset['pop_density'],
                                  y=subset['med_household_income'])
print(f'm={m:.4f}, b={b:.4f}, r^2={r**2:.4f}, p={p:.4f}')
Explanation: Looks like we have some missing values. We'll spatially interpolate them later.
Visually, it appears that these two variables may be negatively correlated? In general, where one is high, the other is low.
End of explanation
# now it's your turn
# look through the list of columns, pick two new variables, and map them
# do they look like they are correlated? would you expect them to be?
Explanation: Every 1 person/km^2 increase in density is associated with a m change in median household income.
End of explanation
# get the tract labels (GEOIDs) and pick one (arbitrarily) to work with throughout
labels = tracts.index.tolist()
label = labels[603]
label
%%time
# calculate rook spatial weights
w_rook = ps.lib.weights.Rook.from_dataframe(tracts, ids=labels, id_order=labels)
# find the neighbors of some tract
# this is a raw contiguity matrix, so weights are binary 1s and 0s meaning neighbor/not
w_rook[label]
Explanation: 3. Spatial weights matrix
Spatial analysis depends on spatial relationships. A spatial weights matrix defines the spatial relationships among our units of analysis (tracts, in this case). It tells how they're spatially connected to one another. These weights can take on many different forms. Pick the right form for your theoretical needs, including:
rook contiguity
queen contiguity
k-nearest neighbors
distance band
3.1. Contiguity-based weights: rook contiguity
Using rook contiguity, two spatial units must share an edge of their boundaries to be considered neighbors. This isn't terribly common in practice since queen is usually more useful, but it's worth understanding as a trivial example.
End of explanation
%%time
# calculate queen spatial weights
w_queen = ps.lib.weights.Queen.from_dataframe(tracts, ids=labels, id_order=labels)
# find the neighbors of some tract
# this is a raw contiguity matrix, so weights are binary 1s and 0s meaning neighbor/not
w_queen[label]
# how many neighbors does this tract have?
w_queen.cardinalities[label]
# convert cardinalities to series and describe data
pd.Series(w_queen.cardinalities).describe()
Explanation: 3.2. Contiguity-based weights: queen contiguity
Using queen contiguity, two spatial units need only share a vertex (a single point) of their boundaries to be considered neighbors.
End of explanation
# min number of neighbors
w_queen.min_neighbors
# max number of neighbors
w_queen.max_neighbors
# islands are observations with no neighbors, disconnected in space (can cause modeling problems)
w_queen.islands
Explanation: How many neighbors does the average tract have?
End of explanation
fig, ax = plt.subplots(figsize=(6, 6))
tracts.plot(ax=ax, facecolor='#666666', edgecolor='w', linewidth=0.5)
# plot some tract of interest in red
tract = tracts.loc[[label]]
tract.plot(ax=ax, facecolor='#ff0000', edgecolor='w', linewidth=2)
# plot the neighbors in blue
neighbors = tracts.loc[w_queen[label]]
neighbors.plot(ax=ax, facecolor='#0033cc', edgecolor='w', linewidth=2)
# zoom to area of interest
xmin, ymin, xmax, ymax = neighbors.unary_union.bounds
ax.axis('equal')
ax.set_xlim(xmin-100, xmax+100)  # +/- 100 meters
ax.set_ylim(ymin, ymax)
ax.set_title('Neighbors of tract {}'.format(label))
_ = ax.axis('off')
%%time
# draw a queen-contiguity graph of the tracts
fig, ax = plt.subplots(figsize=(12, 12), facecolor='#111111')
tracts.plot(ax=ax, facecolor='#333333', edgecolor='k', linewidth=0.3)
# extract centroids of tract and its neighbors, then draw lines between them
for tract, neighbors in w_queen:
    tract_centroid = tracts.loc[tract, 'geometry'].centroid
    for neighbor_centroid in tracts.loc[neighbors, 'geometry'].centroid:
        Xs = [tract_centroid.x, neighbor_centroid.x]
        Ys = [tract_centroid.y, neighbor_centroid.y]
        ax.plot(Xs, Ys, color='r', linewidth=0.3)
_ = ax.axis('off')
Explanation: Plot a census tract of interest, along with its neighbors:
End of explanation
%%time
# k-nearest neighbors finds the closest k tract centroids to each tract centroid
w_knn = ps.lib.weights.KNN.from_dataframe(tracts, k=6)
# they all have exactly k neighbors
w_knn.neighbors[label]
Explanation: 3.3. Distance-based weights: k-nn
Find the k-nearest neighbors of each tract, by centroid.
End of explanation
# calculate maximum nearest neighbor distance so each unit is assured of >=1 neighbor
x = tracts.centroid.x
y = tracts.centroid.y
coords = np.array([x, y]).T
threshold = ps.lib.weights.min_threshold_distance(coords)
threshold
%%time
# calculate linear decay continuous weights
w_dist = ps.lib.weights.distance.DistanceBand.from_dataframe(tracts,
                                                             threshold=threshold,
                                                             binary=False,
                                                             alpha=-1)
# how many distance-band neighbors does our tract have?
len(w_dist.neighbors[label])
# map the neighbors, colored by weight from nearest to furthest
fig, ax = plt.subplots(figsize=(6, 6))
tracts.plot(ax=ax, facecolor='#333333', edgecolor='gray', linewidth=0.1)
# get the tract of interest and its neighbors/weights
tract = tracts.loc[[label]]
weights = pd.Series(w_dist[label])
neighbors = tracts.loc[weights.index, ['geometry']]
neighbors['weights_scaled'] = weights
# plot the tract's neighbors in blues by weight
neighbors.plot(ax=ax,
               column='weights_scaled',
               cmap='Blues_r',
               edgecolor='gray',
               linewidth=0.3,
               scheme='NaturalBreaks')
# plot the tract of interest in red
tract.plot(ax=ax, facecolor='r', edgecolor='r', linewidth=0.1)
# zoom to area of interest
xmin, ymin, xmax, ymax = neighbors.unary_union.bounds
ax.set_xlim(xmin, xmax)
ax.set_ylim(ymin, ymax)
ax.set_title('Neighbors of tract {}'.format(label))
_ = ax.axis('off')
# now it's your turn
# recompute the distance-based spatial weights with a gravity decay
# how and why does this impact the number of neighbors and the map above? why?
Explanation: 3.4. Distance-based weights: distance band
Here, other tracts are considered neighbors of some tract if they are within a given threshold distance of it, by centroid. Distance band weights can be specified to take on continuous values rather than binary (1s and 0s), with these values being the inverse distance between each pair of "neighboring" units.
linear distance-decay exponent is -1, so $w_l=\frac{1}{d}$
gravity model distance-decay exponent is -2, so $w_g=\frac{1}{d^2}$
End of explanation
# inspect the neighbors and weights of our tract
w_queen[label]
# check the current transformation of the weights matrix (O = original)
w_queen.get_transform()
# transform the queen weights
w_queen.set_transform('R')
w_queen[label]
# transform the linear-decay distance-based weights
w_dist.set_transform('R')
#w_dist[label]
Explanation: 3.5. Standardizing weights
A spatial weights matrix with raw values (e.g., binary 1s and 0s for neighbor/not) is not always the best for analysis. Some sort of standardization is useful. Typically, we want to apply a row-based transformation so every row of the matrix sums up to 1. We'll see some examples of why this matters in practice shortly.
End of explanation
# save your matrix to disk
f = ps.lib.io.open('tracts_queen.gal', 'w')
f.write(w_queen)
f.close()
# read a matrix from disk (notice its transformation)
w_queen = ps.lib.io.open('tracts_queen.gal', 'r').read()
w_queen[label]
Explanation: PySAL supports the following transformations:
O: original, returning the object to the initial state
B: binary, with every neighbor having assigned a weight of 1
R: row-based, with all the neighbors of a given observation adding up to 1
V: variance stabilizing, with the sum of all the weights being constrained to the number of observations
It can take a long time to calculate a weights matrix for a large data set. Once you've created yours, you might want to save it to disk to re-use in subsequent analyses.
End of explanation
# how many tracts are missing values for this variable?
col = 'med_household_income'
nulls = tracts[pd.isnull(tracts[col])].index
len(nulls)
# for example, this tract is missing a value
tract = nulls[0]
tract
# local averaging: equal-weighted queen-adjacent tracts
neighbors = w_queen[tract]
tracts.loc[neighbors, col].mean()
# or, calculate inverse distance weighted mean
neighbors = w_dist[tract]
inv_dist_wt = pd.Series(neighbors)
(tracts.loc[neighbors, col] * inv_dist_wt).sum()
# or, interpolate all the missing values across this variable
estimates = {}
for tract in nulls:
    neighbors = w_dist[tract]
    inv_dist_wt = pd.Series(w_dist[tract])
    estimates[tract] = (tracts.loc[neighbors, col] * inv_dist_wt).sum()
pd.Series(estimates).head()
# now it's your turn
# spatially interpolate missing values for median home value
Explanation: 4. Spatial interpolation
Interpolation lets you estimate unobserved values based on observed values. With spatial data, you can perform spatial interpolation by filling in missing data points based on nearby values. This assumes positive spatial autocorrelation exists: more on that in a moment. Remember Tobler's first law of geography.
Nearest neighbor interpolation is perhaps the simplest method: just assign the value of the nearest neighbor
Local averaging assigns missing values by taking the average of adjacent neighbors or neighbors within some radius
Inverse distance weighting assigns missing values using a distance weighted average: that is, the mean weighs nearby values more than it weighs distant values (and your distance decay choice is important!)
Kriging is a sophisticated method that incorporates information about spatial trends and autocorrelation with a variogram
We'll look at an example comparing local averaging to inverse distance weighting.
End of explanation
# pick a variable to investigate and drop null rows
col = 'med_household_income'
tracts_not_null = tracts[[col, 'geometry']].dropna()
y = tracts_not_null[col]
# recompute spatial weights for just these observations then row-standardize
w_queen = ps.lib.weights.Queen.from_dataframe(tracts_not_null)
w_queen.set_transform('R')
# compute spatial lag
y_lag = ps.lib.weights.lag_spatial(w_queen, y)
# is a tract's med income similar to those of its neighbors?
col_lag = f'{col}_lag'
data_lag = pd.DataFrame(data={col:y, col_lag:y_lag}).astype(int)
data_lag
Explanation: 4. Spatial lag
Spatial lag tells us how values locate near other (dis)similar values. While spatial interpolation filled in unobserved (missing) values using nearby values, spatial lag lets us compare observed values to their nearby values.
Here we calculate the spatial lag of a variable. If the spatial weights matrix is row-standardized (important), then the spatial lag is the average value of an observation's neighbors, however "neighbor" is defined in the matrix.
End of explanation
# does household income exhibit spatial autocorrelation?
# let's find out
data_lag.sample(5)
Explanation: 5. Spatial autocorrelation
Spatial autocorrelation is a central question in ESDA. Statistical models typically assume that the observations are independent of each other. This assumption is violated when a variable's value at one location is correlated with its value at nearby locations.
Such spatial autocorrelation is common in the real world due to proximity-based spillover effects. For example, a home's value may be a function of its own characteristics and accessibility, but it is also a function of nearby homes' values. In other words, homes near one another tend to have similar home values.
positive spatial autocorrelation: nearby values tend to be more similar (e.g. income, home values, temperature, rainfall)
negative spatial autocorrelation: nearby values tend to be more dissimilar (e.g. fire stations, libraries)
Substantive spatial autocorrelation can be explained by social or economic theory that describes a spatial relationship. Nuisance spatial autocorrelation stems from data problems.
End of explanation
# calculate the statistic
mi = ps.explore.esda.Moran(data_lag[col], w_queen)
# show the I value
mi.I
# statistical inference: show the p value
mi.p_sim
Explanation: 5.1. Moran's I
Moran's I measures global spatial autocorrelation: do things tend to be near other (dis)similar things. Values > 0 indicate positive spatial autocorrelation, and values < 0 indicate negative spatial autocorrelation.
End of explanation
# now it's your turn
# calculate the moran's I of median home values
# is it statistically significant? what does it tell you?
Explanation: If we generated a large number of maps with the same values but randomly allocated over space, and calculated Moran's I for each of these maps, only 1/1000 of them would display a larger absolute value than the one we computed from the real-world data set. Thus there is a 1/1000 chance of getting the observed value of Moran's I if the spatial distribution of our variable is random. We can conclude that the variable's distribution is statistically significantly postively spatially autocorrelated.
End of explanation
fig, ax = plt.subplots(figsize=(6, 6))
sns.regplot(x=col, y=col_lag, data=data_lag, scatter_kws={'s':1, 'color':'gray'})
plt.show()
Explanation: 5.2. Moran plots
A Moran plot scatter plots the spatially-lagged values (y-axis) vs the original variable's values (x-axis).
End of explanation
# standardize the variable's values (i.e., calculate z-scores)
y_std = (y - y.mean()) / y.std()
y_std.head()
# compute spatial lag of standardized values and save as series with same index
y_std_lag = pd.Series(ps.lib.weights.lag_spatial(w_queen, y_std),
                      index=y_std.index,
                      name=col_lag)
y_std_lag
# estimate a simple linear regression model
m, b, r, p, se = stats.linregress(x=y_std, y=y_std_lag)
print('m={:.4f}, b={:.4f}, r^2={:.4f}, p={:.4f}'.format(m, b, r ** 2, p))
# the slope is the same as moran's I, calculated earlier
mi.I
# standardized moran's plot
fig, ax = plt.subplots(figsize=(4, 4))
ax.scatter(x=y_std, y=y_std_lag, s=1, color='gray')
# draw quadrants and ignore outliers beyond 3 std devs (99.7% of distribution)
plt.axvline(0, c='k', alpha=0.5)
plt.axhline(0, c='k', alpha=0.5)
ax.set_xlim(-3, 3)
ax.set_ylim(-3, 3)
# draw a line with moran's I as the slope
Xs = pd.Series([-3, 3])
Ys = Xs * mi.I
line = ax.plot(Xs, Ys, lw=2)
# now it's your turn
# visualize a standardized moran's plot of median home values
Explanation: Notice the 95% confidence interval shading and the positive slope. Given the p-value of Moran's I that we calculated earlier, we can conclude that the slope of the line is statistically-significantly different from zero.
More useful, however, is a standardized Moran plot. Moran's I is the slope of the line in the standardized Moran plot, which makes this all a bit easier to conceptualize.
End of explanation
# standardized moran's plot again, from above, but labeled this time
fig, ax = plt.subplots(figsize=(6, 6))
ax.scatter(x=y_std, y=y_std_lag, s=1, color='gray')
# draw quadrants and ignore outliers beyond 3 std devs
plt.axvline(0, c='k', alpha=0.5)
plt.axhline(0, c='k', alpha=0.5)
ax.set_xlim(-3, 3)
ax.set_ylim(-3, 3)
# label the quadrants
ax.text(1.25, 1.25, 'HH', fontsize=30)
ax.text(1.25, -1.75, 'HL', fontsize=30)
ax.text(-1.75, 1.25, 'LH', fontsize=30)
ax.text(-1.75, -1.75, 'LL', fontsize=30)
# draw a line with moran's I as the slope
Xs = pd.Series([-3, 3])
Ys = Xs * mi.I
line = ax.plot(Xs, Ys, lw=2)
# calculate LISA values using the queen spatial weights
lisa = ps.explore.esda.Moran_Local(data_lag[col], w_queen)
# set the statistical significance threshold (alpha)
alpha = 0.05
# identify whether each observation is significant or not
# p-value interpretation same as earlier with moran's I
data_lag['significant'] = lisa.p_sim < alpha
data_lag['significant'].value_counts()
# identify the quadrant each observation belongs to
data_lag['quadrant'] = lisa.q
data_lag['quadrant'] = data_lag['quadrant'].replace({1:'HH', 2:'LH', 3:'LL', 4:'HL'})
data_lag['quadrant'].sort_values().value_counts()
# what have we got in the end?
data_lag
Explanation: 5.3. LISAs
Local Indicators of Spatial Autocorrelation: are there specific areas with high concentrations of (dis)similar values?
Moran's I tells us about spatial clustering globally across the data set as a whole. However, it does not tell us where these clusters occur. For that, we need a local measure. Essentially, we will classify the data set's observations into four groups based on the four quadrants of the Moran plot:
HH: high value near other high values (hot spots)
LL: low value near other low values (cold spots)
HL: high value near low values (spatial outliers)
LH: low value near high values (spatial outliers)
Let's see what that looks like, visually.
End of explanation
fig, ax = plt.subplots(figsize=(9, 9))
# merge original tracts and LISA quadrants data together, plot tracts basemap
tracts_lisa = tracts.merge(data_lag, how='left', left_index=True, right_index=True)
tracts_lisa.plot(ax=ax, facecolor='#999999', edgecolor='k', linewidth=0.1)
# plot each quandrant's tracts (if significant LISA) in a different color
quadrant_colors = {'HH':'r', 'LL':'b', 'LH':'skyblue', 'HL':'pink'}
for q, c in quadrant_colors.items():
    mask = tracts_lisa['significant'] & (tracts_lisa['quadrant']==q)
    rows = tracts_lisa.loc[mask]
    rows.plot(ax=ax, color=c, edgecolor='k', linewidth=0.1)
ax.axis('off')
fig.savefig('clusters.png', dpi=600, bbox_inches='tight')
Explanation: Now map the tracts, colored according to their LISA quadrants, to identify clusters:
End of explanation
<END_TASK>