ST_Intersection not returning all intersecting geometries

I have a query that will calculate the proportions of given columns based on the area of intersection. The returns mostly the polygons entirely inside the other polygon, and in some cases returns the ones that intersect. Why is it not returning polygons entirely inside AND those bisected by the boundaries of another?

SELECT
ROW_NUMBER() OVER () AS my_new_id,
st_intersection(a.geom, b.geom)
into newlayer
from mass_pop_censustract as a, mass_regions as b


EDITED: The query should return the bright green plus the areas in the blue hatch (right picture). Essentially, I want to calculate the proportion of the population that lives in the portion of the census tracts (yellow) that falls in my area of interest (light green, left picture.)

• I think you want to use ST_Intersects...? – DPSSpatial Jul 11 '18 at 20:18
• I am using ST_Intersection because I want it to return a geometry. So that line is giving me my <geom> column. ST_Intersects would return true or false right? For instance, if I wanted to check if something intersects. – atlasofcoffee Jul 11 '18 at 20:25
• I think ST_Intersects would return ones that are bisected... ST_Intersection just returns the portion that is 'clipped' by the overlapping features, by ST_Intersects would return features where a portion intersects. In terms of the true/false, I think that's just the way the MS SQL Server version works... – DPSSpatial Jul 11 '18 at 20:31
• So if I want a geometry column that shows both the polygons entirely inside and those bisected, how would I get that with ST_Intersects? – atlasofcoffee Jul 11 '18 at 20:36
• I think in this case you might need a combination of both - ST_Intersects to figure out which polygons actually intersect, and ST_Intersection to calculate the area of those bisected... does that make sense? – DPSSpatial Jul 11 '18 at 22:28

Following my comment above, here's a suggestion to get the areal proportion of those intersections directly, to then be used in further steps, instead of your initial workflow. Might help, if that is actually what you want.

Since you didn't mention the data's CRS I'll assume WGS84 (EPSG:4326) and use a cast to geography type to derive measurements based on spheroidal algebra and in metric units.
I'll also assume the census tract data has an id, and a pop column with each tracts population.

Running

WITH
its AS (
SELECT ct.id,
ct.pop,
ST_Area(ct.geom::geography) AS area,
ST_Area(ST_Intersection(ct.geom, rg.geom)::geography) AS prop_area
FROM <census> AS ct
JOIN <region> AS rg
ON ST_Intersects(ct.geom, rg.geom)
)

SELECT id,
prop_area / 1000 AS prop_area_sqkm,
prop_area / area AS prop_area_pct,
pop * (prop_area / area) AS prop_pop
FROM its;


will calculate prop_pop as population proportion for each census tract id based on percentage of the area intersecting the ROI.

If you need to join this with other data holding population numbers, the query returns the respective intersection area prop_area_sqkm (if you are working with e.g. population density), and that area's proportion prop_area_pct in percentage (for e.g. external total population data).

Does that help?

• let me add: obviously, calculating the population proportion this way can not be more than a purely statistical assumption, and no reliable numbers can come from that (i.e. if all major urban areas of a tract happen to not fall inside your ROI, the above proportion is plain wrong). I felt the need to add that, but, of course, I don't know your usecase... – ThingumaBob Jul 13 '18 at 17:33
• This is very useful, thank you! I just ran it, it's taking a long time because the population data are block groups in California. Will report the outcome, thanks again. – atlasofcoffee Jul 13 '18 at 19:12
• @atlasofcoffee okay, let's see. always make sure you have spatial indexes on both tables! run VACUUM ANALYZE <table_name> to update table stats. it'l take some time, but this shouldn't take ages...if in doubt, run the query with EXPLAIN ANALYZE <QUERY> and post the query plan. there might be ways to speed things up. – ThingumaBob Jul 13 '18 at 19:19
• So the query ran, calculations work but the issue of the geometry persists. It doesn't seem to be capturing the portions of some of the census blocks that overlap with the areas of interest. Is this due to geometry problems? – atlasofcoffee Jul 13 '18 at 21:04
• @atlasofcoffee bit late, sry. well, yes, that's very likely the reason, as most ST functions simply return NULL for invald geometries. you could try SELECT *, ST_IsValidReason(geom) FROM <table> WHERE NOT ST_IsValid(geom) and see if you get results. if indeed there are invalid geometries, you'd get some hints on how to fix them. – ThingumaBob Jul 18 '18 at 13:20

the two queries below should get you what you need. you need to use the st_intersects to spatially filter the records that overlap. (I believe st_intersects first checks for bounding box intersections)

drop table if exists new_layer;
create table new_layer as
SELECT
ROW_NUMBER() OVER () AS my_new_id,
st_intersection(a.geom, b.geom)
from mass_pop_censustract as a join mass_regions as b
on st_intersects(a.geom, b.geom)


or this variation is sometimes useful

drop table if exists new_layer;
create table new_layer as
SELECT
ROW_NUMBER() OVER () AS my_new_id,
st_intersection(a.geom, b.geom)
from mass_pop_censustract as a join mass_regions as b
on st_dwithin(a.geom, b.geom,0.001)

• Thank you for your answer. I will add a screenshot of the result, it still yields only the polygons entirely inside. – atlasofcoffee Jul 11 '18 at 23:57
• I don't think you are understanding what ST_Intersection does: it will ONLY return the shared parts of the geometries from the two tables, it is suppossed to clip it/cut it. if a geometry portion of table a is not directly hit by a geometry of table b it will not show up. and vice versa – ziggy Jul 12 '18 at 1:27