[postgis-users] ST_Value from Polygon

[Andreas Forø Tollefsen]

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/0a1c425a/attachment.html>