diff --git a/src/pg/sql/22_maxp.sql b/src/pg/sql/22_maxp.sql
new file mode 100644
index 0000000..e23ee12
--- /dev/null
+++ b/src/pg/sql/22_maxp.sql
@@ -0,0 +1,16 @@
+-- max-p regionalization
+
+CREATE OR REPLACE FUNCTION
+  CDB_MaxP(
+      subquery TEXT,
+      colnames TEXT[],
+      min_size int default 1,
+      initial int default 99,
+      geom_col TEXT DEFAULT 'the_geom',
+      id_col TEXT DEFAULT 'cartodb_id')
+RETURNS TABLE (region_class text, p_val numeric, rowid bigint)
+AS $$
+  from crankshaft.clustering import MaxP
+  maxp = MaxP()
+  return maxp.maxp(subquery, colnames, floor=min_size)
+$$ LANGUAGE plpythonu;
diff --git a/src/py/crankshaft/crankshaft/analysis_data_provider.py b/src/py/crankshaft/crankshaft/analysis_data_provider.py
index cbc27bc..ae141e8 100644
--- a/src/py/crankshaft/crankshaft/analysis_data_provider.py
+++ b/src/py/crankshaft/crankshaft/analysis_data_provider.py
@@ -65,3 +65,17 @@ class AnalysisDataProvider:
             return data
         except plpy.SPIError, err:
             plpy.error('Analysis failed: %s' % err)
+
+    def get_maxp(self, params):
+        """fetch data for spatial markov"""
+        try:
+            query = pu.construct_neighbor_query('queen', params)
+            data = plpy.execute(query)
+
+            if len(data) == 0:
+                # TODO: replace with better message in PR#157
+                plpy.error('No non-null valued rows')
+
+            return data
+        except plpy.SPIError, err:
+            plpy.error('Analysis failed: %s' % err)
diff --git a/src/py/crankshaft/crankshaft/clustering/__init__.py b/src/py/crankshaft/crankshaft/clustering/__init__.py
index d9682fa..aa334f0 100644
--- a/src/py/crankshaft/crankshaft/clustering/__init__.py
+++ b/src/py/crankshaft/crankshaft/clustering/__init__.py
@@ -2,3 +2,4 @@
 from moran import *
 from kmeans import *
 from getis import *
+from maxp import *
diff --git a/src/py/crankshaft/crankshaft/clustering/maxp.py b/src/py/crankshaft/crankshaft/clustering/maxp.py
new file mode 100644
index 0000000..c135369
--- /dev/null
+++ b/src/py/crankshaft/crankshaft/clustering/maxp.py
@@ -0,0 +1,68 @@
+"""
+  max-p clustering
+"""
+
+import pysal as ps
+import numpy as np
+import random
+
+import crankshaft.pysal_utils as pu
+from crankshaft.analysis_data_provider import AnalysisDataProvider
+
+
+class MaxP:
+    def __init__(self, data_provider=None):
+        if data_provider:
+            self.data_provider = data_provider
+        else:
+            self.data_provider = AnalysisDataProvider()
+
+    def maxp(self, subquery, colnames, floor=1,
+             geom_col='the_geom', id_col='cartodb_id'):
+        """
+            Inputs:
+            @param subquery (text): subquery to expose the data need for the
+                                    analysis. This query needs to expose all
+                                    of the columns in `colnames`, `id_col`, and
+                                    `geom_col`
+            @param colnames (list): list of column names (as strings)
+            @param floor (float): ...
+            @param geom_col (text): geometry column used for calculating the
+                                    spatial neighborhood
+            @param id_col (text): id column used for keeping the identity of
+                                  the data
+            Outputs: a list of tuples with the following columns:
+
+            classification_id: group that the geometry belongs to
+            rowid: identifier from id_col
+        """
+        params = {'subquery': subquery,
+                  'colnames': colnames,
+                  'id_col': id_col,
+                  'geom_col': geom_col}
+
+        resp = self.data_provider.get_maxp(params)
+
+        attr_vals = pu.get_attributes(resp, len(colnames))
+
+        weight = pu.get_weight(resp, w_type='queen')
+        r = ps.Maxp(weight, attr_vals,
+                    floor=floor,
+                    floor_variable=np.ones((weight.n, 1)))
+
+        cluster_classes = get_cluster_classes(weight.id_order, r.regions)
+        r.inference()
+        return zip(cluster_classes, [r.pvalue] * len(weight.id_order),
+                   weight.id_order)
+
+
+def get_cluster_classes(ids, clusters):
+    """
+
+    """
+    cluster_classes = []
+    for i in ids:
+        for r_id, r in enumerate(clusters):
+            if i in r:
+                cluster_classes.append(r_id)
+    return cluster_classes