cartodb4/observatory-extension
1
0
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity

fix non-noded intersection between shoreline clipped and non-shoreline clipped geometries by using a safe_intersection function

Browse Source
This commit is contained in:
John Krauss 2017-03-07 20:38:12 +00:00
parent fd3918b29c
commit deede798e9
7 changed files with 207 additions and 20 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

6
scripts/generate_fixtures.py
View File

@ -214,6 +214,7 @@ FIXTURES = [
('us.census.tiger.fullname', 'us.census.tiger.pointlm_geom', '2016'),
('us.census.tiger.fullname', 'us.census.tiger.prisecroads_geom', '2016'),
('us.census.tiger.name', 'us.census.tiger.county', '2015'),
('us.census.tiger.name', 'us.census.tiger.county_clipped', '2015'),
('us.census.tiger.name', 'us.census.tiger.block_group', '2015'),
]
@ -358,7 +359,10 @@ def main():
dump('*', tablename, 'WHERE geom && ST_MakeEnvelope(-74,40.69,-73.9,40.72, 4326)')
continue
elif 'whosonfirst' in table_id:
where = '(\'85632785\',\'85633051\',\'85633111\',\'85633147\',\'85633253\',\'85633267\')'
where = "('85632785','85633051','85633111','85633147','85633253','85633267')"
compare = 'IN'
elif 'county' in table_id and 'tiger' in table_id:
where = "('48061', '36047')"
compare = 'IN'
else:
where = '\'36047%\''

19
src/pg/sql/40_observatory_utility.sql
View File

@ -231,3 +231,22 @@ CREATE AGGREGATE cdb_observatory.FIRST (
basetype = anyelement,
stype = anyelement
);
-- Attempt to perform intersection, if there's an exception then buffer
-- https://gis.stackexchange.com/questions/50399/how-best-to-fix-a-non-noded-intersection-problem-in-postgis
CREATE OR REPLACE FUNCTION cdb_observatory.safe_intersection(geom_a geometry, geom_b geometry)
RETURNS geometry AS
$$
BEGIN
RETURN ST_Intersection(geom_a, geom_b);
EXCEPTION
WHEN OTHERS THEN
BEGIN
RETURN ST_Intersection(ST_Buffer(geom_a, 0.0000001), ST_Buffer(geom_b, 0.0000001));
EXCEPTION
WHEN OTHERS THEN
RETURN ST_GeomFromText('POLYGON EMPTY');
END;
END
$$
LANGUAGE 'plpgsql' STABLE STRICT;

37
src/pg/sql/41_observatory_augmentation.sql
View File

@ -582,8 +582,8 @@ BEGIN
(normalization IS NULL AND LOWER(denom_reltype) LIKE 'denominator')
THEN ' CASE ' ||
-- denominated point-in-poly or user polygon is same as OBS polygon
' WHEN ST_GeometryType(cdb_observatory.FIRST(_geoms.geom)) = ''ST_Point'' ' ||
' OR cdb_observatory.FIRST(_geoms.geom = ' || geom_tablename || '.' || geom_colname || ')' ||
' WHEN EVERY(ST_GeometryType(_geoms.geom) = ''ST_Point'') ' ||
' OR EVERY(_geoms.geom::TEXT = ' || geom_tablename || '.' || geom_colname || '::TEXT)' ||
' THEN cdb_observatory.FIRST(' || numer_tablename || '.' || numer_colname ||
' / NullIf(' || denom_tablename || '.' || denom_colname || ', 0))' ||
-- denominated polygon interpolation
@ -594,7 +594,7 @@ BEGIN
' THEN ST_Area(_geoms.geom) / Nullif(ST_Area(' || geom_tablename || '.' || geom_colname || '), 0) ' ||
' WHEN ST_Within(' || geom_tablename || '.' || geom_colname || ', _geoms.geom) ' ||
' THEN 1 ' ||
' ELSE (ST_Area(ST_MakeValid(ST_Intersection(_geoms.geom, ' || geom_tablename || '.' || geom_colname || '))) ' ||
' ELSE (ST_Area(ST_MakeValid(cdb_observatory.safe_intersection(_geoms.geom, ' || geom_tablename || '.' || geom_colname || '))) ' ||
' / Nullif(ST_Area(' || geom_tablename || '.' || geom_colname || '), 0))' ||
' END) / '
' NULLIF(SUM(' || denom_tablename || '.' || denom_colname || ' ' ||
@ -602,7 +602,7 @@ BEGIN
' THEN ST_Area(_geoms.geom) / Nullif(ST_Area(' || geom_tablename || '.' || geom_colname || '), 0) ' ||
' WHEN ST_Within(' || geom_tablename || '.' || geom_colname || ', _geoms.geom) ' ||
' THEN 1 ' ||
' ELSE (ST_Area(ST_MakeValid(ST_Intersection(_geoms.geom, ' || geom_tablename || '.' || geom_colname || '))) ' ||
' ELSE (ST_Area(ST_MakeValid(cdb_observatory.safe_intersection(_geoms.geom, ' || geom_tablename || '.' || geom_colname || '))) ' ||
' / Nullif(ST_Area(' || geom_tablename || '.' || geom_colname || '), 0))' ||
' END), 0) ' ||
' / (COUNT(*) / COUNT(distinct ' || geom_tablename || '.' || geom_geomref_colname || ')) ' ||
@ -612,8 +612,8 @@ BEGIN
(normalization IS NULL AND numer_aggregate ILIKE 'sum')
THEN ' CASE ' ||
-- areaNormalized point-in-poly or user polygon is the same as OBS polygon
' WHEN ST_GeometryType(cdb_observatory.FIRST(_geoms.geom)) = ''ST_Point'' ' ||
' OR cdb_observatory.FIRST(_geoms.geom = ' || geom_tablename || '.' || geom_colname || ')' ||
' WHEN EVERY(ST_GeometryType(_geoms.geom) = ''ST_Point'') ' ||
' OR EVERY(_geoms.geom::TEXT = ' || geom_tablename || '.' || geom_colname || '::TEXT)' ||
' THEN cdb_observatory.FIRST(' || numer_tablename || '.' || numer_colname ||
' / (Nullif(ST_Area(' || geom_tablename || '.' || geom_colname || '::Geography), 0)/1000000)) ' ||
-- areaNormalized polygon interpolation
@ -625,7 +625,7 @@ BEGIN
' WHEN ST_Within(' || geom_tablename || '.' || geom_colname || ', _geoms.geom) THEN ' ||
' ST_Area(' || geom_tablename || '.' || geom_colname || ') ' ||
' / Nullif(ST_Area(_geoms.geom), 0) ' ||
' ELSE (ST_Area(ST_MakeValid(ST_Intersection(_geoms.geom, ' || geom_tablename || '.' || geom_colname || '))) ' ||
' ELSE (ST_Area(ST_MakeValid(cdb_observatory.safe_intersection(_geoms.geom, ' || geom_tablename || '.' || geom_colname || '))) ' ||
' / Nullif(ST_Area(_geoms.geom), 0))' ||
' END / (Nullif(ST_Area(' || geom_tablename || '.' || geom_colname || '::Geography), 0) / 1000000)) ' ||
' / (COUNT(*) / COUNT(distinct ' || geom_tablename || '.' || geom_geomref_colname || ')) ' ||
@ -636,8 +636,8 @@ BEGIN
(normalization IS NULL OR LOWER(normalization) LIKE 'pre%')
THEN ' CASE ' ||
-- predenominated point-in-poly or user polygon is the same as OBS- polygon
' WHEN ST_GeometryType(cdb_observatory.FIRST(_geoms.geom)) = ''ST_Point'' ' ||
' OR cdb_observatory.FIRST(_geoms.geom = ' || geom_tablename || '.' || geom_colname || ')' ||
' WHEN EVERY(ST_GeometryType(_geoms.geom) = ''ST_Point'') ' ||
' OR EVERY(_geoms.geom::TEXT = ' || geom_tablename || '.' || geom_colname || '::TEXT)' ||
' THEN cdb_observatory.FIRST(' || numer_tablename || '.' || numer_colname || ') ' ||
' ELSE ' ||
-- predenominated polygon interpolation weighted by universe
@ -650,7 +650,7 @@ BEGIN
' THEN ST_Area(_geoms.geom) / Nullif(ST_Area(' || geom_tablename || '.' || geom_colname || '), 0) ' ||
' WHEN ST_Within(' || geom_tablename || '.' || geom_colname || ', _geoms.geom) ' ||
' THEN 1 ' ||
' ELSE (ST_Area(ST_MakeValid(ST_Intersection(_geoms.geom, ' || geom_tablename || '.' || geom_colname || '))) ' ||
' ELSE (ST_Area(ST_MakeValid(cdb_observatory.safe_intersection(_geoms.geom, ' || geom_tablename || '.' || geom_colname || '))) ' ||
' / Nullif(ST_Area(' || geom_tablename || '.' || geom_colname || '), 0)) ' ||
' END) ' ||
' / Nullif(SUM(' || denom_tablename || '.' || denom_colname ||
@ -658,7 +658,7 @@ BEGIN
' THEN ST_Area(_geoms.geom) / Nullif(ST_Area(' || geom_tablename || '.' || geom_colname || '), 0) ' ||
' WHEN ST_Within(' || geom_tablename || '.' || geom_colname || ', _geoms.geom) ' ||
' THEN 1 ' ||
' ELSE (ST_Area(ST_MakeValid(ST_Intersection(_geoms.geom, ' || geom_tablename || '.' || geom_colname || '))) ' ||
' ELSE (ST_Area(ST_MakeValid(cdb_observatory.safe_intersection(_geoms.geom, ' || geom_tablename || '.' || geom_colname || '))) ' ||
' / Nullif(ST_Area(' || geom_tablename || '.' || geom_colname || '), 0))' ||
' END), 0) ' ||
' / (COUNT(*) / COUNT(distinct ' || geom_tablename || '.' || geom_geomref_colname || ')) ' ||
@ -668,8 +668,8 @@ BEGIN
(normalization IS NULL OR LOWER(normalization) LIKE 'pre%')
THEN ' CASE ' ||
-- predenominated point-in-poly or user polygon is the same as OBS- polygon
' WHEN ST_GeometryType(cdb_observatory.FIRST(_geoms.geom)) = ''ST_Point'' ' ||
' OR cdb_observatory.FIRST(_geoms.geom = ' || geom_tablename || '.' || geom_colname || ')' ||
' WHEN EVERY(ST_GeometryType(_geoms.geom) = ''ST_Point'') ' ||
' OR EVERY(_geoms.geom::TEXT = ' || geom_tablename || '.' || geom_colname || '::TEXT)' ||
' THEN cdb_observatory.FIRST(' || numer_tablename || '.' || numer_colname || ') ' ||
' ELSE ' ||
-- predenominated polygon interpolation
@ -679,15 +679,16 @@ BEGIN
' THEN ST_Area(_geoms.geom) / Nullif(ST_Area(' || geom_tablename || '.' || geom_colname || '), 0) ' ||
' WHEN ST_Within(' || geom_tablename || '.' || geom_colname || ', _geoms.geom) ' ||
' THEN 1 ' ||
' ELSE (ST_Area(ST_MakeValid(ST_Intersection(_geoms.geom, ' || geom_tablename || '.' || geom_colname || '))) ' ||
' ELSE (ST_Area(ST_MakeValid(cdb_observatory.safe_intersection(_geoms.geom, ' ||
geom_tablename || '.' || geom_colname || '))) ' ||
' / Nullif(ST_Area(' || geom_tablename || '.' || geom_colname || '), 0))' ||
' END) ' ||
' / (COUNT(*) / COUNT(distinct ' || geom_tablename || '.' || geom_geomref_colname || ')) ' ||
'END '
-- Everything else. Point only!
ELSE ' CASE ' ||
' WHEN ST_GeometryType(cdb_observatory.FIRST(_geoms.geom)) = ''ST_Point'' ' ||
' OR cdb_observatory.FIRST(_geoms.geom = ' || geom_tablename || '.' || geom_colname || ')' ||
' WHEN EVERY(ST_GeometryType(_geoms.geom) = ''ST_Point'') ' ||
' OR EVERY(_geoms.geom::TEXT = ' || geom_tablename || '.' || geom_colname || '::TEXT)' ||
' THEN cdb_observatory.FIRST(' || numer_tablename || '.' || numer_colname || ') ' ||
' ELSE cdb_observatory._OBS_RaiseNotice(''Cannot perform calculation over polygon for ' ||
numer_id || '/' || coalesce(denom_id, '') || '/' || geom_id || '/' || numer_timespan || ''')::Numeric ' ||
@ -705,8 +706,8 @@ BEGIN
'.' || geom_colname || ')::TEXT'
-- code below will return the intersection of the user's geom and the
-- OBS geom
--'"value": "'' || ' || 'cdb_observatory.FIRST(ST_MakeValid(ST_Intersection(_geoms.geom, ' || geom_tablename ||
-- '.' || geom_colname || ')))::TEXT || ''"'''
--'''value'', ' || 'ST_Union(ST_MakeValid(cdb_observatory.safe_intersection(_geoms.geom, ' || geom_tablename ||
-- '.' || geom_colname || ')))::TEXT'
ELSE ''
END || ')', ', ')
AS colspecs,

4
src/pg/test/expected/40_observatory_utility_test.out
View File

@ -21,3 +21,7 @@ t
obs_dumpversion_notnull
t
(1 row)
ERROR: Error performing intersection: TopologyException: found non-noded intersection between LINESTRING (-97.1968 25.9574, -97.1971 25.9576) and LINESTRING (-97.197 25.9575, -97.1972 25.9576) at -97.19699802694231 25.957551976080605
complex_safe_intersection_works
t
(1 row)

1
src/pg/test/fixtures/drop_fixtures.sql vendored
View File

@ -26,6 +26,7 @@ DROP TABLE IF EXISTS observatory.obs_6c1309a64d8f3e6986061f4d1ca7b57743e75e74;
DROP TABLE IF EXISTS observatory.obs_0310c639744a2014bb1af82709228f05b59e7d3d;
DROP TABLE IF EXISTS observatory.obs_87a814e485deabe3b12545a537f693d16ca702c2;
DROP TABLE IF EXISTS observatory.obs_e32f8e59c7c8861ee5ee4029b3ace2af9a5c9caf;
DROP TABLE IF EXISTS observatory.obs_23cb5063486bd7cf36f17e89e5e65cd31b331f6e;
DROP TABLE IF EXISTS observatory.obs_1ea93bbc109c87c676b3270789dacf7a1430db6c;
DROP TABLE IF EXISTS observatory.obs_b393b5b88c6adda634b2071a8005b03c551b609a;
DROP TABLE IF EXISTS observatory.obs_8e30e6b3792430b410ba5b9e49cdc6a0d404d48f;

146
src/pg/test/fixtures/load_fixtures.sql vendored
View File

File diff suppressed because one or more lines are too long

12
src/pg/test/sql/40_observatory_utility_test.sql
View File

@ -47,3 +47,15 @@ SELECT cdb_observatory._OBS_StandardizeMeasureName('test 343 %% 2 qqq }}{{}}') =
SELECT cdb_observatory.OBS_DumpVersion()
IS NOT NULL AS OBS_DumpVersion_notnull;
-- Should fail to perform intersection
SELECT ST_IsValid(ST_Intersection(
cdb_observatory.OBS_GetBoundaryByID('48061', 'us.census.tiger.county'),
cdb_observatory.OBS_GetBoundaryByID('48061', 'us.census.tiger.county_clipped')
)) AS complex_intersection_fails;
-- Should succeed in intersecting
SELECT ST_IsValid(cdb_observatory.safe_intersection(
cdb_observatory.OBS_GetBoundaryByID('48061', 'us.census.tiger.county'),
cdb_observatory.OBS_GetBoundaryByID('48061', 'us.census.tiger.county_clipped')
)) AS complex_safe_intersection_works;