Last active
August 5, 2019 17:32
-
-
Save darrell/4f6918e5e12e1e68753d2e1339435f75 to your computer and use it in GitHub Desktop.
How to calculate the median population center of Oregon using PostGIS
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
-- 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; |
And before you ask, yes this is very close to the outlet mall in Wilsonville.
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment
The finished product.