democall.sql

-- SELECT Faraday_JoyDivision('households'::text, 'year_built'::text) AS the_geom

explain analyze WITH

-- Snap the disparate points to a regular grid where theres data for the given filter
grid AS (
  SELECT
    -- note transform here
    ST_SnapToGrid(ST_Transform(wkb_geometry, 4326), 0.002) the_geom,
    year_built
  FROM
    households
  WHERE
    year_built IS NOT NULL AND
    LEFT(id, 2) = '77'
),
-- Assign Group IDs and average filter value based on grid position
-- Another consequence is that the whole dataset has to be loaded in memory
-- before starting to compute the average ("stats" CTE).
-- Computing the aggregate directly in a single CTE would help,
-- as well as directly filtering for intersection with bounding box of an area of interest.

stats AS (
  SELECT
    the_geom,
    avg(year_built) avg_field,
    ST_Y(the_geom) pulserow,
    ST_X(the_geom) pulseorder
  FROM
    grid
  GROUP BY
    the_geom
),
-- Get the range of filter values in the dataset
varbounds AS (
  SELECT 
    min(avg_field) AS minfield, 
    max(avg_field) AS maxfield
  FROM stats 
),
-- Offset the points to the North based on the mean filter value and a multiplier
offsetsup AS (
  SELECT
    stats.pulserow,
    stats.pulseorder,
    ST_Translate(stats.the_geom, 0, (((stats.avg_field - varbounds.minfield)/(varbounds.maxfield - varbounds.minfield)) * 0.01)) the_geom
  FROM
    stats, varbounds
),
-- Bring the point rows back down to a level where the lowest value is not offset North at all
offsets AS (
  SELECT
    stats.pulserow,
    stats.pulseorder,
    ST_Translate(offsetsup.the_geom, 0, (stats.pulserow - ST_Y(offsetsup.the_geom))) the_geom
  FROM
    stats, offsetsup
),

-- Build horizontal lines on the variably-offset point rows
lines AS (
  SELECT
    pulserow,
    ST_MakeLine(the_geom ORDER BY pulseorder ASC) the_geom
  FROM
    offsets
  GROUP BY
    pulserow
),
-- Define bounding polygons around the active data
hulls AS (
  SELECT 
    round(pulserow::numeric,1)::text || round(pulseorder::numeric,1)::text groupHULL, 
    ST_Buffer(ST_ConcaveHull(ST_Collect(the_geom), 0.99),0.01) the_geom  -- ST_Buffer is expensive
  FROM 
    offsets
  GROUP BY 
    groupHull
)
-- Clip the horizontal lines to the extent of the bounding polygons and return that geometry
-- one important thing to keep in mind is that CTE [ with cte_name as (...) ] will be fully computed by postgresql,
-- without applying filters that has been specified outside of them. 
-- this means that the ST_Intersects filter will not be effective
-- and you'll always be selecting all geometries from the original table in the first CTE (grid).

SELECT 
  st_intersection(hulls.the_geom, lines.the_geom) the_geom
FROM 
  hulls, lines
WHERE 
  st_intersects(lines.the_geom, hulls.the_geom)

joydivision.sql

CREATE OR REPLACE FUNCTION Faraday_JoyDivision(tablename TEXT,field TEXT,grid_tolerance FLOAT DEFAULT 0.002,offset_mult FLOAT DEFAULT 0.01) RETURNS TABLE(the_geom GEOMETRY) AS $$
DECLARE
  query TEXT;
BEGIN
  query := '
  WITH
 
  -- Snap the disparate points to a regular grid where theres data for the given filter
  grid AS (
    SELECT
      -- note transform here
      ST_SnapToGrid(ST_Transform(wkb_geometry, 4326), $1) the_geom,
      '|| field ||'
    FROM
      '|| tablename ||'
    WHERE
      '|| field ||' IS NOT NULL
  ),
 -- Assign Group IDs and average filter value based on grid position
  stats AS (
    SELECT
      the_geom,
      avg('|| field ||') avg_field,
      ST_Y(the_geom) pulserow,
      ST_X(the_geom) pulseorder
    FROM
      grid
    GROUP BY
      the_geom
  ),
-- Get the range of filter values in the dataset
  varbounds AS (
    SELECT 
      min(avg_field) AS minfield, 
      max(avg_field) AS maxfield
    FROM stats 
  ),
 -- Offset the points to the North based on the mean filter value and a multiplier
  offsetsup AS (
    SELECT
      stats.pulserow,
      stats.pulseorder,
      ST_Translate(stats.the_geom, 0, (((stats.avg_field - varbounds.minfield)/(varbounds.maxfield - varbounds.minfield)) * $2)) the_geom
    FROM
      stats, varbounds
  ),
-- Bring the point rows back down to a level where the lowest value is not offset North at all
  offsets AS (
    SELECT
      stats.pulserow,
      stats.pulseorder,
      ST_Translate(offsetsup.the_geom, 0, (stats.pulserow - ST_Y(offsetsup.the_geom))) the_geom
    FROM
      stats, offsetsup
  ),

-- Build horizontal lines on the variably-offset point rows
  lines AS (
    SELECT
      pulserow,
      ST_MakeLine(the_geom ORDER BY pulseorder ASC) the_geom
    FROM
      offsets
    GROUP BY
      pulserow
  ),
 -- Define bounding polygons around the active data
  hulls AS (
    SELECT 
      round(pulserow::numeric,1)::text || round(pulseorder::numeric,1)::text groupHULL, 
      ST_Buffer(ST_ConcaveHull(ST_Collect(the_geom), 0.99),0.01) the_geom
    FROM 
      offsets
    GROUP BY 
      groupHull
  )
 -- Clip the horizontal lines to the extent of the bounding polygons and return that geometry
  SELECT 
    st_intersection(a.the_geom, b.the_geom) the_geom
  FROM 
    hulls a, lines b
  WHERE 
    st_intersects(b.the_geom, a.the_geom)';
  RETURN QUERY EXECUTE query USING grid_tolerance, offset_mult;
END
$$ LANGUAGE 'plpgsql' STABLE;