SHARE
TWEET

Untitled

a guest Aug 17th, 2019 67 Never
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
  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.
  9. UPDATE 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
  18. DROP TABLE IF EXISTS pop_pts, foo;
  19. CREATE 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.
  30. CREATE 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.  
  36. DROP TABLE IF EXISTS median_pt;
  37. CREATE 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.
  41. INSERT 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.  
  53. DROP TABLE IF EXISTS med_lines;
  54.   CREATE TABLE med_lines as
  55. WITH 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.   )
  66. SELECT st_setsrid(st_makeline(ARRAY[st_makepoint(max_x,med_y),st_makepoint(min_x,med_y)]),2992) as geom
  67.   FROM m
  68. UNION
  69. SELECT st_setsrid(st_makeline(ARRAY[st_makepoint(med_x,max_y),st_makepoint(med_x,min_y)]),2992)
  70.     from m;
RAW Paste Data
We use cookies for various purposes including analytics. By continuing to use Pastebin, you agree to our use of cookies as described in the Cookies Policy. OK, I Understand
 
Top