-- this assumes you have loaded the census blockgroup table
-- into postgres already. You should also reproject into a state plane system
-- in this case, I used EPSG:2992 which is the preferred state wide projection
-- for Oregon.


-- do some tidying on the block group files, to make the boundaries the same
-- as the state boundary data I'm using.
UPDATE acs_bg_2017 a
  UPDATE acs_bg_2017 a
  SET geom=st_multi(st_buffer(st_intersection(a.geom,s.geom),0.000))
   FROM or_state_boundary s
   WHERE s.feature=3;

-- next we're going to generate a point for each person in the block group
-- and place it randomly within that block group.
-- this creates one MultiPoint feature per block group
DROP TABLE IF EXISTS pop_pts, foo;
CREATE TABLE  pop_pts as
  SELECT
    geoid_data,
    st_area(geom) as area
    ,pop17
    ,case pop17 when 0 then 0 else pop17::numeric/(st_area(geom)/(5280*5280)) end as density
    ,st_generatepoints(geom,pop17) as geom
  FROM acs_bg_2017;

-- Next we want to dump all of those points from all of those blockgroups
-- into a single table with one row per point.
CREATE TEMP TABLE foo as
  SELECT st_x(geom) as lon -- not technically lon or lat, sue me.
  ,st_y(geom) lat
  FROM (SELECT (st_dumppoints(geom)).geom FROM pop_pts) a ;


DROP TABLE IF EXISTS median_pt;
CREATE TABLE median_pt(geom geometry(Point,2992)); -- to keep QGIS happy by defining the projection

-- We calculate the median lat/northing and lon/easting for all the points
-- and create a single point at those coordinates.
INSERT INTO median_pt
  SELECT
  st_setsrid(
    st_makepoint(
      (select percentile_disc(0.5) within group (order by lon) from foo),
      (select percentile_disc(0.5) within group (order by lat) from foo)
    )
  ,2992);

-- Now we want to generate lines N-S and E-W extending from the bounding box
-- of the state through the median point we found above.

DROP TABLE IF EXISTS med_lines;
  CREATE TABLE med_lines as
WITH m as (
    WITH p as (select (st_dumppoints(st_extent(geom))).geom from acs_bg_2017)
    SELECT
    max(st_x(p.geom)) as max_x,
    min(st_x(p.geom)) as min_x,
    max(st_y(p.geom)) as max_y,
    min(st_y(p.geom)) as min_y,
    avg(st_x(m.geom)) as med_x,
    avg(st_y(m.geom)) as med_y
    from p,median_pt as m
  )
SELECT st_setsrid(st_makeline(ARRAY[st_makepoint(max_x,med_y),st_makepoint(min_x,med_y)]),2992) as geom
  FROM m
UNION
SELECT st_setsrid(st_makeline(ARRAY[st_makepoint(med_x,max_y),st_makepoint(med_x,min_y)]),2992)
    from m;