Skip to content

Instantly share code, notes, and snippets.

@fjcerignoni

fjcerignoni/opi_run_overlap.sql

Last active February 20, 2024 14:11
Show Gist options

  • Select an option

    Learn more about clone URLs
  • Save fjcerignoni/245afa10be026be2b0a3aff7d05b4ea9 to your computer and use it in GitHub Desktop.

Select an option

Learn more about clone URLs
Save fjcerignoni/245afa10be026be2b0a3aff7d05b4ea9 to your computer and use it in GitHub Desktop.
Download ZIP
Raw
opi_run_overlap.sql
-- public.st_safe_intersection()
CREATE OR REPLACE FUNCTION ST_Safe_Intersection(geom_a geometry, geom_b geometry)
RETURNS geometry
LANGUAGE plpgsql
IMMUTABLE STRICT
AS $function$
BEGIN
RETURN ST_Intersection(geom_a, geom_b);
EXCEPTION
WHEN others THEN
BEGIN
RETURN ST_Intersection(ST_Buffer(geom_a, 0.000001), ST_Buffer(geom_b, 0.000001));
EXCEPTION
WHEN others THEN
BEGIN
RETURN ST_Intersection(ST_SnapToGrid(geom_a, 0.000001), ST_SnapToGrid(geom_b, 0.000001));
EXCEPTION
WHEN others THEN
BEGIN
RETURN 'MULTIPOLYGON EMPTY';
END;
END;
END;
END;
$function$;
----------------------------------------------------------------------------------------------------------------------------------------
-- layers.run_overlap(text);
CREATE OR REPLACE PROCEDURE layers.run_overlap(layer_table text DEFAULT NULL::text)
LANGUAGE plpgsql
AS $procedure$
DECLARE
l_id int;
l_table text;
l_type text;
debug_log text;
BEGIN
IF layer_table IS NOT NULL THEN
RAISE NOTICE 'Running overlap for layer: % ...', layer_table;
RAISE NOTICE 'Check if layers temp table exists...';
DROP TABLE IF EXISTS layers;
RAISE NOTICE '----------------------';
CREATE TEMP TABLE layers ON COMMIT DROP AS (
SELECT layer_id, layer_table_name, "type"
FROM layers."_config"
WHERE overlap IS TRUE
AND layer_table_name = layer_table
);
ELSE
RAISE NOTICE 'Running overlap for ALL layers ...';
RAISE NOTICE 'Check if layers aux table exists...';
DROP TABLE IF EXISTS layers;
RAISE NOTICE '----------------------';
CREATE TEMP TABLE layers ON COMMIT DROP AS (
SELECT layer_id, layer_table_name, "type"
FROM layers."_config"
WHERE overlap IS TRUE
);
END IF;
FOR l_id, l_table, l_type IN SELECT layer_id, layer_table_name, "type" FROM layers
LOOP
IF l_type = 'polygon' THEN
RAISE NOTICE 'Overlaping Layer: % ...', l_table;
-- debug_log =
EXECUTE format($polygon$
INSERT INTO results.overlap_results
(pi_area, layer_id, layer_gid, value, unit, geom, valid)
SELECT i_id, %s::int, l_gid, st_area(st_transform(geom,97823))/ 10000, 'ha'::text,
CASE WHEN st_isvalid(geom) IS false THEN st_collectionextract(st_makevalid(geom),3) ELSE geom END as geom,
-- CASE WHEN coalesce(st_area(st_transform(geom,97823))::numeric / nullif(st_area(st_transform(geom_b,97823))::numeric,0),0)::numeric >= 0.05
-- THEN 1 ELSE 0 END
1
FROM (
SELECT a.gid as i_id, b.gid as l_gid,
CASE WHEN st_within(b.geom,a.geom) THEN b.geom
ELSE st_multi(st_collectionextract(st_safe_intersection(a.geom,b.geom),3))
END as geom, b.geom as geom_b
FROM layers.te_amz_povos_isolados_area a
JOIN layers.%I b ON st_intersects(a.geom,b.geom)
AND NOT st_touches(a.geom,b.geom)
AND NOT st_equals(a.geom,b.geom)
) foo
WHERE NOT EXISTS (
SELECT 1 FROM results.overlap_results c
WHERE c.pi_area = foo.i_id
AND c.layer_id = %s::int
AND c.layer_gid = foo.l_gid
)
$polygon$,
l_id,l_table,l_id);
-- RAISE NOTICE '%', debug_log;
RAISE NOTICE 'Layer: % Done.', l_table;
ELSIF l_type = 'line' THEN null;
RAISE NOTICE 'Overlaping Layer: % ...', l_table;
-- debug_log =
EXECUTE format($line$
INSERT INTO results.overlap_results
(pi_area, layer_id, layer_gid, value, unit, geom)
SELECT i_id, %s::int, l_gid, st_length(st_transform(geom,97823)) / 1000, 'km'::text,
CASE WHEN st_isvalid(geom) IS false THEN st_collectionextract(st_makevalid(geom),2) ELSE geom END as geom
FROM (
SELECT a.gid as i_id, b.gid as l_gid,
CASE WHEN st_within(b.geom,a.geom) THEN b.geom
ELSE st_safe_intersection(a.geom,b.geom)
END as geom
FROM layers.te_amz_povos_isolados_area a
JOIN layers.%I b ON st_intersects(a.geom,b.geom)
AND NOT st_touches(a.geom,b.geom)
) foo
WHERE NOT EXISTS (
SELECT 1 FROM results.overlap_results c
WHERE c.pi_area = foo.i_id
AND c.layer_id = %s::int
AND c.layer_gid = foo.l_gid
)
$line$,
l_id,l_table,l_id);
-- RAISE NOTICE '%', debug_log;
RAISE NOTICE 'Layer: % Done.', l_table;
ELSIF l_type = 'point' THEN
RAISE NOTICE 'Overlaping: % ...', l_table;
-- debug_log =
EXECUTE format($point$
INSERT INTO results.overlap_results
(pi_area, layer_id, layer_gid, value, unit, geom)
SELECT a.gid, %s::int, b.gid, 1::numeric, 'un'::text, b.geom
FROM layers.te_amz_povos_isolados_area a
JOIN layers.%I b ON st_within(b.geom,a.geom)
WHERE NOT EXISTS (
SELECT 1 FROM results.overlap_results c
WHERE c.pi_area = a.gid
AND c.layer_id = %s::int
AND c.layer_gid = b.gid
)
$point$,
l_id, l_table, l_id);
-- RAISE NOTICE '%', debug_log;
RAISE NOTICE 'Layer: % Done.', l_table;
ELSIF l_type = 'raster' THEN null;
--SCRIPT
ELSE RAISE NOTICE 'Geometry Type not identifyied';
END IF;
RAISE NOTICE '----------------------';
END LOOP;
RAISE NOTICE 'Process complete!';
END;
$procedure$
;
----------------------------------------------------------------------------------------------------------------------------------------
layers.run_overlap_update();
CREATE OR REPLACE PROCEDURE layers.run_overlap_update(pi_areas integer[], layer_table text DEFAULT NULL::text)
LANGUAGE plpgsql
AS $procedure$
DECLARE
l_id int;
l_table text;
l_type text;
debug_log text;
BEGIN
IF layer_table IS NOT NULL THEN
RAISE NOTICE 'Running overlap for layer: % ...', layer_table;
RAISE NOTICE 'Check if layers temp table exists...';
DROP TABLE IF EXISTS layers;
RAISE NOTICE '----------------------';
CREATE TEMP TABLE layers ON COMMIT DROP AS (
SELECT layer_id, layer_table_name, "type"
FROM layers."_config"
WHERE overlap IS TRUE
AND layer_table_name = layer_table
);
ELSE
RAISE NOTICE 'Running overlap for ALL layers ...';
RAISE NOTICE 'Check if layers aux table exists...';
DROP TABLE IF EXISTS layers;
RAISE NOTICE '----------------------';
CREATE TEMP TABLE layers ON COMMIT DROP AS (
SELECT layer_id, layer_table_name, "type"
FROM layers."_config"
WHERE overlap IS TRUE
);
END IF;
FOR l_id, l_table, l_type IN SELECT layer_id, layer_table_name, "type" FROM layers
LOOP
IF l_type = 'polygon' THEN
RAISE NOTICE 'Overlaping Layer: % ...', l_table;
-- debug_log =
EXECUTE format($polygon$
INSERT INTO results.overlap_results
(pi_area, layer_id, layer_gid, value, unit, geom, valid)
SELECT i_id, %s::int, l_gid, st_area(st_transform(geom,97823))/ 10000, 'ha'::text,
CASE WHEN st_isvalid(geom) IS false THEN st_collectionextract(st_makevalid(geom),3) ELSE geom END as geom,
-- CASE WHEN coalesce(st_area(st_transform(geom,97823))::numeric / nullif(st_area(st_transform(geom_b,97823))::numeric,0),0)::numeric >= 0.05
-- THEN 1 ELSE 0 END
1
FROM (
SELECT a.gid as i_id, b.gid as l_gid,
CASE WHEN st_within(b.geom,a.geom) THEN b.geom
ELSE st_multi(st_collectionextract(st_safe_intersection(a.geom,b.geom),3))
END as geom, b.geom as geom_b
FROM layers.te_amz_povos_isolados_area a
JOIN layers.%I b ON st_intersects(a.geom,b.geom)
AND NOT st_touches(a.geom,b.geom)
AND NOT st_equals(a.geom,b.geom)
WHERE a.gid = ANY(%L)
) foo
WHERE NOT EXISTS (
SELECT 1 FROM results.overlap_results c
WHERE c.pi_area = foo.i_id
AND c.layer_id = %s::int
AND c.layer_gid = foo.l_gid
)
$polygon$,
l_id,l_table,pi_areas,l_id);
-- RAISE NOTICE '%', debug_log;
RAISE NOTICE 'Layer: % Done.', l_table;
ELSIF l_type = 'line' THEN null;
RAISE NOTICE 'Overlaping Layer: % ...', l_table;
-- debug_log =
EXECUTE format($line$
INSERT INTO results.overlap_results
(pi_area, layer_id, layer_gid, value, unit, geom)
SELECT i_id, %s::int, l_gid, st_length(st_transform(geom,97823)) / 1000, 'km'::text,
CASE WHEN st_isvalid(geom) IS false THEN st_collectionextract(st_makevalid(geom),2) ELSE geom END as geom
FROM (
SELECT a.gid as i_id, b.gid as l_gid,
CASE WHEN st_within(b.geom,a.geom) THEN b.geom
ELSE st_safe_intersection(a.geom,b.geom)
END as geom
FROM layers.te_amz_povos_isolados_area a
JOIN layers.%I b ON st_intersects(a.geom,b.geom)
AND NOT st_touches(a.geom,b.geom)
WHERE a.gid = ANY(%L)
) foo
WHERE NOT EXISTS (
SELECT 1 FROM results.overlap_results c
WHERE c.pi_area = foo.i_id
AND c.layer_id = %s::int
AND c.layer_gid = foo.l_gid
)
$line$,
l_id,l_table,pi_areas,l_id);
-- RAISE NOTICE '%', debug_log;
RAISE NOTICE 'Layer: % Done.', l_table;
ELSIF l_type = 'point' THEN
RAISE NOTICE 'Overlaping: % ...', l_table;
-- debug_log =
EXECUTE format($point$
INSERT INTO results.overlap_results
(pi_area, layer_id, layer_gid, value, unit, geom)
SELECT a.gid, %s::int, b.gid, 1::numeric, 'un'::text, b.geom
FROM layers.te_amz_povos_isolados_area a
JOIN layers.%I b ON st_within(b.geom,a.geom)
WHERE a.gid = ANY(%L)
AND NOT EXISTS (
SELECT 1 FROM results.overlap_results c
WHERE c.pi_area = a.gid
AND c.layer_id = %s::int
AND c.layer_gid = b.gid
)
$point$,
l_id, l_table, pi_areas, l_id);
-- RAISE NOTICE '%', debug_log;
RAISE NOTICE 'Layer: % Done.', l_table;
ELSIF l_type = 'raster' THEN null;
--SCRIPT
ELSE RAISE NOTICE 'Geometry Type not identifyied';
END IF;
RAISE NOTICE '----------------------';
END LOOP;
RAISE NOTICE 'Process complete!';
END;
$procedure$
;
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment