Written and maintained by CASRAI Editorial Board
Last updated
GeoDa is an open-source software program designed for exploratory spatial data analysis (ESDA), spatial autocorrelation statistics, and spatial regression modeling. Developed by Luc Anselin and the Center for Spatial Data Science, GeoDa enables researchers to construct spatial weights matrices, compute Moran’s I and LISA cluster statistics, and explore geographic patterns through dynamically linked statistical graphics and map views.
Definition in One Paragraph
GeoDa is a specialized software environment for exploratory spatial data analysis (ESDA) that integrates statistical graphics with geographic mapping. Created by Luc Anselin to make spatial statistics accessible without command-line code, GeoDa links statistical plots with geographic maps: selecting observations in a histogram or scatter plot immediately highlights them across choropleth maps. The software provides dedicated tools for constructing spatial weights matrices, calculating global Moran’s I to evaluate overall spatial clustering, generating Local Indicators of Spatial Association (LISA) to detect hot spots and spatial outliers, and diagnosing spatial dependence in spatial data analysis.
Formula and Mathematical Foundations
GeoDa models spatial dependence using spatial weights matrices, global autocorrelation, and local decomposition statistics.
1. Spatial Weights Matrix and Row-Standardization
GeoDa models neighbor topology across n spatial units via an n × n spatial weights matrix W. For units i and j:
• wij = 1 if units i and j share a boundary or fall within a distance threshold; otherwise wij = 0.
• Diagonal elements are zero: wii = 0.
GeoDa standardizes rows so weights sum to 1, producing normalized weights wij*:
wij* = wij / ∑j=1n wij
Row-standardization transforms the spatial lag into a local weighted average of neighboring attribute values.
2. Global Moran’s I Formula
Global Moran’s I measures overall spatial dependence across the study area:
I = (n / S0) × [ ∑i=1n ∑j=1n wij (xi − x̄)(xj − x̄) ] / [ ∑i=1n (xi − x̄)2 ]
Where n is the observation count, xi and xj are observed attributes, x̄ is the sample mean, zi = xi − x̄ is the centered deviation, wij is the spatial weight, and S0 = ∑i=1n ∑j=1n wij is the sum of weights. With row-standardized weights, S0 = n, simplifying Moran’s I to:
I = [ ∑i=1n zi [Wz]i ] / [ ∑i=1n zi2 ]
Where [Wz]i = ∑j=1n wij* zj is the spatial lag. Moran’s I equals the slope of a linear regression of [Wz] on z.
3. Local Moran’s I (LISA)
To identify localized clusters, GeoDa computes Anselin’s (1995) Local Moran statistic for each unit i:
Ii = (zi / m2) × ∑j=1n wij* zj
Where m2 = (1/n) ∑i=1n zi2 is the sample variance. The mean of all local indices equals the global index: (1/n) ∑i=1n Ii = I.
Worked Example with Real Numbers and Arithmetic Shown
Consider four contiguous geographic zones arranged along a line (1 — 2 — 3 — 4).
1. Attribute Values and Centered Deviations:
Observed values: x1 = 10, x2 = 15, x3 = 25, x4 = 30.
Mean: x̄ = (10 + 15 + 25 + 30) / 4 = 80 / 4 = 20.0.
Deviations (zi = xi − 20.0): z1 = −10, z2 = −5, z3 = +5, z4 = +10.
2. Spatial Weights and Lags:
Shared boundaries define row-standardized weights:
• Zone 1: w1,2* = 1.0; [Wz]1 = 1.0(−5) = −5.0.
• Zone 2: w2,1* = 0.5, w2,3* = 0.5; [Wz]2 = 0.5(−10) + 0.5(5) = −2.5.
• Zone 3: w3,2* = 0.5, w3,4* = 0.5; [Wz]3 = 0.5(−5) + 0.5(10) = +2.5.
• Zone 4: w4,3* = 1.0; [Wz]4 = 1.0(5) = +5.0.
3. Numerator and Denominator:
Cross-products zi × [Wz]i: Zone 1: (−10)(−5.0) = 50.0; Zone 2: (−5)(−2.5) = 12.5; Zone 3: (5)(2.5) = 12.5; Zone 4: (10)(5.0) = 50.0.
Numerator sum = 50.0 + 12.5 + 12.5 + 50.0 = 125.0.
Squared deviations zi2: 100 + 25 + 25 + 100 = 250.0 (Denominator).
4. Moran’s I and Local LISA:
Global Moran’s I = 125.0 / 250.0 = 0.5000.
Local Moran values with m2 = 250.0 / 4 = 62.5:
• I1 = 50.0 / 62.5 = 0.80; I2 = 12.5 / 62.5 = 0.20; I3 = 12.5 / 62.5 = 0.20; I4 = 50.0 / 62.5 = 0.80.
Mean local index: (0.80 + 0.20 + 0.20 + 0.80) / 4 = 0.50, confirming (1/n) ∑ Ii = I.
Interpretation and Reading Guidance
Interpreting GeoDa output requires assessing both global metrics and local cluster maps.
Global Moran’s I Benchmarks
Under Complete Spatial Randomness (CSR), expected Moran’s I is E[I] = −1 / (n − 1), approaching 0 in large samples.
• I > E[I] (Positive Autocorrelation): Similar values cluster together geographically (I = 0.50 represents strong clustering).
• I ≈ E[I] (Spatial Randomness): Values show no spatial patterning.
• I < E[I] (Negative Autocorrelation): High and low values alternate in dispersed patterns.
LISA Cluster Quadrants
GeoDa maps each observation into one of four quadrants on the Moran scatter plot:
| LISA Quadrant | Local Value (zi) | Spatial Lag ([Wz]i) | Classification | GeoDa Color |
|---|---|---|---|---|
| High-High (Quadrant I) | Positive (> 0) | Positive (> 0) | Hot Spot Cluster | Dark Red |
| Low-High (Quadrant II) | Negative (< 0) | Positive (> 0) | Spatial Outlier | Light Blue |
| Low-Low (Quadrant III) | Negative (< 0) | Negative (< 0) | Cold Spot Cluster | Dark Blue |
| High-Low (Quadrant IV) | Positive (> 0) | Negative (< 0) | Spatial Outlier | Light Red |
| Not Significant | Any | Any | Random Pattern | Light Gray |
Permutation Pseudo P-Values
GeoDa computes empirical pseudo p-values by reshuffling observed attributes across locations while preserving spatial weights:
pseudo p = (M + 1) / (R + 1)
Where R is the permutation count and M is the count of permutations producing a statistic at least as extreme as observed.
Assumptions and Common Mistakes
Researchers using GeoDa should monitor five key analytical considerations:
1. Arbitrary Spatial Weights: Switching between Queen and Rook contiguity or setting an improper distance threshold can isolate observations with zero neighbors. Verify connectivity histograms to ensure network continuity.
2. Modifiable Areal Unit Problem (MAUP): Autocorrelation metrics depend on geographic scale and boundary definitions. Aggregating from tract to county level changes variance and spatial lag relationships.
3. Multiple Testing Bias in LISA: Testing every areal unit inflates false discovery rates. Apply False Discovery Rate (FDR) adjustments or Bonferroni bounds to avoid spurious cluster detections.
4. Conflating Clustering with Causality: Spatial clustering detected in GeoDa does not demonstrate causal influence. Clustered attributes frequently reflect shared socioeconomic or environmental covariates, requiring formal spatial regression analysis.
5. Unsmoothed Rare Event Rates: Mapping crude rates across areas with small populations creates artificial outliers. Apply Empirical Bayes rate smoothing in GeoDa before computing spatial autocorrelation.
How to Report It (APA Example)
When reporting results from GeoDa in an academic manuscript, specify the software version, weights matrix, permutation count, test statistics, and cluster distributions.
APA 7th Example:
“Exploratory spatial data analysis was conducted using GeoDa (Version 1.22; Anselin et al., 2006) to examine spatial autocorrelation in median household income across 120 urban census tracts. Spatial proximity was defined using a first-order Queen contiguity spatial weights matrix with row standardization. Global Moran’s I indicated significant positive spatial autocorrelation, I = .462, pseudo p = .001 (based on 999 permutations). Local Indicators of Spatial Association (LISA) identified two statistically significant High-High clusters (hot spots, p < .01) encompassing 19 tracts in the northern district and one Low-Low cluster (cold spot, p < .01) of 14 tracts in the southern corridor.”
How to Calculate It in Excel, R, or Python
Spatial autocorrelation workflows developed in GeoDa can be reproduced programmatically across standard environments.
1. Python Implementation
In Python, use GeoPandas and PySAL (libpysal, esda):
import geopandas as gpd
from libpysal.weights import Queen
from esda.moran import Moran, Moran_Local
gdf = gpd.read_file("tracts.geojson")
w = Queen.from_dataframe(gdf, use_index=False)
w.transform = 'R'
lm = Moran(gdf['income'], w, permutations=999)
print(f"Moran's I: {lm.I:.4f}, pseudo p: {lm.p_sim:.4f}")
lisa = Moran_Local(gdf['income'], w, permutations=999)
gdf['quadrant'] = lisa.q
2. R Implementation
In R, use the spdep package:
library(sf)
library(spdep)
tracts <- st_read("tracts.geojson")
nb <- poly2nb(tracts, queen = TRUE)
w <- nb2listw(nb, style = "W", zero.policy = TRUE)
moran_mc <- moran.mc(tracts$income, w, nsim = 999, zero.policy = TRUE)
print(moran_mc)
lisa_local <- localmoran_perm(tracts$income, w, nsim = 999, zero.policy = TRUE)
3. Excel Formula Logic
For small spatial matrices, implement the arithmetic directly in Excel:
• Observations x: Cells B2:B5 contain 10, 15, 25, 30.
• Centered Deviations zi: Enter =B2-AVERAGE($B$2:$B$5) in C2 and copy down.
• Spatial Lag [Wz]i: Enter row-standardized neighbor sums (for Zone 2: =0.5*C2+0.5*C4).
• Cross-Products: Enter =C2*D2 in E2 and copy down.
• Global Moran’s I: Enter =SUM(E2:E5)/DEVSQ(B2:B5), evaluating to 0.5000.
Frequently Asked Questions
What is the difference between GeoDa and standard GIS software?
Traditional GIS programs emphasize data management, map production, and geometric overlays. GeoDa focuses specifically on exploratory spatial data analysis (ESDA), interactive brushing between charts and maps, spatial weights construction, Moran’s I and LISA testing, and spatial regression diagnostics.
What is the difference between Queen and Rook contiguity in GeoDa?
Rook contiguity defines two polygons as neighbors only if they share a common linear border of positive length. Queen contiguity defines polygons as neighbors if they share either a border segment or a single common vertex point. Queen contiguity is commonly preferred because it prevents disconnects caused by minor digitizing differences.
How does GeoDa calculate p-values for spatial autocorrelation?
GeoDa generates pseudo p-values through conditional random spatial permutations rather than assuming normality. By repeatedly permuting observed values across spatial locations while holding weights fixed, GeoDa builds an empirical reference distribution: pseudo p = (M + 1) / (R + 1), where R is the permutation count and M is the count of permuted statistics meeting or exceeding the observed index.
Can GeoDa analyze point geometries?
Yes. GeoDa supports both point and polygon geometries. Users can construct spatial weights for points using distance bands, k-nearest neighbors (k-NN), Delaunay triangulations, or Thiessen (Voronoi) polygon conversions.
Sources
• Anselin, L. (1995). Local indicators of spatial association—LISA. Geographical Analysis, 27(2), 93–115. https://doi.org/10.1111/j.1538-4632.1995.tb00338.x (Accessed 2026-11-08)
• Anselin, L. (2024). An Introduction to Spatial Data Science with GeoDa. Volume 1: Exploring Spatial Data. CRC Press / Chapman & Hall. ISBN: 9781032341446. https://geodacenter.github.io/ (Accessed 2026-11-08)
• Anselin, L., Syabri, I., & Kho, Y. (2006). GeoDa: An introduction to spatial data analysis. Geographical Analysis, 38(1), 5–22. https://doi.org/10.1111/j.0016-7363.2005.00671.x (Accessed 2026-11-08)
• Center for Spatial Data Science (CSDS), University of Chicago. (2024). GeoDa Software Documentation and Tutorials: Global Spatial Autocorrelation & LISA. https://geodacenter.github.io/documentation.html (Accessed 2026-11-08)
