diff --git a/src/pg/sql/50_contours.sql b/src/pg/sql/50_contours.sql index 545f518..5501160 100644 --- a/src/pg/sql/50_contours.sql +++ b/src/pg/sql/50_contours.sql @@ -21,10 +21,12 @@ CREATE OR REPLACE FUNCTION bandwidth NUMERIC DEFAULT 0.0001, levels NUMERIC[] DEFAULT null ) -RETURNS table (level Numeric, geom text ) +RETURNS table (level Numeric, geom geometry ) AS $$ BEGIN + RETURN QUERY - select l as level, ST_GeomFromText(geom_text, 4326) as geom from _CDB_Contours(subquery,grid_size,bandwidth,levels); + 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; + diff --git a/src/py/crankshaft/crankshaft/contours/contours.py b/src/py/crankshaft/crankshaft/contours/contours.py index feb249d..14c55f8 100644 --- a/src/py/crankshaft/crankshaft/contours/contours.py +++ b/src/py/crankshaft/crankshaft/contours/contours.py @@ -5,15 +5,20 @@ from sklearn.neighbors import KernelDensity from skimage.measure import find_contours import plpy - def cdb_generate_contours(query, grid_size, bandwidth, levels): - data = plpy.execute( 'select ST_X(the_geom) as x , ST_Y(the_geom) as y from ({query})') - xs, ys = [d['x'], d['y'] from data] - return generate_contours(xs,xy,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): - return [coord[0]*(x_range[1]-x_range[0])/grid_size+x_range[0], - coord[1]*(y_range[1]-y_range[0])/grid_size+y_range[0]] + 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]) @@ -23,15 +28,17 @@ 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.0001, levels=None): - max_y, min_y = ys.max(), ys.min() - max_x, min_x = xs.max(), xs.min() +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. - kde = KernelDensity(bandwidth=0.0001, metric='haversine', + 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)) @@ -39,7 +46,7 @@ def generate_contours(xs,ys, grid_res=100, bandwidth=0.0001, levels=None): 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 = [] @@ -47,5 +54,5 @@ def generate_contours(xs,ys, grid_res=100, bandwidth=0.0001, levels=None): 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