updating to run in crankshaft. Output still wrong but getting there

src/pg/sql/50_contours.sql

@@ -0,0 +1,32 @@
+
+CREATE OR REPLACE FUNCTION
+  _CDB_Contours (
+      subquery TEXT,
+      grid_size NUMERIC DEFAULT 100,
+      bandwidth NUMERIC DEFAULT 0.0001,
+      levels NUMERIC[] DEFAULT null
+      )
+RETURNS table (level Numeric, geom_text text )
+AS $$
+  from crankshaft.contours import cdb_generate_contours
+  # TODO: use named parameters or a dictionary
+  return cdb_generate_contours(subquery, grid_size, bandwidth, levels)
+$$ LANGUAGE plpythonu;
+
+
+CREATE OR REPLACE FUNCTION
+  CDB_Contours (
+    subquery TEXT,
+    grid_size NUMERIC DEFAULT 100,
+    bandwidth NUMERIC DEFAULT 0.0001,
+    levels NUMERIC[] DEFAULT null
+    )
+RETURNS table (level Numeric, geom geometry )
+AS $$
+BEGIN
+
+  RETURN QUERY
+    select cont.level as level, ST_GeomFromText(cont.geom_text, 4326)::geometry as geom from _CDB_Contours(subquery,grid_size,bandwidth,levels) as cont;
+END;
+$$ LANGUAGE plpgsql;
+

src/py/crankshaft/crankshaft/__init__.py

@@ -1,2 +1,3 @@
 import random_seeds
 import clustering
+import contours

src/py/crankshaft/crankshaft/contours/__init__.py

@@ -0,0 +1 @@
+from contours import *

src/py/crankshaft/crankshaft/contours/contours.py

@@ -0,0 +1,58 @@
+from scipy.stats import gaussian_kde
+from scipy.interpolate import griddata
+import numpy as np 
+from sklearn.neighbors import KernelDensity
+from skimage.measure import find_contours
+import plpy
+
+def cdb_generate_contours(query, grid_size, bandwidth, levels):
+    plpy.notice('one')
+    data   = plpy.execute( 'select ST_X(the_geom) as x , ST_Y(the_geom) as y from ({0}) as a '.format(query))
+    plpy.notice('two')
+
+    xs = [d['x'] for d in data]
+    ys = [d['y'] for d in data]
+    plpy.notice('three')
+    return generate_contours(xs,ys,grid_size,bandwidth,levels)
+  
+def scale_coord(coord, x_range,y_range,grid_size):
+    plpy.notice('ranges %,  % ', x_range, y_range)
+    return [coord[0]*(x_range[1]-x_range[0])/float(grid_size)+x_range[0],
+            coord[1]*(y_range[1]-y_range[0])/float(grid_size)+y_range[0]]
+    
+def make_wkt(data,x_range, y_range, grid_size):
+    joined = ','.join([' '.join(map(str,scale_coord(coord_pair, x_range, y_range, grid_size))) for coord_pair in data])
+    return '({0})'.format(joined)
+    
+def make_multi_line(data,x_range,y_range, grid_size):
+    joined = ','.join([ make_wkt(ring,x_range,y_range,grid_size)  for ring in data ])
+    return 'MULTILINESTRING({0})'.format(joined)
+
+def generate_contours(xs,ys, grid_res=100, bandwidth=0.001, levels=None):
+    plpy.notice("HERE")
+    max_y, min_y = np.max(ys), np.min(ys)
+    max_x, min_x = np.max(xs), np.min(xs)
+    positions = np.vstack([ys,xs]).T
+    grid_x,grid_y = np.meshgrid(np.linspace(min_x, max_x , grid_res), np.linspace(min_y, max_y, grid_res))
+    xy = np.vstack([grid_y.ravel(), grid_x.ravel()]).T
+    xy *= np.pi / 180.
+
+    plpy.notice(" Generating kernel density")
+    kde = KernelDensity(bandwidth=bandwidth, metric='haversine',
+                        kernel='gaussian', algorithm='ball_tree')
+    kde.fit(positions*np.pi/180.)
+    results = np.exp(kde.score_samples(xy))
+    results = results.reshape((grid_x.shape[0], grid_y.shape[0]))
+    
+    if not levels:
+        levels = np.linspace(results.min(), results.max(),60)
+    plpy.notice(' finding contours')
+    CS = [find_contours(results, level) for level in levels]
+    
+    vertices = []
+    for contours,level in zip(CS,levels):
+        if len(contours)>0:
+            multiline = make_multi_line(contours, (min_x,max_x), (min_y, max_y), grid_res)
+            vertices.append([level, multiline ])
+    plpy.notice('generated vertices retunring ?')
+    return vertices

src/py/crankshaft/setup.py

@@ -42,7 +42,7 @@ setup(
     # provisioned in the production servers.
     install_requires=['pysal==1.9.1'],
 
-    requires=['pysal', 'numpy' ],
+    requires=['pysal', 'numpy', 'sklearn', 'scikit-image'],
 
     test_suite='test'
 )