Finding differences between polygons in PostGIS

Question

I have another redistricting question. The datasets and queries available on my GitHub, but I will include the queries here, too.

I am working on two complete datasets:

1) A PLAN - composed of individual Proposed Districts 2) Official Census Congressional Districts Boundaries - composed of individual Official District boundaries

What there are two things that I am looking for with these:

1) The INTERSECTION of each Proposed District and Official District 2) The DIFFERENCE of each Proposed District and Official District that it intersects. 2a) The part of the Proposed District that does not overlap with the Official District. 2B) The part of the Official District that does no overlap with the plan.

So these are really three queries that I am looking for, organized into two broad categories. 2a and 2b are essentially mirror images of each other, so if I can see one, I can get the other.

I asked a question about intersections about a year ago and am using that answer to generate (1), which so far seems correct.

-- 1a Congressional Districts

DROP TABLE IF EXISTS
    plan_cong_all_int_src;

CREATE TABLE plan_cong_all_int_src AS (
    SELECT ROW_NUMBER() OVER ()                                              AS pid,
           plan.gid                                                          AS plan_id,
           census.gid                                                        as cong_id,
           plan.district                                                     AS plan_name,
           census.namelsad                                                   AS cong_name,
           plan.geom                                                         AS plan_geom,
           census.geom                                                       AS cong_geom,
           ST_MULTI(ST_BUFFER(ST_INTERSECTION(plan.geom, census.geom), 0.0)) AS intersection_geom
    FROM ca_five_districts AS plan
             INNER JOIN tl_2019_us_cd_116 AS census
                        ON (ST_INTERSECTS(plan.geom, census.geom))
    WHERE NOT ST_ISEMPTY(ST_BUFFER(ST_INTERSECTION(plan.geom, census.geom), 0.0)));

CREATE INDEX pln_cong_sdx
    ON
        plan_cong_all_int_src
            USING GIST (intersection_geom);

DROP TABLE IF EXISTS plan_cong_all_int_rep;

CREATE TABLE plan_cong_all_int_rep AS
    (SELECT src.pid                                                 AS pid,
            src.plan_name                                           AS plan_id,
            src.cong_name                                           AS cong_id,
            src.plan_name || ' - ' || src.cong_name                 AS composing_districts,
            ST_AREA(src.intersection_geom) / ST_AREA(src.plan_geom) AS overlap_ratio_for_plan,
            ST_AREA(src.intersection_geom) / ST_AREA(src.cong_geom) AS overlap_ratio_for_congress,
            src.intersection_geom                                   AS geom
     FROM plan_cong_all_int_src src);

But when I try something similar to get the part of the Proposed District that does not intersect with the Official District, the remainder, if you will, I get something that looks wrong to me.

-- Get all parts of the districts that aren't contained by the intersecting CD

DROP TABLE IF EXISTS
    plan_cong_noncntned_int_src;

CREATE TABLE plan_cong_noncntned_int_src AS (
    SELECT ROW_NUMBER() OVER ()                                            AS pid,
           plan.gid                                                        AS plan_id,
           census.gid                                                      as cong_id,
           plan.district                                                   AS plan_name,
           census.namelsad                                                 AS cong_name,
           plan.geom                                                       AS plan_geom,
           census.geom                                                     AS cong_geom,
           ST_MULTI(ST_BUFFER(ST_DIFFERENCE(plan.geom, census.geom), 0.0)) AS difference_geom
    FROM ca_five_districts AS plan,
         tl_2019_us_cd_116 AS census
    WHERE ST_INTERSECTS(plan.geom, census.geom));

CREATE INDEX pln_cong_noncntned_sdx
    ON
        plan_cong_noncntned_int_src
            USING GIST (difference_geom);

DROP TABLE IF EXISTS plan_cong_noncntned_int_rep;

CREATE TABLE plan_cong_noncntned_int_rep AS
    (SELECT src.pid                                               AS pid,
            src.plan_name                                         AS plan_id,
            src.cong_name                                         AS cong_id,
            src.plan_name || ' - ' || src.cong_name               AS composing_districts,
            ST_AREA(src.difference_geom) / ST_AREA(src.plan_geom) AS overlap_ratio_for_plan,
            ST_AREA(src.difference_geom) / ST_AREA(src.cong_geom) AS overlap_ratio_for_congress,
            src.difference_geom                                   AS geom
     FROM plan_cong_noncntned_int_src src);-- Get all parts of the districts that aren't contained by the intersecting CD

DROP TABLE IF EXISTS
    plan_cong_noncntned_int_src;

CREATE TABLE plan_cong_noncntned_int_src AS (
    SELECT ROW_NUMBER() OVER ()                                            AS pid,
           plan.gid                                                        AS plan_id,
           census.gid                                                      as cong_id,
           plan.district                                                   AS plan_name,
           census.namelsad                                                 AS cong_name,
           plan.geom                                                       AS plan_geom,
           census.geom                                                     AS cong_geom,
           ST_MULTI(ST_BUFFER(ST_DIFFERENCE(plan.geom, census.geom), 0.0)) AS difference_geom
    FROM ca_five_districts AS plan,
         tl_2019_us_cd_116 AS census
    WHERE ST_INTERSECTS(plan.geom, census.geom));

CREATE INDEX pln_cong_noncntned_sdx
    ON
        plan_cong_noncntned_int_src
            USING GIST (difference_geom);

DROP TABLE IF EXISTS plan_cong_noncntned_int_rep;

CREATE TABLE plan_cong_noncntned_int_rep AS
    (SELECT src.pid                                               AS pid,
            src.plan_name                                         AS plan_id,
            src.cong_name                                         AS cong_id,
            src.plan_name || ' - ' || src.cong_name               AS composing_districts,
            ST_AREA(src.difference_geom) / ST_AREA(src.plan_geom) AS overlap_ratio_for_plan,
            ST_AREA(src.difference_geom) / ST_AREA(src.cong_geom) AS overlap_ratio_for_congress,
            src.difference_geom                                   AS geom
     FROM plan_cong_noncntned_int_src src);

For what I am trying to do, I would like to get it such that INT + DIFF = DISTRICT.

What I am getting is that INT = DIFF.

What is wrong?

erasing in postgis is a bit tricky, I have found that you need to st_union the geometry you want to cut out. check out this question https://gis.stackexchange.com/questions/250674/postgis-st-difference-similar-to-arcgis-erase — ziggy, Jan 16 '20 at 20:42
@ziggy is this giving you correct results? I tried something similar and got an erroneous answer. — Dheeraj Chand, Jan 17 '20 at 02:10

score 1 · Answer 1 · answered Jan 21 '20 at 15:29

Maybe you can take some inspiration from this post from Paul Ramsey: http://blog.cleverelephant.ca/2019/07/postgis-overlays.html

To summarise: he start by creating independant polygons from the intersections so there is no overlapp, and after, for each of these polygons, he get the number of original polygons that overlapp this one.

You could do something similar: start by creating polygons without overlapp (by cutting at the intersection) and then, for each polygon, get which category it belongs to (only original, new, or both) by intersecting with the first 2 table.

Finding differences between polygons in PostGIS

1 Answers1