inital attempt at contours

Adapt Makefile of the extension for some PGXS versions

pg/Makefile

@@ -28,3 +28,6 @@ REGRESS_OPTS = --inputdir='$(TEST_DIR)' --outputdir='$(TEST_DIR)'
 PG_CONFIG = pg_config
 PGXS := $(shell $(PG_CONFIG) --pgxs)
 include $(PGXS)
+
+# This seems to be needed at least for PG 9.3.11
+all: $(DATA)

pg/doc/05_cdb_union_adjacent.md

@@ -1,25 +0,0 @@
-### Union Adjacent
-
-This is an aggregate function that will take a set of polygons and return a geometry array
-of regions where the polygons are continuous. Basically it combines polygons
-which are touching in to single polygons.
-
-It takes a single value:
-
-* `geometry` a list of geometries to be clustered and joined
-
-and returns
-
-* `geometry[]` an array of the joined geometries.
-
-An example usage would be something like:
-
-```postgresql
-  with joined_polygons as (
-    select cdb_union_adjacent(the_geom) regions from some_table
-  )
-  select unnest(region) the_geom from joined_polygons
-```
-
-which will produce a table with regions of continuous polygons from the original
-table.

pg/sql/0.0.1/05_cdb_union_adjacent.sql

@@ -1,43 +0,0 @@
-CREATE OR REPLACE FUNCTION _cdb_final_union_adjacent( joined_geoms geometry[] )
-RETURNS geometry[] AS $$
-BEGIN
-    RETURN joined_geoms;
-END
-$$ LANGUAGE plpgsql;
-
-
-CREATE OR REPLACE FUNCTION _cdb_state_update_union_adjacent(clusters geometry[], new_geom  geometry)
-RETURNS geometry[] AS $$
-DECLARE
-  joins  geometry[] :='{}';
-  unjoined geometry[] :='{}';
-  i integer;
-  combined geometry;
-BEGIN
-  joins := (select array_agg(g)
-            from unnest(clusters) a(g)
-            where ST_TOUCHES(g, new_geom));
-
-  unjoined := (select array_agg(g)
-               from unnest(clusters) a(g)
-               where ST_TOUCHES(g, new_geom) = false);
-
-  IF array_length(joins, 1) > 0 THEN
-    joins := array_append(joins, new_geom);
-    combined := ST_UNION(joins);
-  ELSE
-    combined := new_geom;
-  END IF;
-
-  unjoined := array_append(unjoined, combined);
-  RETURN unjoined;
-END
-$$
-LANGUAGE plpgsql;
-
-CREATE AGGREGATE cdb_union_adjacent(geometry)(
-  SFUNC=_cdb_state_update_union_adjacent,
-  STYPE=geometry[],
-  FINALFUNC=_cdb_final_union_adjacent,
-  INITCOND='{}'
-);

pg/sql/0.0.1/07_contours.sql

@@ -0,0 +1,12 @@
+
+
+CREATE OR REPLACE FUNCTION
+  cdb_contours_count (
+      query TEXT,
+      levels NUMERIC[]
+  )
+RETURNS TABLE (the_geom geometry , level Numeric)
+AS $$
+  from crankshaft.contours import create_countours_count
+  return create_countours_count(query,levels)
+$$ LANGUAGE plpythonu;

pg/test/0.0.1/expected/05_cdb_union_adjacent_test.out

@@ -1,22 +0,0 @@
-\i test/fixtures/touching_polygons.sql
--- test table (polygons, some of which touch and some which dont)
-CREATE TABLE touching_polygons(cartodb_id integer, the_geom geometry);
-INSERT INTO  touching_polygons VALUES
-(1, ST_GeomFromText('POLYGON ((0 0, 1 0,1 1, 0 1, 0 0 ))')),
-(2, ST_GeomFromText('POLYGON ((1 0, 2 0, 2 1, 1 1, 1 0))')),
-(1, ST_GeomFromText('POLYGON ((0 1, 1 1,1 2, 0 2, 0 1 ))')),
-(4, ST_GeomFromText('POLYGON ((3 0, 4 0, 4 1, 3 1, 3 0))')),
-(5, ST_GeomFromText('POLYGON ((3 1, 4 1, 4 2, 3 2, 3 1))'));
-WITH joined_polygons AS (
-  SELECT cdb_crankshaft.cdb_union_adjacent(the_geom) the_geom FROM touching_polygons
-),
-unnested_polygons as (
-  select unnest(joined_polygons.the_geom) the_geom from joined_polygons
-)
-select ST_ASTEXT(unnested_polygons.the_geom) from unnested_polygons;
-                   st_astext                    
-------------------------------------------------
- POLYGON((1 0,0 0,0 1,0 2,1 2,1 1,2 1,2 0,1 0))
- POLYGON((4 1,4 0,3 0,3 1,3 2,4 2,4 1))
-(2 rows)
-

pg/test/0.0.1/sql/05_cdb_union_adjacent_test.sql

@@ -1,9 +0,0 @@
-\i test/fixtures/touching_polygons.sql
-
-WITH joined_polygons AS (
-  SELECT cdb_crankshaft.cdb_union_adjacent(the_geom) the_geom FROM touching_polygons
-),
-unnested_polygons as (
-  select unnest(joined_polygons.the_geom) the_geom from joined_polygons
-)
-select ST_ASTEXT(unnested_polygons.the_geom) from unnested_polygons;

pg/test/fixtures/touching_polygons.sql

@@ -1,8 +0,0 @@
--- test table (polygons, some of which touch and some which dont)
-CREATE TABLE touching_polygons(cartodb_id integer, the_geom geometry);
-INSERT INTO  touching_polygons VALUES
-(1, ST_GeomFromText('POLYGON ((0 0, 1 0,1 1, 0 1, 0 0 ))')),
-(2, ST_GeomFromText('POLYGON ((1 0, 2 0, 2 1, 1 1, 1 0))')),
-(1, ST_GeomFromText('POLYGON ((0 1, 1 1,1 2, 0 2, 0 1 ))')),
-(4, ST_GeomFromText('POLYGON ((3 0, 4 0, 4 1, 3 1, 3 0))')),
-(5, ST_GeomFromText('POLYGON ((3 1, 4 1, 4 2, 3 2, 3 1))'));

python/crankshaft/crankshaft/__init__.py

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

python/crankshaft/crankshaft/contours/__init__.py

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

python/crankshaft/crankshaft/contours/contours.py

@@ -0,0 +1,47 @@
+
+import matplotlib.pyplot as plt
+import numpy as np
+import plpy
+
+def contour_to_polygon(contour):
+    plpy.notice('appending contour ')
+    c = np.append(contour, [contour[0]], axis=0)
+    points =','.join( [  " ".join(str(a) for a in b) for b in c])
+
+    return "POLYGON(({points}))::geometry".format(points=points)
+
+def create_countours_count(query,levels,mesh_size=20):
+    qresult = plpy.execute( "select ST_X(the_geom)::Numeric as x, ST_Y(the_geom)::Numeric as y from ({query}) a ".format(query=query))
+    x =[]
+    y =[]
+    for a in qresult:
+        if a['x'] and a['y']:
+            x.append(float(a['x']))
+            y.append(float(a['y']))
+
+    plpy.notice(np.shape(x))
+    plpy.notice(np.shape(y))
+
+    if None in x:
+        plpy.notice("NULL IN LIST X ")
+    if None in y:
+        plpy.notice("NULL IN LIST Y ")
+
+    x_min,x_max = np.min(x), np.max(x)
+    y_min,y_max = np.min(y), np.max(y)
+    plpy.notice(x_min)
+    plpy.notice(x_max)
+    plpy.notice(y_min)
+    plpy.notice(y_max)
+    plpy.notice(mesh_size)
+
+    x_grid = np.linspace(x_min,x_max, mesh_size)
+    y_grid = np.linspace(y_min,y_max, mesh_size)
+    range  = [[x_min,x_max],[y_min,y_max]]
+    a, xedges, yedges= np.histogram2d(x,y,bins=(mesh_size,mesh_size), range=range)
+    a = np.swapaxes(a,0,1)
+    plpy.notice("here about to create the contours")
+
+    CS = plt.contour(xedges[1:],yedges[1:] ,a,4,linewidths=0.5, colors='b')
+    plpy.notice(levels)
+    return[(contour_to_polygon(CS.Cntr.trace((level))[0]), float(level)) for level in levels]

python/crankshaft/setup.py

@@ -10,7 +10,7 @@ from setuptools import setup, find_packages
 setup(
     name='crankshaft',
 
-    version='0.0.01',
+    version='0.0.1',
 
     description='CartoDB Spatial Analysis Python Library',
 
@@ -40,9 +40,9 @@ setup(
 
     # The choice of component versions is dictated by what's
     # provisioned in the production servers.
-    install_requires=['pysal==1.11.0','numpy==1.6.1','scipy==0.17.0'],
+    install_requires=['pysal==1.11.0','numpy==1.10.1','scipy==0.17.0', 'matplotlib==1.4.3'],
 
-    requires=['pysal', 'numpy'],
+    requires=['pysal', 'numpy', 'matplotlib'],
 
     test_suite='test'
 )