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