[postgis-users] ST_Value from Polygon

Mon Feb 28 12:34:33 PST 2011

By reducing the tile size from 100x100 to 50x50 I get much better performance. But still vectorizing is faster than counting pixels...

From: postgis-users-bounces at postgis.refractions.net [mailto:postgis-users-bounces at postgis.refractions.net] On Behalf Of Pierre Racine
Sent: 28 février 2011 10:32
To: Andreas Forø Tollefsen; PostGIS Users Discussion
Subject: Re: [postgis-users] ST_Value from Polygon

Andreas,

Actually, since your raster coverage is tiled you should have used ST_Intersects(priogrid_land.cell, globshort.rast) from the beginning.

Understand that to do this query the system has to determine the value for each single pixel (they are 7 231 680 000 !!!!). Having precomputed statistics, like the histogram for each tile, would not help as your polygon grid do not necessarily correspond to the tile grid. Does it?

The real solution, I think, lies in using a lower resolution (overview) of the coverage to compute the statistics. The resulting numbers would be a bit different but not really far from the one you get at highest resolution. You would have to import your raster using the -l option and launch your query on the overview table instead.

I'm looking for a solution actually counting pixel values but it is still slower than using ST_Intersection() (to my own surprise!). That means, for now, it is faster to convert each tile to a vector representation (that'a what ST_Intersection() is doing) and compute the areas of those polygons than to count and summarize pixel values. Work in progress...

Pierre

From: Andreas Forø Tollefsen [mailto:andreasft at gmail.com]
Sent: 28 février 2011 04:48
To: PostGIS Users Discussion
Cc: Pierre Racine; Paragon Corporation
Subject: Re: [postgis-users] ST_Value from Polygon

I have reduced search time. It takes 534402 ms for 100 polygons now by using a && operator in the intersection query.
Having 64818 polygons this would take about 96 hours to complete.

DROP TABLE IF EXISTS globshortrel;

SELECT
gid,
((foo.geomval).val),
CAST (SUM(ST_Area((foo.geomval).geom)) AS decimal(8,5)) as area,
CAST(SUM(ST_Area((foo.geomval).geom))/0.25*100 AS decimal(6,4)) as percentarea,
CAST (SUM(ST_Area((foo.geomval).geom))/0.0000077160617284 AS int) AS npixels
INTO globshortrel
FROM (SELECT priogrid_land.gid, ST_Intersection(globshort.rast, priogrid_land.cell) AS geomval FROM globshort, priogrid_land WHERE priogrid_land.cell && globshort.rast) AS foo
WHERE gid >= 139300 AND gid <= 139399
GROUP BY gid, (foo.geomval).val;

Total query runtime: 534402 ms.

2011/2/28 Andreas Forø Tollefsen <andreasft at gmail.com<mailto:andreasft at gmail.com>>
If I limit to only 10 polygons it uses 221285 ms.
2011/2/28 Andreas Forø Tollefsen <andreasft at gmail.com<mailto:andreasft at gmail.com>>

Had the query running for over 48 hours but now it had crashed. Any idea on how to optimize for speed and stability?

DROP TABLE IF EXISTS globshortrel;

SELECT
gid,
((foo.geomval).val),
CAST (SUM(ST_Area((foo.geomval).geom)) AS decimal(7,5)) as area,
CAST(SUM(ST_Area((foo.geomval).geom))/0.25*100 AS decimal(6,4)) as percentarea,
CAST (SUM(ST_Area((foo.geomval).geom))/0.0000077160617284 AS int) AS npixels
INTO globshortrel
FROM (SELECT globshort.rid, priogrid_land.cell, priogrid_land.gid, ST_Intersection(globshort.rast, priogrid_land.cell) AS geomval FROM globshort, priogrid_land) AS foo
GROUP BY gid, (foo.geomval).val
ORDER BY gid;

-- SELECT * FROM globshortrel;

-- SELECT DISTINCT ON(gid) gid, val, percentarea
-- FROM globshortrel
-- GROUP BY gid, val, percentarea
-- ORDER BY gid, percentarea DESC
ERROR:  out of memory
DETAIL:  Failed on request of size 1753647.
CONTEXT:  SQL function "st_dumpaspolygons" statement 1
PL/pgSQL function "st_intersection" line 9 at RETURN QUERY
SQL function "st_intersection" statement 1

********** Error **********

ERROR: out of memory
SQL state: 53200
Detail: Failed on request of size 1753647.
Context: SQL function "st_dumpaspolygons" statement 1
PL/pgSQL function "st_intersection" line 9 at RETURN QUERY
SQL function "st_intersection" statement 1

2011/2/25 Andreas Forø Tollefsen <andreasft at gmail.com<mailto:andreasft at gmail.com>>
I wrote this query. It seems to give me the area, percentage of area and the pixelcount within a 0.5x0.5 cell.

SELECT gid, SUM(ST_Area((foo.geomval).geom)) as area, SUM(ST_Area((foo.geomval).geom))/0.25*100 as percentarea, SUM(ST_Area((foo.geomval).geom))/7.716049382716049382716049382716e-6 AS Npixels
FROM (SELECT onecell.rid, priogrid_land.cell, priogrid_land.gid, ST_Intersection(onecell.rast, priogrid_land.cell) AS geomval FROM onecell, priogrid_land) AS foo
WHERE gid = 139358
GROUP BY gid, ((foo.geomval).val)
ORDER BY area, ((foo.geomval).val)

Result:
gid; area; percentarea; npixels
139358;0.000123456790106502;0.0493827160426008;15.9999999978027
139358;0.000563271604818283;0.225308641927313;72.9999999844495
139358;0.0966820987653989;38.6728395061596;12529.9999999957
139358;0.152631172839676;61.0524691358705;19781.0000000221

2011/2/25 Andreas Forø Tollefsen <andreasft at gmail.com<mailto:andreasft at gmail.com>>

Hi Pierre. I did tile it when i processed it.
I used this gdal2wktraster syntax: python gdal2wktraster.py -r c:\prio_grid\source\globcover\globshort.tif -t globshort -s 4326 -k 720x360 -I -o globshort.sql

The size of my raster is X: 129600 Y: 55800
2011/2/24 Pierre Racine <Pierre.Racine at sbf.ulaval.ca<mailto:Pierre.Racine at sbf.ulaval.ca>>

You should just divide the area of your polygon by the area of one of your pixel.

I want to investigate this example a little bit further. What is the size of your raster in pixels? I understand that you did not tile it. Did you?

Pierre

From: postgis-users-bounces at postgis.refractions.net<mailto:postgis-users-bounces at postgis.refractions.net> [mailto:postgis-users-bounces at postgis.refractions.net<mailto:postgis-users-bounces at postgis.refractions.net>] On Behalf Of Andreas Forø Tollefsen
Sent: 24 février 2011 09:53
To: Paragon Corporation
Cc: PostGIS Users Discussion

Subject: Re: [postgis-users] ST_Value from Polygon

Great.
Thank you so much. I should have noticed that these were unioned numbers.
Btw. are there any way of getting the count of pixels within a polygon without actually aggregating them or union them?

Andreas
2011/2/24 Paragon Corporation <lr at pcorp.us<mailto:lr at pcorp.us>>
Andreas,
Sorry should have recognized what you're doing.  The intersection returns a polygon which is a union of the clipped raster pixel squares.  So you need to use Sum of area instead and then divide by the area of a pixel to get the equivalent of your count.

So

SELECT gid, SUM(ST_Area((foo.geomval).geom))/ [put your pixel area size here]  as ct
FROM (SELECT globshort.rid, priogrid_land.cell, priogrid_land.gid, ST_Intersection(globshort.rast, priogrid_land.cell) AS geomval FROM globshort, priogrid_land) AS foo
WHERE gid >= 139358 AND gid <= 139365
GROUP BY gid
ORDER BY gid

________________________________
From: Andreas Forø Tollefsen [mailto:andreasft at gmail.com<mailto:andreasft at gmail.com>]
Sent: Thursday, February 24, 2011 8:33 AM
To: PostGIS Users Discussion
Cc: Paragon Corporation

Subject: Re: [postgis-users] ST_Value from Polygon

I am a bit unsure whether my results are actually correct. According to a total count using the below query, I get very different results between the cells.
Since the raster does actually cover the whole vector cell, i would assume that the count should be similar in all cells. Meaning, the pixel count should be the same.
What i get is different, and it seems that the query is not providing me with the number of pixels within the grid cell.
Any idea why this is so different?

SELECT gid, count((foo.geomval).val) as ct
FROM (SELECT globshort.rid, priogrid_land.cell, priogrid_land.gid, ST_Intersection(globshort.rast, priogrid_land.cell) AS geomval FROM globshort, priogrid_land) AS foo
WHERE gid >= 139358 AND gid <= 139365
GROUP BY gid
ORDER BY gid

Result:
139358;632
139359;1030
139360;912
139361;731
139362;760
139363;1230
139364;1314
139365;1014

The attached image shows the raster pixels within one cell.

2011/2/24 Andreas Forø Tollefsen <andreasft at gmail.com<mailto:andreasft at gmail.com>>
Thanks!
That solved it.

This will probably take a lot of time. I have 259200 polygons measuring 0.5 x 0.5 decimal degrees while the raster dataset is of global cover and has a pixelsize of 0.00277777777777778x0.00277777777777778.

Andreas

2011/2/23 Paragon Corporation <lr at pcorp.us<mailto:lr at pcorp.us>>
Andrea,

Try

SELECT DISTINCT ON(gid) gid, (foo.geomval).val, COUNT((foo.geomval).val) AS ct
FROM (SELECT globshort.rid, priogrid_land.cell, priogrid_land.gid, ST_Intersection(globshort.rast, priogrid_land.cell) AS geomval FROM globshort, priogrid_land) AS foo
WHERE gid > 151000 AND gid < 151010
GROUP BY gid, (foo.geomval).val
ORDER BY gid, ct DESC

________________________________
From: postgis-users-bounces at postgis.refractions.net<mailto:postgis-users-bounces at postgis.refractions.net> [mailto:postgis-users-bounces at postgis.refractions.net<mailto:postgis-users-bounces at postgis.refractions.net>] On Behalf Of Andreas Forø Tollefsen
Sent: Wednesday, February 23, 2011 4:05 AM
To: PostGIS Users Discussion
Subject: Re: [postgis-users] ST_Value from Polygon
Hi. Thanks Regina and Leo,
I have been testing the raster and geom intersection a bit. I guess what i need is to use the ST_Intersection together with a max(count) function.
So my result will be the rastervalue with the highest count within each of the grid cells.
However, as far as i know, there is now Max(COUNT) function in postgresql.

Any idea how i can modify the below query to only return the rastervalue within the grid cell occuring most frequently?
Consequently i want only one row for each gid, and the maximum occuring rastervalue.

SELECT gid, (foo.geomval).val, COUNT((foo.geomval).val) AS ct
FROM (SELECT globshort.rid, priogrid_land.cell, priogrid_land.gid, ST_Intersection(globshort.rast, priogrid_land.cell) AS geomval FROM globshort, priogrid_land) AS foo
WHERE gid > 151000 AND gid < 151010
GROUP BY gid, (foo.geomval).val;

gid; val; ct
151001;14;381
151001;150;9
151001;50;7
151001;140;91
151001;40;1
151001;70;2
151001;130;4
151001;200;48
151001;100;3
151001;;0
151001;190;1
151001;20;203
151001;11;111
151001;210;16
151001;30;105

2011/2/23 Paragon Corporation <lr at pcorp.us<mailto:lr at pcorp.us>>
Have you looked at ST_Intersection.  I'm not sure how large your grids are so might still be a bit too slow.

http://www.postgis.org/documentation/manual-svn/RT_ST_Intersection.html

Below is a link to our slides from our North Carolina GIS meeting that may answer some of your questions (shows some Raster examples) as well as the 3D ones people have asked.

http://www.postgis.us/presentations

Hope that helps,
Regina and Leo
________________________________
From: postgis-users-bounces at postgis.refractions.net<mailto:postgis-users-bounces at postgis.refractions.net> [mailto:postgis-users-bounces at postgis.refractions.net<mailto:postgis-users-bounces at postgis.refractions.net>] On Behalf Of Andreas Forø Tollefsen
Sent: Tuesday, February 22, 2011 4:28 AM
To: PostGIS Users Discussion
Subject: [postgis-users] ST_Value from Polygon
Hi all,

I am working with a large raster dataset that i want to aggregate into vector grids.
The raster dataset is a landcover dataset, and i want to find which of the raster values are the most dominant within each of the vector grid cells.

I have been looking at the ST_Value function, but this is not usable together with the cell polygon.

I have written a script that gives me the raster value of the centroid of each cell, but i want to find which raster class is the largest.
Hence i need to calculate the area of each raster class within each cell and select the largest class.

Any idea? So far i have only come this far:

DROP TABLE IF EXISTS globshortpoly;
SELECT priogrid_land.cell, ST_Value(rast, ST_Centroid(cell))
INTO globshortpoly
FROM priogrid_land, globshort
WHERE rast && priogrid_land.cell
LIMIT 1000

_______________________________________________
postgis-users mailing list
postgis-users at postgis.refractions.net<mailto:postgis-users at postgis.refractions.net>
http://postgis.refractions.net/mailman/listinfo/postgis-users

_______________________________________________
postgis-users mailing list
postgis-users at postgis.refractions.net<mailto:postgis-users at postgis.refractions.net>
http://postgis.refractions.net/mailman/listinfo/postgis-users

_______________________________________________
postgis-users mailing list
postgis-users at postgis.refractions.net<mailto:postgis-users at postgis.refractions.net>
http://postgis.refractions.net/mailman/listinfo/postgis-users

-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.osgeo.org/pipermail/postgis-users/attachments/20110228/5a2681d5/attachment.html>