Learning_raster_and_geom_intersect.worklog.log

##
## 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

Using patch1!

Using patch2!

Testing classification! Works pretty nice!

https://github.com/Komzpa, thank you very much for all your help! I learned a lot.
https://trac.osgeo.org/postgis/ticket/4132