Skip to content

Instantly share code, notes, and snippets.

@mapbutcher
Created March 5, 2014 00:41
Show Gist options
  • Save mapbutcher/9358937 to your computer and use it in GitHub Desktop.
Save mapbutcher/9358937 to your computer and use it in GitHub Desktop.
--shp2pgsql -s 4326 -d -I 'SSC_2011_AUST.shp' 'SSC' | psql -d gis -h localhost -U simonhope
--shp2pgsql -s 4326 -d -I 'POA_2011_AUST.shp' 'POA' | psql -d gis -h localhost -U simonhope
--update the data types
ALTER TABLE SSC ALTER COLUMN SSC_CODE TYPE NUMERIC USING SSC_CODE::INTEGER;
ALTER TABLE POA ALTER COLUMN POA_CODE TYPE NUMERIC USING POA_CODE::INTEGER;
--explode multipolygons to single parts
CREATE TABLE POA_EXPLODE AS
SELECT
poa_code,
poa_name,
(ST_Dump(geom) ).geom
FROM
POA;
CREATE TABLE SSC_EXPLODE AS
SELECT
ssc_code,
ssc_name,
(ST_Dump(geom) ).geom
FROM
SSC;
--combine geometries into a single table
DROP TABLE POA_SSC_COMBO;
CREATE TABLE POA_SSC_COMBO AS
SELECT ssc_code AS ssc_code, ssc_name as ssc_name, null AS poa_code, null as poa_name, geom as geom FROM SSC_EXPLODE
UNION ALL
SELECT null AS ssc_code, null as ssc_name, poa_code AS poa_code, poa_name as poa_name, geom as geom FROM POA_EXPLODE;
--creates a single table of non-overalapping polygons
--warning takes a long time to execute..
DROP TABLE POA_SSC_OVERLAY;
CREATE TABLE POA_SSC_OVERLAY AS SELECT
geom
FROM
ST_Dump (
(
SELECT
ST_Polygonize (the_geom) AS the_geom
FROM
(
SELECT
ST_Union (the_geom) AS the_geom
FROM
(
SELECT
ST_ExteriorRing (geom) AS the_geom
FROM
POA_SSC_COMBO
) AS lines
) AS noded_lines
)
);
-- points on surface
DROP TABLE POA_SSC_OVERLAY_PTS;
CREATE TABLE POA_SSC_OVERLAY_PTS AS SELECT
ST_PointOnSurface (geom) AS geom
FROM
POA_SSC_OVERLAY;
-- group by geom and aggregate original ids by point overlap
-- Replicates an ArcGIS-style Union
DROP TABLE POA_SSC_UNION;
CREATE TABLE POA_SSC_UNION AS (
SELECT
NEW .geom AS geom,
MAX (orig.SSC_CODE) AS SSC_CODE,
MIN (orig.POA_CODE) AS POA_CODE
FROM
POA_SSC_COMBO AS orig,
POA_SSC_OVERLAY_PTS AS pt,
POA_SSC_OVERLAY AS NEW
WHERE
orig.geom && pt.geom
AND NEW.geom && pt.geom
AND ST_Intersects (orig.geom, pt.geom)
AND ST_Intersects (NEW .geom, pt.geom)
GROUP BY
NEW .geom
);
-- Join with the original tables to pull in attributes
-- This is still single part geometry
DROP TABLE POA_SSC_UNIONJOIN;
CREATE TABLE POA_SSC_UNIONJOIN AS SELECT
G.geom AS geom,
S.SSC_CODE,
S.SSC_NAME,
P.POA_CODE,
P.POA_NAME
FROM
POA_SSC_UNION AS G
LEFT JOIN SSC S ON S.SSC_CODE = G.SSC_CODE
LEFT JOIN POA P ON P.POA_CODE = G.POA_CODE;
DROP TABLE POA_SSC_FINALUNION;
CREATE TABLE POA_SSC_FINALUNION AS SELECT
U.SSC_CODE,
MIN (U.SSC_NAME) AS SSC_NAME,
U.POA_CODE,
MIN (U.POA_NAME) AS POA_NAME,
ST_Multi (ST_Union(U.geom)) AS GEOM
FROM
POA_SSC_UNIONJOIN AS U
GROUP BY
SSC_CODE,
POA_CODE;
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment