[postgis-users] ST_Value from Polygon

Mon Feb 28 00:59:38 PST 2011

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