Re: [postgis-users] How to set up a db for editing in postgis?
Hi, Thanks for the reply. I am aware that I need to set up the spatial_re and geometry columns in order to store spatial information. My question was set in the context of specifically trying to edit the geometry. I have set up a postgresql database with postgis and I can import a shapefile into it. For some reason though I cannot edit the data through wfs-t and I am wondering if it to do with permissions/restrictions etc. If I try to edit the shapefile directly through geoserver it works ok. which permissions and rights should be adjusted to enable editing in the database? yours, Robert Von: Jean-François Gigand An: PostGIS Users Discussion Gesendet: Montag, den 28. Februar 2011, 19:24:30 Uhr Betreff: Re: [postgis-users] How to set up a db for editing in postgis? Hi, To use a PG database with PostGIS, it needs the "geometry_columns" and "spatial_ref_sys" tables, and a set of PostGIS functions. To set these up, you need to run the SQL files : /usr/share/postgresql/8.4/contrib/postgis-1.5/postgis.sql /usr/share/postgresql/8.4/contrib/postgis-1.5/spatial_ref_sys.sql (the path may vary, if your PG version is different - these are those of my Debian install). You need to run these SQL files on your database, as an administrator user. JF 2011/2/28 Robert Buckley : > Hi, > > I have been trying to edit with a postgis databank for ages now and can´t > seem to get it working. > > I have geoserver setup in tomcat6 on ubuntu 10.04. > I have postgres8.4 running with postgis. > I can connect to postgis through pgadmin > I have set up an editing web interface with geoext/openlayers > > I have tried editing a geoserver layer from a shapefile and it works ok!!! > I had to play with the permissions and found that I could only edit when the > file permissions were 0777 > > To edit in postgis I first created a db in postgis called "zgb" the owner > is my postgresql user"geoadmin1" > The tables in postgis give all privilages to public > In geoserver the files are owned by "tomcat6" except "www" which is owned by > "geoadmin1" - this is so that I can´t edit the .js and.html in winscp > > I am really confused as to how to set permissions to be able to edit a > postgisdb. How should the privalages be set to be able to edit geometries? > > If it turns ot that this is not really an issue, what else could be > restricting the editing ability of my app? > > What should I generally watch out for to be able to edit postgis? > > thanks for any help, > > Robert > > > ___ > 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
Re: [postgis-users] ST_Value from Polygon
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-boun...@postgis.refractions.net [mailto:postgis-users-boun...@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:andrea...@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.077160617284 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 mailto:andrea...@gmail.com>> If I limit to only 10 polygons it uses 221285 ms. 2011/2/28 Andreas Forø Tollefsen mailto: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.077160617284 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 mailto: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.78027 139358;0.000563271604818283;0.225308641927313;72.999844495 139358;0.0966820987653989;38.6728395061596;12529.57 139358;0.152631172839676;61.0524691358705;19781.000221 2011/2/25 Andreas Forø Tollefsen mailto:andrea...@gmail.com>> Hi P
Re: [postgis-users] How to set up a db for editing in postgis?
Hi, To use a PG database with PostGIS, it needs the "geometry_columns" and "spatial_ref_sys" tables, and a set of PostGIS functions. To set these up, you need to run the SQL files : /usr/share/postgresql/8.4/contrib/postgis-1.5/postgis.sql /usr/share/postgresql/8.4/contrib/postgis-1.5/spatial_ref_sys.sql (the path may vary, if your PG version is different - these are those of my Debian install). You need to run these SQL files on your database, as an administrator user. JF 2011/2/28 Robert Buckley : > Hi, > > I have been trying to edit with a postgis databank for ages now and can´t > seem to get it working. > > I have geoserver setup in tomcat6 on ubuntu 10.04. > I have postgres8.4 running with postgis. > I can connect to postgis through pgadmin > I have set up an editing web interface with geoext/openlayers > > I have tried editing a geoserver layer from a shapefile and it works ok!!! > I had to play with the permissions and found that I could only edit when the > file permissions were 0777 > > To edit in postgis I first created a db in postgis called "zgb" the owner > is my postgresql user"geoadmin1" > The tables in postgis give all privilages to public > In geoserver the files are owned by "tomcat6" except "www" which is owned by > "geoadmin1" - this is so that I can´t edit the .js and.html in winscp > > I am really confused as to how to set permissions to be able to edit a > postgisdb. How should the privalages be set to be able to edit geometries? > > If it turns ot that this is not really an issue, what else could be > restricting the editing ability of my app? > > What should I generally watch out for to be able to edit postgis? > > thanks for any help, > > Robert > > > ___ > 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] How to set up a db for editing in postgis?
Hi, I have been trying to edit with a postgis databank for ages now and can´t seem to get it working. I have geoserver setup in tomcat6 on ubuntu 10.04. I have postgres8.4 running with postgis. I can connect to postgis through pgadmin I have set up an editing web interface with geoext/openlayers I have tried editing a geoserver layer from a shapefile and it works ok!!! I had to play with the permissions and found that I could only edit when the file permissions were 0777 To edit in postgis I first created a db in postgis called "zgb" the owner is my postgresql user"geoadmin1" The tables in postgis give all privilages to public In geoserver the files are owned by "tomcat6" except "www" which is owned by "geoadmin1" - this is so that I can´t edit the .js and.html in winscp I am really confused as to how to set permissions to be able to edit a postgisdb. How should the privalages be set to be able to edit geometries? If it turns ot that this is not really an issue, what else could be restricting the editing ability of my app? What should I generally watch out for to be able to edit postgis? thanks for any help, Robert ___ postgis-users mailing list postgis-users@postgis.refractions.net http://postgis.refractions.net/mailman/listinfo/postgis-users
Re: [postgis-users] Infinite loop in st_intersects - because of incorrect data out of st_transform?
On 28/02/11 15:18, Magnus Hagander wrote:
Hi!
Running the following query locks up postgis completely (in
geos::algorithm::RobustDeterminant):
SELECT st_intersects(somegeometry,
'010320E61001000500737979F3DDCC2CC0F92154F9E7534540F07FF07F8F806E993F7E55C0304B29FFEA8554400634E8D1DD424540B5FEE6A37FCD4540737979F3DDCC2CC0F92154F9E7534540'::geometry)
I believe this is because there are infinite values in that geometry:
# select
ST_AsText('010320E61001000500737979F3DDCC2CC0F92154F9E7534540F07FF07F8F806E993F7E55C0304B29FFEA8554400634E8D1DD424540B5FEE6A37FCD4540737979F3DDCC2CC0F92154F9E7534540'::geometry);
st_astext
---
POLYGON((-14.4001308522972 42.6555167828373,inf inf,-85.9726317957995
82.0924680617579,42.5223944076352 43.6054577711015,-14.4001308522972
42.6555167828373))
(1 row)
ISTM that this should either be rejected as an invalid geometry, or at
least not hang
Hi Magnus,
H - I can definitely reproduce this on trunk with GEOS 3.2 series.
The backtrace from inside GEOS looks like this:
(gdb) bt
#0 0x7fb17c672525 in
geos::algorithm::RobustDeterminant::signOfDet2x2
(x1=-nan(0x8), y1=-nan(0x8),
x2=-nan(0x8), y2=-nan(0x8))
at RobustDeterminant.cpp:175
#1 0x7fb17c6670a0 in geos::algorithm::CGAlgorithms::isCCW
(ring=) at CGAlgorithms.cpp:163
#2 0x7fb17c6a0dae in geos::geomgraph::GeometryGraph::addPolygonRing
(this=0x1a27970, lr=0x1a26560, cwLeft=2, cwRight=0) at GeometryGraph.cpp:262
#3 0x7fb17c6a0f8e in geos::geomgraph::GeometryGraph::addPolygon
(this=0x1a27970, p=0x1a265e0) at GeometryGraph.cpp:288
#4 0x7fb17c6a16b5 in geos::geomgraph::GeometryGraph::add
(this=0x1a27970, g=0x1a265e0) at GeometryGraph.cpp:183
#5 0x7fb17c6a1a54 in GeometryGraph (this=0x1a27970, newArgIndex=1,
newParentGeom=0x1a265e0, bnr=...) at GeometryGraph.cpp:522
#6 0x7fb17c6d74bb in GeometryGraphOperation (this=out>, g0=0x1a26480, g1=0x1a265e0) at GeometryGraphOperation.cpp:59
#7 0x7fb17c701e16 in RelateOp (this=0xfff8, g0=0xfff8,
g1=0x400) at RelateOp.cpp:56
#8 0x7fb17c701e7f in geos::operation::relate::RelateOp::relate
(a=, b=) at RelateOp.cpp:42
#9 0x7fb17c67e549 in geos::geom::Geometry::intersects
(this=0x1a26480, g=0x1a265e0) at Geometry.cpp:371
#10 0x7fb17cefa047 in GEOSIntersects_r (extHandle=0x1a1d480,
g1=0x7ff, g2=0x400) at geos_ts_c.cpp:284
Obviously the robust determinant is not as robust as it's name would
suggest ;) AFAICT taking a quick look at the GEOS code, I don't think
it is able to detect the NaNs within x1/y1/x2/y2 as an exit condition
for the main loop, but someone more familiar with the GEOS code should
verify this.
Now, this geometry actually comes back as the result of st_transform,
which suggests to me that there is *also* a bug in st_transform, which
should never output an invalid geometry. It is the result of:
st_transform('010320CD0B0100050099E6673CA2FC2DC1AD7BF5ED45CC534199E6673CA2FC2DC1D7BDFAF6D28960415A06E670D7E94B41D7BDFAF6D28960415A06E670D7E94B41AD7BF5ED45CC534199E6673CA2FC2DC1AD7BF5ED45CC5341',
4326)
And that geometry is valid before it goes into st_transform, from what
I can tell:
# select
st_astext('010320CD0B0100050099E6673CA2FC2DC1AD7BF5ED45CC534199E6673CA2FC2DC1D7BDFAF6D28960415A06E670D7E94B41D7BDFAF6D28960415A06E670D7E94B41AD7BF5ED45CC534199E6673CA2FC2DC1AD7BF5ED45CC5341'::geometry);
st_astext
--
--
POLYGON((-982609.1179802 5189911.7181081,-982609.1179802
8670871.7181081,3658670.8820198 8670871.7181081,3658670.8820198
5189911.
7181081,-982609.1179802 5189911.7181081))
(1 row)
To me this looks like bugs "below" postgis - it's certainly below the
PostgreSQL level that I know well :-) Not sure if it can be worked
around/fixed at the postgis level or if it needs to be done in the
lower level libraries, but either way it shouldn't behave this way :-)
Yeah. The second question here is how did the (Inf Inf) point get into
PostGIS anyway, which I think is a separate bug related to this.
Looking here
(http://trac.osgeo.org/postgis/browser/trunk/postgis/lwgeom_transform.c)
at transform_point(), if PROJ.4 returns an error then *pj_errno_ref
should be set to non-zero but even though the projection result is a
point at infinity, the error code is still returning success.
I think the key thing to knowing whether this part is a PROJ.4 bug or a
PostGIS bug would be to clarify whether PROJ.4's *pj_errno_ref should be
non-zero if a projected point appears at infinity - if yes, it's a
PROJ.4 bug and
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:andrea...@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.077160617284 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 mailto:andrea...@gmail.com>> If I limit to only 10 polygons it uses 221285 ms. 2011/2/28 Andreas Forø Tollefsen mailto: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.077160617284 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 mailto: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.78027 139358;0.000563271604818283;0.225308641927313;72.999844495 139358;0.0966820987653989;38.6728395061596;12529.57 139358;0.152631172839676;61.0524691358705;19781.000221 2011/2/25 Andreas Forø Tollefsen mailto: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 mailto: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 i
[postgis-users] Infinite loop in st_intersects - because of incorrect data out of st_transform?
Hi!
Running the following query locks up postgis completely (in
geos::algorithm::RobustDeterminant):
SELECT st_intersects(somegeometry,
'010320E61001000500737979F3DDCC2CC0F92154F9E7534540F07FF07F8F806E993F7E55C0304B29FFEA8554400634E8D1DD424540B5FEE6A37FCD4540737979F3DDCC2CC0F92154F9E7534540'::geometry)
I believe this is because there are infinite values in that geometry:
# select
ST_AsText('010320E61001000500737979F3DDCC2CC0F92154F9E7534540F07FF07F8F806E993F7E55C0304B29FFEA8554400634E8D1DD424540B5FEE6A37FCD4540737979F3DDCC2CC0F92154F9E7534540'::geometry);
st_astext
---
POLYGON((-14.4001308522972 42.6555167828373,inf inf,-85.9726317957995
82.0924680617579,42.5223944076352 43.6054577711015,-14.4001308522972
42.6555167828373))
(1 row)
ISTM that this should either be rejected as an invalid geometry, or at
least not hang
Now, this geometry actually comes back as the result of st_transform,
which suggests to me that there is *also* a bug in st_transform, which
should never output an invalid geometry. It is the result of:
st_transform('010320CD0B0100050099E6673CA2FC2DC1AD7BF5ED45CC534199E6673CA2FC2DC1D7BDFAF6D28960415A06E670D7E94B41D7BDFAF6D28960415A06E670D7E94B41AD7BF5ED45CC534199E6673CA2FC2DC1AD7BF5ED45CC5341',
4326)
And that geometry is valid before it goes into st_transform, from what
I can tell:
# select
st_astext('010320CD0B0100050099E6673CA2FC2DC1AD7BF5ED45CC534199E6673CA2FC2DC1D7BDFAF6D28960415A06E670D7E94B41D7BDFAF6D28960415A06E670D7E94B41AD7BF5ED45CC534199E6673CA2FC2DC1AD7BF5ED45CC5341'::geometry);
st_astext
--
--
POLYGON((-982609.1179802 5189911.7181081,-982609.1179802
8670871.7181081,3658670.8820198 8670871.7181081,3658670.8820198
5189911.
7181081,-982609.1179802 5189911.7181081))
(1 row)
To me this looks like bugs "below" postgis - it's certainly below the
PostgreSQL level that I know well :-) Not sure if it can be worked
around/fixed at the postgis level or if it needs to be done in the
lower level libraries, but either way it shouldn't behave this way :-)
--
Magnus Hagander
Me: http://www.hagander.net/
Work: http://www.redpill-linpro.com/
___
postgis-users mailing list
postgis-users@postgis.refractions.net
http://postgis.refractions.net/mailman/listinfo/postgis-users
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.077160617284 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 > If I limit to only 10 polygons it uses 221285 ms. > > 2011/2/28 Andreas Forø Tollefsen > > 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.077160617284 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 >> >>> 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.78027 >>> 139358;0.000563271604818283;0.225308641927313;72.999844495 >>> 139358;0.0966820987653989;38.6728395061596;12529.57 >>> 139358;0.152631172839676;61.0524691358705;19781.000221 >>> >>> >>> 2011/2/25 Andreas Forø Tollefsen >>> >>> 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 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 > > 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
Re: [postgis-users] ST_Value from Polygon
If I limit to only 10 polygons it uses 221285 ms. 2011/2/28 Andreas Forø Tollefsen > 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.077160617284 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 > >> 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.78027 >> 139358;0.000563271604818283;0.225308641927313;72.999844495 >> 139358;0.0966820987653989;38.6728395061596;12529.57 >> 139358;0.152631172839676;61.0524691358705;19781.000221 >> >> >> 2011/2/25 Andreas Forø Tollefsen >> >> 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 >>> >>> 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 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
Re: [postgis-users] ST_Value from Polygon
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.077160617284 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 > 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.78027 > 139358;0.000563271604818283;0.225308641927313;72.999844495 > 139358;0.0966820987653989;38.6728395061596;12529.57 > 139358;0.152631172839676;61.0524691358705;19781.000221 > > > 2011/2/25 Andreas Forø Tollefsen > > 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 >> >> 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 >>> >>> 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?
