[postgis-users] ST_Value from Polygon

[Andreas Forø Tollefsen]

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>

> If I limit to only 10 polygons it uses 221285 ms.
>
> 2011/2/28 Andreas Forø Tollefsen <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>
>>
>>> 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>
>>>
>>> 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>
>>>>
>>>> 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] *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>
>>>>>
>>>>> 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]
>>>>> *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>
>>>>>
>>>>> 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>
>>>>>
>>>>> 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] *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>
>>>>>
>>>>> 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] *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
>>>>> http://postgis.refractions.net/mailman/listinfo/postgis-users
>>>>>
>>>>>
>>>>>
>>>>>
>>>>> _______________________________________________
>>>>> postgis-users mailing list
>>>>> postgis-users at postgis.refractions.net
>>>>> http://postgis.refractions.net/mailman/listinfo/postgis-users
>>>>>
>>>>>
>>>>>
>>>>>
>>>>>
>>>>>
>>>>>
>>>>> _______________________________________________
>>>>> postgis-users mailing list
>>>>> 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/474ba04e/attachment.html>