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 <andrea...@gmail.com> > 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
