5EdbsdyK

1-- this assumes you have loaded the census blockgroup table 
2-- into postgres already. You should also reproject into a state plane system
3-- in this case, I used EPSG:2992 which is the preferred state wide projection
4-- for Oregon.
5
6
7-- do some tidying on the block group files, to make the boundaries the same
8-- as the state boundary data I'm using.
9UPDATE acs_bg_2017 a 
10  UPDATE acs_bg_2017 a 
11  SET geom=st_multi(st_buffer(st_intersection(a.geom,s.geom),0.000)) 
12   FROM or_state_boundary s 
13   WHERE s.feature=3;
14
15-- next we're going to generate a point for each person in the block group
16-- and place it randomly within that block group.
17-- this creates one MultiPoint feature per block group
18DROP TABLE IF EXISTS pop_pts, foo;
19CREATE TABLE  pop_pts as  
20  SELECT  
21    geoid_data,
22    st_area(geom) as area
23    ,pop17
24    ,case pop17 when 0 then 0 else pop17::numeric/(st_area(geom)/(5280*5280)) end as density
25    ,st_generatepoints(geom,pop17) as geom 
26  FROM acs_bg_2017;
27
28-- Next we want to dump all of those points from all of those blockgroups
29-- into a single table with one row per point.
30CREATE TEMP TABLE foo as  
31  SELECT st_x(geom) as lon -- not technically lon or lat, sue me.
32  ,st_y(geom) lat 
33  FROM (SELECT (st_dumppoints(geom)).geom FROM pop_pts) a ;
34
35
36DROP TABLE IF EXISTS median_pt;
37CREATE TABLE median_pt(geom geometry(Point,2992)); -- to keep QGIS happy by defining the projection
38
39-- We calculate the median lat/northing and lon/easting for all the points
40-- and create a single point at those coordinates.
41INSERT INTO median_pt 
42  SELECT 
43  st_setsrid(
44    st_makepoint(
45      (select percentile_disc(0.5) within group (order by lon) from foo),
46      (select percentile_disc(0.5) within group (order by lat) from foo)
47    )
48  ,2992);
49  
50-- Now we want to generate lines N-S and E-W extending from the bounding box 
51-- of the state through the median point we found above.
52
53DROP TABLE IF EXISTS med_lines;
54  CREATE TABLE med_lines as
55WITH m as (
56    WITH p as (select (st_dumppoints(st_extent(geom))).geom from acs_bg_2017)
57    SELECT 
58    max(st_x(p.geom)) as max_x,
59    min(st_x(p.geom)) as min_x,
60    max(st_y(p.geom)) as max_y,
61    min(st_y(p.geom)) as min_y,
62    avg(st_x(m.geom)) as med_x,
63    avg(st_y(m.geom)) as med_y
64    from p,median_pt as m
65  )
66SELECT st_setsrid(st_makeline(ARRAY[st_makepoint(max_x,med_y),st_makepoint(min_x,med_y)]),2992) as geom
67  FROM m
68UNION
69SELECT st_setsrid(st_makeline(ARRAY[st_makepoint(med_x,max_y),st_makepoint(med_x,min_y)]),2992)
70    from m;