## ## Trying to learn something from here: ## http://profmarcello.blogspot.com/search/label/PostGIS%20Raster ## ## Bottom (past) -> Top(now) worklog chronology -- Using 2nd patch! mygeo=# CREATE TABLE uso_ocupacao_ibitinga_4 AS mygeo-# SELECT mygeo-# gid, mygeo-# (gv).val, mygeo-# (gv).geom mygeo-# FROM ( mygeo(# SELECT mygeo(# gid, mygeo(# ST_Intersection(rast, geom) AS gv mygeo(# FROM mygeo(# municipios, mygeo(# uso_ocupacao_solo mygeo(# WHERE mygeo(# ST_Intersects(rast, geom) mygeo(# AND nome = 'Araraquara' mygeo(# ) mygeo-# AS tmp; SELECT 5966 Time: 177003.414 ms -- Using 1st patch! mygeo=# \timing Timing is on. mygeo=# CREATE TABLE uso_ocupacao_ibitinga_3 AS mygeo-# SELECT mygeo-# gid, mygeo-# (gv).val, mygeo-# (gv).geom mygeo-# FROM ( mygeo(# SELECT mygeo(# gid, mygeo(# ST_Intersection(rast, geom) AS gv mygeo(# FROM mygeo(# municipios, mygeo(# uso_ocupacao_solo mygeo(# WHERE mygeo(# ST_Intersects(rast, geom) mygeo(# AND nome = 'Araraquara' mygeo(# ) mygeo-# AS tmp; SELECT 5966 Time: 296770.638 ms mygeo=# \e CREATE FUNCTION mygeo=# CREATE TABLE uso_ocupacao_ibitinga_2 AS mygeo-# SELECT mygeo-# gid, mygeo-# (gv).val, mygeo-# (gv).geom mygeo-# FROM ( mygeo(# SELECT mygeo(# gid, mygeo(# ST_Intersection(rast, geom) AS gv mygeo(# FROM mygeo(# municipios, mygeo(# uso_ocupacao_solo mygeo(# WHERE mygeo(# ST_Intersects(rast, geom) mygeo(# AND nome = 'Araraquara' mygeo(# ) mygeo-# AS tmp; SELECT 5966 mygeo=# Updating to: as advised here: https://gist.github.com/Komzpa/88784cfa212b2aa82aa37ae261554753 (thanks for your help!) CREATE OR REPLACE FUNCTION _st_intersects(geom geometry, rast raster, nband integer DEFAULT NULL) RETURNS boolean AS $$ DECLARE hasnodata boolean := TRUE; _geom geometry; BEGIN IF ST_SRID(rast) != ST_SRID(geom) THEN RAISE EXCEPTION 'Raster and geometry do not have the same SRID'; END IF; _geom := ST_ConvexHull(rast); IF nband IS NOT NULL THEN SELECT CASE WHEN bmd.nodatavalue IS NULL THEN FALSE ELSE NULL END INTO hasnodata FROM ST_BandMetaData(rast, nband) AS bmd; END IF; IF ST_Intersects(geom, _geom) IS NOT TRUE THEN RETURN FALSE; ELSEIF nband IS NULL OR hasnodata IS FALSE THEN RETURN TRUE; END IF; SELECT ST_Union(t.geom) INTO _geom FROM ST_PixelAsPolygons(rast, nband) AS t; RETURN ST_Intersects(geom, _geom); END; $$ LANGUAGE 'plpgsql' IMMUTABLE PARALLEL SAFE COST 1000; Original function: CREATE OR REPLACE FUNCTION public._st_intersects(geom geometry, rast raster, nband integer DEFAULT NULL::integer) RETURNS boolean LANGUAGE plpgsql IMMUTABLE PARALLEL SAFE COST 1000 AS $function$ DECLARE hasnodata boolean := TRUE; _geom geometry; BEGIN IF ST_SRID(rast) != ST_SRID(geom) THEN RAISE EXCEPTION 'Raster and geometry do not have the same SRID'; END IF; _geom := ST_ConvexHull(rast); IF nband IS NOT NULL THEN SELECT CASE WHEN bmd.nodatavalue IS NULL THEN FALSE ELSE NULL END INTO hasnodata FROM public.ST_BandMetaData(rast, nband) AS bmd; END IF; IF ST_Intersects(geom, _geom) IS NOT TRUE THEN RETURN FALSE; ELSEIF nband IS NULL OR hasnodata IS FALSE THEN RETURN TRUE; END IF; SELECT public.ST_Collect(t.geom) INTO _geom FROM public.ST_PixelAsPolygons(rast, nband) AS t; RETURN public.ST_Intersects(geom, _geom); END; $function$ mygeo=# \df _ST_Intersects List of functions Schema | Name | Result data type | Argument data types | Type --------+----------------+------------------+-----------------------------------------------------------------+-------- public | _st_intersects | boolean | geom1 geometry, geom2 geometry | normal public | _st_intersects | boolean | geom geometry, rast raster, nband integer DEFAULT NULL::integer | normal public | _st_intersects | boolean | rast1 raster, nband1 integer, rast2 raster, nband2 integer | normal (3 rows) mygeo=# \df _ST_Intersects_k List of functions Schema | Name | Result data type | Argument data types | Type --------+------------------+------------------+-----------------------------------------------------------------+-------- public | _st_intersects_k | boolean | geom geometry, rast raster, nband integer DEFAULT NULL::integer | normal (1 row) mygeo=# \df ST_Intersects List of functions Schema | Name | Result data type | Argument data types | Type --------+---------------+------------------+-----------------------------------------------------------------+-------- public | st_intersects | boolean | geography, geography | normal public | st_intersects | boolean | geom1 geometry, geom2 geometry | normal public | st_intersects | boolean | geom geometry, rast raster, nband integer DEFAULT NULL::integer | normal public | st_intersects | boolean | rast1 raster, nband1 integer, rast2 raster, nband2 integer | normal public | st_intersects | boolean | rast1 raster, rast2 raster | normal public | st_intersects | boolean | rast raster, geom geometry, nband integer DEFAULT NULL::integer | normal public | st_intersects | boolean | rast raster, nband integer, geom geometry | normal public | st_intersects | boolean | text, text | normal (8 rows) #running again mygeo=# CREATE TABLE uso_ocupacao_ibitinga_2 AS mygeo-# SELECT mygeo-# gid, mygeo-# (gv).val, mygeo-# (gv).geom mygeo-# FROM ( mygeo(# SELECT mygeo(# gid, mygeo(# ST_Intersection(rast, geom) AS gv mygeo(# FROM mygeo(# municipios, mygeo(# uso_ocupacao_solo mygeo(# WHERE mygeo(# _st_intersects_k(geom,rast) mygeo(# AND nome = 'Araraquara' mygeo(# ) mygeo-# AS tmp; mygeo-# AS tmp; ERROR: GEOSIntersects: TopologyException: side location conflict at -48.37406343009765 -21.794323111743079 CONTEXT: PL/pgSQL function _st_intersects(geometry,raster,integer) line 22 at RETURN PL/pgSQL function st_intersection(geometry,raster,integer) line 5 at assignment SQL function "st_intersection" statement 1 # CREATE OR REPLACE FUNCTION _st_intersects_k(geom geometry, rast raster, nband integer DEFAULT NULL) RETURNS boolean AS $$ DECLARE hasnodata boolean := TRUE; _geom geometry; BEGIN IF ST_SRID(rast) != ST_SRID(geom) THEN RAISE EXCEPTION 'Raster and geometry do not have the same SRID'; END IF; _geom := ST_ConvexHull(rast); IF nband IS NOT NULL THEN SELECT CASE WHEN bmd.nodatavalue IS NULL THEN FALSE ELSE NULL END INTO hasnodata FROM public.ST_BandMetaData(rast, nband) AS bmd; END IF; IF ST_Intersects(geom, _geom) IS NOT TRUE THEN RETURN FALSE; ELSEIF nband IS NULL OR hasnodata IS FALSE THEN RETURN TRUE; END IF; SELECT public.ST_Union(t.geom) INTO _geom FROM public.ST_PixelAsPolygons(rast, nband) AS t; RETURN public.ST_Intersects(geom, _geom); END; $$ LANGUAGE 'plpgsql' IMMUTABLE PARALLEL SAFE COST 1000; CREATE FUNCTION # Creating it to fix.. https://gist.github.com/Komzpa/88784cfa212b2aa82aa37ae261554753 mygeo=# select nome from municipios where st_isvalid(geom) = false; nome ------ (0 rows) https://gis.stackexchange.com/questions/165151/postgis-update-multipolygon-with-st-makevalid-gives-error mygeo=# CREATE TABLE uso_ocupacao_ibitinga_2 AS SELECT gid, (gv).val, (gv).geom FROM ( SELECT gid, ST_Intersection(rast, geom) AS gv FROM municipios, uso_ocupacao_solo WHERE ST_Intersects(rast, geom) AND nome = 'Araraquara' ) AS tmp; ERROR: GEOSIntersects: TopologyException: side location conflict at -48.265729627714023 -21.544296308489137 CONTEXT: PL/pgSQL function _st_intersects(geometry,raster,integer) line 22 at RETURN PL/pgSQL function st_intersection(geometry,raster,integer) line 5 at assignment SQL function "st_intersection" statement 1