[postgis-users] ST_Value from Polygon

[Andreas Forø Tollefsen]

Ok thanks. I got a bit confused regarding the vectorizing. Isnt that what i
am doing with this query? I have now reduced the resolution to increase
performance. Each cell is now 0.02222222222222224x0.02222222222222224.
(previous X 8)


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
INTO globlowrel
FROM (SELECT priogrid_land.gid, ST_Intersection(globlow.rast,
priogrid_land.cell) AS geomval FROM globlow, priogrid_land WHERE
priogrid_land.cell && globlow.rast) AS foo
WHERE (foo).gid >= 139300 AND (foo).gid <= 139301
GROUP BY gid, (foo.geomval).val;

2011/2/28 Pierre Racine <Pierre.Racine at sbf.ulaval.ca>

> 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>
>
> 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/20110301/5f400f9f/attachment.html>