If I limit to only 10 polygons it uses 221285 ms. 2011/2/28 Andreas Forø Tollefsen <andrea...@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 <andrea...@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 <andrea...@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.rac...@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-boun...@postgis.refractions.net [mailto: >>>> postgis-users-boun...@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 <l...@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:andrea...@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 <andrea...@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 <l...@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-boun...@postgis.refractions.net [mailto: >>>> postgis-users-boun...@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 <l...@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-boun...@postgis.refractions.net [mailto: >>>> postgis-users-boun...@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@postgis.refractions.net >>>> http://postgis.refractions.net/mailman/listinfo/postgis-users >>>> >>>> >>>> >>>> >>>> _______________________________________________ >>>> postgis-users mailing list >>>> postgis-users@postgis.refractions.net >>>> http://postgis.refractions.net/mailman/listinfo/postgis-users >>>> >>>> >>>> >>>> >>>> >>>> >>>> >>>> _______________________________________________ >>>> postgis-users mailing list >>>> postgis-users@postgis.refractions.net >>>> http://postgis.refractions.net/mailman/listinfo/postgis-users >>>> >>>> >>> >> >
_______________________________________________ postgis-users mailing list postgis-users@postgis.refractions.net http://postgis.refractions.net/mailman/listinfo/postgis-users
