[postgis-users] ST_Value from Polygon

Fri Feb 25 00:49:13 PST 2011

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/20110225/7751709f/attachment.html>