[postgis-users] Questions related to some tests I did in PostGIS

[DanielFranco]

Hi, I'm kind of investigating some PostGIS functions to apply in one of the
projects i'm being charged of
and I would like to understand, if possible, how these functions really
works inside the database. The story is very long, so I hope I'm not
bothering you people.

About the project, basically, it uses a land coverage shapefile (generated
by geoprocessing a raster) and 
it calculates the intersection between the shapefile and some polygons. So,
what I'm trying to do is to intersect the polygons with the raster used to
generate the shapefile in order to gain some performance, since the
shapefile would not be required anymore. The raster size is about
40000x80000 and it was tiled in database to 100x100 pixels each tile (pixel
scale is 5 in x and y axis).

As I studied the PostGIS Raster functions, I realized that the area
calculations by intersecting a polygon to the raster is done by
vectorizating the raster in query time and then doing the intersection,
which showed me that working with raster format is not suitable for this
case, since I could use the shapefile as a "pre-processing stage" and then
run queries within it, besides that, the shapefile contains some fields
(obtained by geoprocessing the raster) that I could not get by simply
loading the raster and I would have to analyze it to get the same effect.

I run queries with 2 differents subsets of polygons. One of them is the
result of the
"SELECT geom FROM vector WHERE gridcode = 18" query and the other one is a
simple select of 
the geometry of an unique big polygon. The first subset has a bigger
intersection area.

/* 1st subsect */

-- vector/vector approach --
SELECT sum(ST_Area(ST_Intersection(V.geom,V2.geom)))
FROM 
	(SELECT geom FROM vector WHERE gridcode = 18) as V
	INNER JOIN 
	mosaic_vector as V2
	ON (ST_Intersects(V.geom,V2.geom)) -- 49610419.7890016 -- 18 046 ms --

--raster/vector approach --
SELECT  sum(ST_Area((gv).geom)) 
FROM
	(SELECT (ST_Intersection(rast,1,geom)) as gv
	FROM 	(SELECT geom FROM vector WHERE gridcode = 18) as V
		INNER JOIN 
		mosaic as R
		ON ST_Intersects(R.rast, V.geom)
	) as intersection -- 49 610 419.7890017 -- 129 985 ms --

/* 2nd subsect */

-- vector/vector approach --
SELECT sum(ST_Area(ST_Intersection(V.geom,V2.geom)))
FROM 
	(SELECT geom FROM polygon) as V
	INNER JOIN 
	mosaic_vector as V2
	ON (ST_Intersects(V.geom,V2.geom)) -- 4710228419.85417 -- 5 613 ms --

--raster/vector approach --
SELECT  sum(ST_Area((gv).geom)) 
FROM
	(SELECT (ST_Intersection(rast,1,geom)) as gv
	FROM 	(SELECT geom FROM polygon) as V
		INNER JOIN 
		mosaic as R
		ON ST_Intersects(R.rast, V.geom)
	) as intersection -- 4 710 228 419.85417 -- 80 696 ms --


Another thing that I tried to do was to count the pixels of the intersection
of the raster and the polygons
(in this case I used the ST_Clip function to clip the raster and then count
the pixels) and as expected
the error of area calculation was big due to ST_Clip problem of
including/not including pixels in the boundaries of the polygon. However,
the query time was faster than the vector/vector approach in the first
subset and not in the second subset for some reason that I don't know how to
explain.

/* Pixels count first subset*/

-- By Dumping as polygons --

SELECT  sum(ST_Area((gv).geom)) 
FROM
	(SELECT (ST_DumpAsPolygons(ST_Clip(rast,1,geom))) as gv
	FROM 	(SELECT geom FROM mosaic_vector WHERE gridcode = 18) as V
		INNER JOIN 
		mosaic as R
		ON ST_Intersects(R.rast, V.geom)) as intersection -- 49 596 500 -- 14 716
ms --

-- By counting pixels --

SELECT (soma * abs(scale_x) * abs(scale_y)) 
FROM
	(SELECT sum(ST_Count(ST_Clip(rast,1,geom))) as soma
	FROM 	(SELECT geom FROM mosaic_vector WHERE gridcode = 18) as V 
		INNER JOIN 
		mosaic as R
		ON ST_Intersects(r.rast, V.geom)
	) as intersection 

	CROSS JOIN

	(SELECT ST_ScaleX(rast) as scale_x, ST_ScaleY(rast) as scale_y 
	FROM mosaic LIMIT 1
	) as scales; -- 49 596 500 -- 11 562 ms --

/* Pixels count second subset */

-- Clip --

SELECT  sum(ST_Area((gv).geom)) 
FROM
	(SELECT (ST_DumpAsPolygons(ST_Clip(rast,1,geom))) as gv
	FROM 	(SELECT geom FROM polygon) as V
		INNER JOIN 
		mosaic as R
		ON ST_Intersects(R.rast, V.geom)) as intersection -- 4 710 221 900 -- 540
761 ms --

-- Pixels --

SELECT (soma * abs(scale_x) * abs(scale_y)) 
FROM
	(SELECT sum(ST_Count(ST_Clip(rast,1,geom))) as soma
	FROM 	(SELECT geom FROM polygon) as V 
		INNER JOIN 
		mosaic as R
		ON ST_Intersects(r.rast, V.geom)
	) as intersection 

	CROSS JOIN

	(SELECT ST_ScaleX(rast) as scale_x, ST_ScaleY(rast) as scale_y 
	FROM mosaic LIMIT 1
	) as scales; -- 4 710 221 900 -- 473 387 ms --


So, I resampled the raster (keeping its real size, but changing the number
of pixels) to make it have pixels in scale 1 in x and y axis in order to get
better precision and faster query time.
By doing that, I've indeed gotten better precision, but the query time
increased a lot by increasing the number of pixels and the vector/vector
approach was better in that aspect.

So, my questions are:
Why counting pixels was not faster in the second subset with pixels of 5x5
meters?

How ST_Count works? Does it uses the ST_Value function in each pixel of the
pixel matrix row by row in each tile to keep track of the nodata value
pixels? 

Is the raster support in postgis, in general, only suitable to get pixel
statistics of some region in the raster?

Whats the advantages of using postgis raster to analyze rasters compared to
some of the arcgis tools (for example)? I ask that because is much more
easier to visualize the data in arcgis compared to postgis tables in QGIS
plugin (wkt raster)

Thanks for letting me question this.



--
View this message in context: http://postgis.17.x6.nabble.com/Questions-related-to-some-tests-I-did-in-PostGIS-tp5007007.html
Sent from the PostGIS - User mailing list archive at Nabble.com.