[postgis-users] ST_Boundary(geometry) for polyon to line conversion?
Hi, ST_Boundary(geometry) Is this the function I would use to convert a polygon to a polyline? yours, Rob___ postgis-users mailing list postgis-users@postgis.refractions.net http://postgis.refractions.net/mailman/listinfo/postgis-users
Re: [postgis-users] Invalid SRID for geography type: ee4
Hi, Geography type supports only EPSG:4326. First load your data as geometry, then change SRS and cast them to geography. Nicolas On 10 April 2012 22:50, Philippe watche...@gmail.com wrote: Hello, I'm running into an issue loading data from the Belgium IGN. I can transform a shapefile into Postgis loadable data with Geometry but not Geography. I can't imagine why. Any ideas ? Thanks The following works /usr/lib/postgresql/9.1/bin/shp2pgsql -WLATIN1 -s3812 -d -I /home/postgres/warehouse/shapes/BE/vector/Shape_L08/commune_polygon.shp be_ign /tmp/be_ign.sql The following fails /usr/lib/postgresql/9.1/bin/shp2pgsql -WLATIN1 -s3812 -G -d -I /home/postgres/warehouse/shapes/BE/vector/Shape_L08/commune_polygon.shp be_ign /tmp/be_ign.sql with Shapefile type: Polygon Postgis type: MULTIPOLYGON[2] Invalid SRID for geography type: ee4 Although 0xee4=3812 which is the EPSG for the Lambert 2008 projection used by that file. And my spatial ref table contains it: public=# select * from spatial_ref_sys where auth_srid=3812; srid | 3812 auth_name | EPSG auth_srid | 3812 srtext | PROJCS[ETRS89 / Belgian Lambert 2008,GEOGCS[ETRS89,DATUM[European_Terrestrial_Reference_System_1989,SPHEROID[GRS 1980,6378137,298.257222101,AUTHORITY[EPSG,7019]],AUTHORITY[EPSG,6258]],PRIMEM[Greenwich,0,AUTHORITY[EPSG,8901]],UNIT[degree,0.01745329251994328,AUTHORITY[EPSG,9122]],AUTHORITY[EPSG,4258]],UNIT[metre,1,AUTHORITY[EPSG,9001]],PROJECTION[Lambert_Conformal_Conic_2SP],PARAMETER[standard_parallel_1,49.84],PARAMETER[standard_parallel_2,51.16],PARAMETER[latitude_of_origin,50.797815],PARAMETER[central_meridian,4.3592158],PARAMETER[false_easting,649328],PARAMETER[false_northing,665262],AUTHORITY[EPSG,3812],AXIS[X,EAST],AXIS[Y,NORTH]] proj4text | +proj=lcc +lat_1=49.84 +lat_2=51.16 +lat_0=50.797815 +lon_0=4.3592158 +x_0=649328 +y_0=665262 +ellps=GRS80 +units=m +no_defs ___ 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_Boundary(geometry) for polyon to line conversion?
Hi, ST_Boundary(geometry) Is this the function I would use to convert a polygon to a polyline? yours, Rob Hi, Yes, but depends on the shape of the polygon. If it has holes, st_boundary will return several linestring, for the exterior and for each hole. If you need only the exterior shape of the polygon, st_exteriorRing may suit your needs better. Nicolas ___ postgis-users mailing list postgis-users@postgis.refractions.net http://postgis.refractions.net/mailman/listinfo/postgis-users
Re: [postgis-users] Simplifying a multipolygon layer, keeping polygon connection
On 10 April 2012 12:04, Sandro Santilli s...@keybit.net wrote: On Mon, Apr 09, 2012 at 07:33:41PM +0200, Nicolas Ribot wrote: Hmm much easier with Postgis Topology... (thanks Brent for the Topology example) By the way, is it the right way to use topology to simplify a layer ? It's one way, there may be multiple ones. Add this new way to to wiki ? Done: http://trac.osgeo.org/postgis/wiki/UsersWikiSimplifyWithTopologyExt It'd be interesting to see some banchmarks of the two methods you reported. Just som rough figures from my test: Using first method, with a plpgsql function and spatial operators: • it takes ~ 7.5 second to process 96 multipolygons containing ~ 47000 vertices. • it takes ~ 9 minutes to process the 36610 communes (=municipalities), containing more than 512000 vertices With Topology: • treating departements took 10.2 seconds on average. • I canceled the CreateTopoGeo after about 25 minutes. I'd point out that the toTopoGeom function does the association with you, directly giving you a TopoGeometry object back. Brent didn't use TopoGeometry at all, which is legit, but that's the intended mean to associate attributes to topologically-defined locations. I've been actually thinking that an ST_Simplify() override taking a TopoGeometry input might be a good idea. You would be simplifying the same up to twice each, but the usage would be simper. I will play a bit with this object. Thank you for the advices. Nicolas --strk; create table new_dept as ( select gid, code_dept, (st_dump(geom)).* from departement ); create index new_dept_geom_gist on new_dept using gist(geom); -- create new empty topology structure select CreateTopology('topo1',2154,0); -- add all departements polygons to topology in one operation as a collection select ST_CreateTopoGeo('topo1',ST_Collect(geom)) from departement; -- create a new topo based on the simplification of existing one -- should not be the right way to do it, but calling ST_ChangeEdgeGeom failed with a simplify linestring select CreateTopology('topo2',2154,0); select ST_CreateTopoGeo('topo2', geom) from ( select ST_Collect(st_simplifyPreserveTopology(geom, 1)) as geom from topo1.edge_data ) as foo; -- associates new simplified faces with points on surface, in the new_dept table alter table new_dept add column simple_geom geometry(POLYGON, 2154); -- retrieves polygons by comparing surfaces (pip is not enough for odd-shaped polygons) with simple_face as ( select st_getFaceGeometry('topo2', face_id) as geom from topo2.face where face_id 0 ) update new_dept d set simple_geom = sf.geom from simple_face sf where st_intersects(d.geom, sf.geom) and st_area(st_intersection(sf.geom, d.geom))/st_area(sf.geom) 0.5; Nicolas On 9 April 2012 01:02, Nicolas Ribot nicolas.ri...@gmail.com wrote: For those interested, I've added an entry to the postgis wiki showing an example of polygon layer simplification: http://trac.osgeo.org/postgis/wiki/UsersWikiSimplifyPreserveTopology Nicolas ___ postgis-users mailing list postgis-users@postgis.refractions.net http://postgis.refractions.net/mailman/listinfo/postgis-users -- ,--o-. | __/ | Delivering high quality PostGIS 2.0 ! | / 2.0 | http://strk.keybit.net - http://vizzuality.com `-o--' ___ 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] Simplifying a multipolygon layer, keeping polygon connection
On Wed, Apr 11, 2012 at 10:01:12AM +0200, Nicolas Ribot wrote: On 10 April 2012 12:04, Sandro Santilli s...@keybit.net wrote: On Mon, Apr 09, 2012 at 07:33:41PM +0200, Nicolas Ribot wrote: Hmm much easier with Postgis Topology... (thanks Brent for the Topology example) By the way, is it the right way to use topology to simplify a layer ? It's one way, there may be multiple ones. Add this new way to to wiki ? Done: http://trac.osgeo.org/postgis/wiki/UsersWikiSimplifyWithTopologyExt Great, thanks. It'd be interesting to see some banchmarks of the two methods you reported. Just som rough figures from my test: Using first method, with a plpgsql function and spatial operators: • it takes ~ 7.5 second to process 96 multipolygons containing ~ 47000 vertices. • it takes ~ 9 minutes to process the 36610 communes (=municipalities), containing more than 512000 vertices With Topology: • treating departements took 10.2 seconds on average. • I canceled the CreateTopoGeo after about 25 minutes. What you mean by treating departments ? Is that the 96 multipolygons with ~47000 vertices ? Can you do the comparison using the same (small) input for both ? I'm interested in speeding up the case using topology. Note that ST_CreateTopoGeo needs to do all in memory which would be inappropriate for large datasets, the toTopoGeom function (and the underlying TopoGeo_addGeometry functions) would instead allow you to populate a topology incrementally. For sure I can tell you that topology isn't meant to be used as a temporary storage. I can imagine, for example, that you want multiple resolutions of your data. You could build the topology once and then generate multiple simplified versions of it, or even simplify on demand with an ST_Simplify(TopoGeometry) [ to be implemented ]. --strk; ___ postgis-users mailing list postgis-users@postgis.refractions.net http://postgis.refractions.net/mailman/listinfo/postgis-users
Re: [postgis-users] PostGIS let me remove an edge shared by a topogeometry
On Tue, Apr 10, 2012 at 09:50:05PM +0200, Jose Carlos Martinez Llario wrote: Dear List, Dont know why PostGIS let me remove an edge which is needed by a topogeometry (see point c) ) This is a bug or it is the expected behavior? ST_RemEdgeModFace (which you're calling) is required to _modify_ one face to take up the space previously occupied by the faces this edge was separating. By convention, when one of the faces is the universal face, the _other_ face is kept, while if both faces are not universal, the face on the _right_ is kept and the face on the _left_ is remved. Now, if you have a TopoGeometry defined by only the face that is removed you get an exception, while if you have a TopoGeometry defined by the face that is modified, or by both faces, you get no exception. At least this is the current intended behavior. Not sure if it makes sense. b) when I try to remove an edge which is needed by a topogeometry but one of the side is the universal polygon (id 0) PostGIS does not throw an error: s6=# select ST_RemEdgeModFace ('t3', 6); NOTICE: Updating next_{right,left}_face of ring edges... NOTICE: Deletion of edge 6 joins faces 2 and 0 st_remedgemodface --- 0 Given the above intended behavior this does indedd sound like a bug. c) The topogeometry is an multiopolygon empty s6=# select st_astext(topogeom::geometry) from t3.parcelas where gid = 2; st_astext MULTIPOLYGON EMPTY I guess the topogeometry is now defined by face 0 ? Try something like: SELECT r.* FROM t3.relation r, t3.parcelas p WHERE r.topogeo_id = id(p.topogeom) AND r.layer_id = layer_id(p.topogeom) AND p.gid = 2; --strk; ,--o-. | __/ |Delivering high quality PostGIS 2.0 ! | / 2.0 |http://strk.keybit.net - http://vizzuality.com `-o--' ___ postgis-users mailing list postgis-users@postgis.refractions.net http://postgis.refractions.net/mailman/listinfo/postgis-users
Re: [postgis-users] PostGIS let me remove an edge shared by a topogeometry
Hi Sandro, thanx for the detailed answer. The query returns 0 rows: s6=# SELECT r.* FROM t3.relation r, t3.parcelas p WHERE r.topogeo_id = id(p.topogeom) AND r.layer_id = layer_id(p.topogeom) AND p.gid = 2; topogeo_id | layer_id | element_id | element_type +--++-- (0 rows) s6=# select * from t3.relation ; topogeo_id | layer_id | element_id | element_type +--++-- 1 |1 | 1 |3 3 |1 | 3 |3 Even though the topogeometry exists: s6=# select gid,(topogeom).* from t3.parcelas; gid | topology_id | layer_id | id | type -+-+--++-- 1 | 30 |1 | 1 |3 2 | 30 |1 | 2 |3 3 | 30 |1 | 3 |3 On 11/04/2012 10:41, Sandro Santilli wrote: On Tue, Apr 10, 2012 at 09:50:05PM +0200, Jose Carlos Martinez Llario wrote: Dear List, Dont know why PostGIS let me remove an edge which is needed by a topogeometry (see point c) ) This is a bug or it is the expected behavior? ST_RemEdgeModFace (which you're calling) is required to _modify_ one face to take up the space previously occupied by the faces this edge was separating. By convention, when one of the faces is the universal face, the _other_ face is kept, while if both faces are not universal, the face on the _right_ is kept and the face on the _left_ is remved. Now, if you have a TopoGeometry defined by only the face that is removed you get an exception, while if you have a TopoGeometry defined by the face that is modified, or by both faces, you get no exception. At least this is the current intended behavior. Not sure if it makes sense. b) when I try to remove an edge which is needed by a topogeometry but one of the side is the universal polygon (id 0) PostGIS does not throw an error: s6=# select ST_RemEdgeModFace ('t3', 6); NOTICE: Updating next_{right,left}_face of ring edges... NOTICE: Deletion of edge 6 joins faces 2 and 0 st_remedgemodface --- 0 Given the above intended behavior this does indedd sound like a bug. c) The topogeometry is an multiopolygon empty s6=# select st_astext(topogeom::geometry) from t3.parcelas where gid = 2; st_astext MULTIPOLYGON EMPTY I guess the topogeometry is now defined by face 0 ? Try something like: SELECT r.* FROM t3.relation r, t3.parcelas p WHERE r.topogeo_id = id(p.topogeom) AND r.layer_id = layer_id(p.topogeom) AND p.gid = 2; --strk; ,--o-. | __/ |Delivering high quality PostGIS 2.0 ! | / 2.0 |http://strk.keybit.net - http://vizzuality.com `-o--' ___ 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] PostGIS let me remove an edge shared by a topogeometry
On Wed, Apr 11, 2012 at 10:57:57AM +0200, Jose Carlos Martinez wrote: Hi Sandro, thanx for the detailed answer. The query returns 0 rows: s6=# SELECT r.* FROM t3.relation r, t3.parcelas p WHERE r.topogeo_id = id(p.topogeom) AND r.layer_id = layer_id(p.topogeom) AND p.gid = 2; topogeo_id | layer_id | element_id | element_type +--++-- (0 rows) s6=# select * from t3.relation ; topogeo_id | layer_id | element_id | element_type +--++-- 1 |1 | 1 |3 3 |1 | 3 |3 Even though the topogeometry exists: s6=# select gid,(topogeom).* from t3.parcelas; gid | topology_id | layer_id | id | type -+-+--++-- 1 | 30 |1 | 1 |3 2 | 30 |1 | 2 |3 3 | 30 |1 | 3 |3 Bad bug indeed. Could you please file it on trac ? If you feel like getting your hands dirtier you could also look at providing a patch for the testcases: topology/test/regress/st_remedgemodface.sql topology/test/regress/st_remedgenewface.sql I'm assuming ST_RemEdgeNewFace is also affected, surely worth adding a testcase for it as well, no matter if it is currently affected or not. --strk; ,--o-. | __/ |Delivering high quality PostGIS 2.0 ! | / 2.0 |http://strk.keybit.net - http://vizzuality.com `-o--' ___ postgis-users mailing list postgis-users@postgis.refractions.net http://postgis.refractions.net/mailman/listinfo/postgis-users
Re: [postgis-users] PostGIS let me remove an edge shared by a topogeometry
On Tue, Apr 10, 2012 at 10:13:07PM +0200, Andrea Peri wrote: I guess you should not write directly on topology schema using create table, and insert directly. He's not writing directly, he's using toTopoGeom for the population and ST_RemEdgeModFace for editing. --strk; ,--o-. | __/ |Delivering high quality PostGIS 2.0 ! | / 2.0 |http://strk.keybit.net - http://vizzuality.com `-o--' ___ postgis-users mailing list postgis-users@postgis.refractions.net http://postgis.refractions.net/mailman/listinfo/postgis-users
Re: [postgis-users] Simplifying a multipolygon layer, keeping polygon connection
On 11 April 2012 10:15, Sandro Santilli s...@keybit.net wrote: On Wed, Apr 11, 2012 at 10:01:12AM +0200, Nicolas Ribot wrote: On 10 April 2012 12:04, Sandro Santilli s...@keybit.net wrote: On Mon, Apr 09, 2012 at 07:33:41PM +0200, Nicolas Ribot wrote: Hmm much easier with Postgis Topology... (thanks Brent for the Topology example) By the way, is it the right way to use topology to simplify a layer ? It's one way, there may be multiple ones. Add this new way to to wiki ? Done: http://trac.osgeo.org/postgis/wiki/UsersWikiSimplifyWithTopologyExt Great, thanks. It'd be interesting to see some banchmarks of the two methods you reported. Just som rough figures from my test: Using first method, with a plpgsql function and spatial operators: • it takes ~ 7.5 second to process 96 multipolygons containing ~ 47000 vertices. • it takes ~ 9 minutes to process the 36610 communes (=municipalities), containing more than 512000 vertices With Topology: • treating departements took 10.2 seconds on average. • I canceled the CreateTopoGeo after about 25 minutes. What you mean by treating departments ? Is that the 96 multipolygons with ~47000 vertices ? Can you do the comparison using the same (small) input for both ? Yes, here are some stats: Mac OS X Snow Leopard, SSD drive, 8 gb ram shared_buffers = 5120MB work_mem = 1024MB POSTGIS=2.0.0 r9605 GEOS=3.3.2-CAPI-1.7.2 PROJ=Rel. 4.7.1, 23 September 2009 GDAL=GDAL 1.8.0, released 2011/01/12 LIBXML=2.7.8 TOPOLOGY RASTER Number of objects in departement layer: 96 Num points 47036, 1°) Using spatial functions: First run: 8953.493 ms, second run: 8808.133 ms 2°) Using Topology First run: 13743 ms, second run: 14028 ms (server stopped and restarted between runs. I'm interested in speeding up the case using topology. Note that ST_CreateTopoGeo needs to do all in memory which would be inappropriate for large datasets, the toTopoGeom function (and the underlying TopoGeo_addGeometry functions) would instead allow you to populate a topology incrementally. I will test this function and let you know. For sure I can tell you that topology isn't meant to be used as a temporary storage. I can imagine, for example, that you want multiple resolutions of your data. You could build the topology once and then generate multiple simplified versions of it, or even simplify on demand with an ST_Simplify(TopoGeometry) [ to be implemented ]. Nicolas ___ postgis-users mailing list postgis-users@postgis.refractions.net http://postgis.refractions.net/mailman/listinfo/postgis-users
Re: [postgis-users] Simplifying a multipolygon layer, keeping polygon connection
On Wed, Apr 11, 2012 at 11:14:48AM +0200, Nicolas Ribot wrote: Yes, here are some stats: Mac OS X Snow Leopard, SSD drive, 8 gb ram shared_buffers = 5120MB work_mem = 1024MB POSTGIS=2.0.0 r9605 GEOS=3.3.2-CAPI-1.7.2 PROJ=Rel. 4.7.1, 23 September 2009 GDAL=GDAL 1.8.0, released 2011/01/12 LIBXML=2.7.8 TOPOLOGY RASTER Number of objects in departement layer: 96 Num points 47036, 1°) Using spatial functions: First run: 8953.493 ms, second run: 8808.133 ms 2°) Using Topology First run: 13743 ms, second run: 14028 ms (server stopped and restarted between runs. Which steps involved by the two methods ? I can see some steps are in common: 1: Loading of input 2: Linking of faces to attribute (area_orig/area_intersection) Could the first step be taken out and the second be profiled separately ? --strk; ,--o-. | __/ |Delivering high quality PostGIS 2.0 ! | / 2.0 |http://strk.keybit.net - http://vizzuality.com `-o--' ___ postgis-users mailing list postgis-users@postgis.refractions.net http://postgis.refractions.net/mailman/listinfo/postgis-users
Re: [postgis-users] Invalid SRID for geography type: ee4
You can't load projected coordinates into geography, you'll need to load to geometry, then re-project to 4326 and cast to geography if you really want to use geography (why do you want to if your data is all neatly contained in Belgium and already projected?) P. On Tue, Apr 10, 2012 at 4:50 PM, Philippe watche...@gmail.com wrote: Hello, I'm running into an issue loading data from the Belgium IGN. I can transform a shapefile into Postgis loadable data with Geometry but not Geography. I can't imagine why. Any ideas ? Thanks The following works /usr/lib/postgresql/9.1/bin/shp2pgsql -WLATIN1 -s3812 -d -I /home/postgres/warehouse/shapes/BE/vector/Shape_L08/commune_polygon.shp be_ign /tmp/be_ign.sql The following fails /usr/lib/postgresql/9.1/bin/shp2pgsql -WLATIN1 -s3812 -G -d -I /home/postgres/warehouse/shapes/BE/vector/Shape_L08/commune_polygon.shp be_ign /tmp/be_ign.sql with Shapefile type: Polygon Postgis type: MULTIPOLYGON[2] Invalid SRID for geography type: ee4 Although 0xee4=3812 which is the EPSG for the Lambert 2008 projection used by that file. And my spatial ref table contains it: public=# select * from spatial_ref_sys where auth_srid=3812; srid | 3812 auth_name | EPSG auth_srid | 3812 srtext | PROJCS[ETRS89 / Belgian Lambert 2008,GEOGCS[ETRS89,DATUM[European_Terrestrial_Reference_System_1989,SPHEROID[GRS 1980,6378137,298.257222101,AUTHORITY[EPSG,7019]],AUTHORITY[EPSG,6258]],PRIMEM[Greenwich,0,AUTHORITY[EPSG,8901]],UNIT[degree,0.01745329251994328,AUTHORITY[EPSG,9122]],AUTHORITY[EPSG,4258]],UNIT[metre,1,AUTHORITY[EPSG,9001]],PROJECTION[Lambert_Conformal_Conic_2SP],PARAMETER[standard_parallel_1,49.84],PARAMETER[standard_parallel_2,51.16],PARAMETER[latitude_of_origin,50.797815],PARAMETER[central_meridian,4.3592158],PARAMETER[false_easting,649328],PARAMETER[false_northing,665262],AUTHORITY[EPSG,3812],AXIS[X,EAST],AXIS[Y,NORTH]] proj4text | +proj=lcc +lat_1=49.84 +lat_2=51.16 +lat_0=50.797815 +lon_0=4.3592158 +x_0=649328 +y_0=665262 +ellps=GRS80 +units=m +no_defs ___ 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] extraction of pixel from raster
hello can anybody plz help me in solving my problem -- I hv imported a raster in POSTGIS2.0 in a table. now i want to extract all pixels with only particular color like blue or some specific from that. how can i do that? Jyoti ___ postgis-users mailing list postgis-users@postgis.refractions.net http://postgis.refractions.net/mailman/listinfo/postgis-users
Re: [postgis-users] Simplifying a multipolygon layer, keeping polygon connection
On 11 April 2012 11:23, Sandro Santilli s...@keybit.net wrote: On Wed, Apr 11, 2012 at 11:14:48AM +0200, Nicolas Ribot wrote: Yes, here are some stats: Mac OS X Snow Leopard, SSD drive, 8 gb ram shared_buffers = 5120MB work_mem = 1024MB POSTGIS=2.0.0 r9605 GEOS=3.3.2-CAPI-1.7.2 PROJ=Rel. 4.7.1, 23 September 2009 GDAL=GDAL 1.8.0, released 2011/01/12 LIBXML=2.7.8 TOPOLOGY RASTER Number of objects in departement layer: 96 Num points 47036, 1°) Using spatial functions: First run: 8953.493 ms, second run: 8808.133 ms 2°) Using Topology First run: 13743 ms, second run: 14028 ms (server stopped and restarted between runs. Which steps involved by the two methods ? I can see some steps are in common: 1: Loading of input 2: Linking of faces to attribute (area_orig/area_intersection) Could the first step be taken out and the second be profiled separately ? Here are the tests I ran: a working table is created from the initial departement table, by dumping out the polygon: create table expl_dep as ( select gid, code_dept, (st_dump(geom)).* from departement ); This table is then used as input. ___ Method 1: using one SQL query: select (st_dump(st_polygonize(distinct geom))).geom as geom from ( select (st_dump(st_simplifyPreserveTopology(st_linemerge(st_union(geom)), 1))).geom as geom from ( select st_exteriorRing((st_dumpRings(geom)).geom) as geom from expl_dep ) as foo ) as bar; -- takes 6.1 to 6.4 seconds Face linking: select d.gid, d.code_dept, s.geom from departement d, simple_dep s where st_intersects(d.geom, s.geom) and st_area(st_intersection(d.geom, s.geom))/st_area(s.geom) 0.5; -- takes 1.2 to 1.3 s Method 2, using topology wrapped in a function: (topo1 and topo2 are created outside the function): create or replace function testSimplTopo() returns text as $$ select ST_CreateTopoGeo('topo1',ST_Collect(geom)) from expl_dep; select ST_CreateTopoGeo('topo2', geom) from ( select ST_Collect(st_simplifyPreserveTopology(geom, 1)) as geom from topo1.edge_data ) as foo; $$ language sql; -- takes 12.0 to 12.3 seconds. Face linking: with simple_face as ( select st_getFaceGeometry('topo2', face_id) as geom from topo2.face where face_id 0 ) select d.gid, d.code_dept, s.geom from simple_face s, departement d where st_intersects(d.geom, s.geom) and st_area(st_intersection(s.geom, d.geom))/st_area(s.geom) 0.5; -- takes 1.2 to 1.3, like the other query By the way, I tried to put calls to select dropTopology() into the functions I wrote, but it fails. Is it because it runs in the same transaction as the one that created the topo ? Nicolas ___ postgis-users mailing list postgis-users@postgis.refractions.net http://postgis.refractions.net/mailman/listinfo/postgis-users
Re: [postgis-users] How to design a database for continents, countries, regions, cities and POIs?
Hi, A good start is to check what has been already done in the domain of gazetteers (http://en.wikipedia.org/wiki/Gazetteer). You will probably find across the list some data and schemas (with hierachical structure) that should fit your needs. Marc-André Morin De : postgis-users-boun...@postgis.refractions.net [mailto:postgis-users-boun...@postgis.refractions.net] De la part de pcr...@pcreso.com Envoyé : 10 avril 2012 02:22 À : Michal Kubenka Cc : postgis-users@postgis.refractions.net Objet : Re: [postgis-users] How to design a database for continents,countries, regions, cities and POIs? Hi Michal, One suggestion... There are two ways (at least :-) to do this in a RDBMS. You can have the spatial relationship implicit in the feature geometries, so a spatial query is used, for example, to determine the cities within a country: select * from polygons a, polygons b where a.type = 'city' and b.type='country' and b.name='Italy'; While flexible effective, relying on spatial queries for quick searches with polygons with many thousands, or even millions of records may not be ideal. The other approach is to explicitly predefine these relationships, so a column for each polygon feature stores the parent id. Simplistically assuming the parent of a city is the country containing it, rather than navigating the hierarchy, the above query becomes: select * from polygons where type=city and parent_id = (select id from polygons where type = 'country' and name = 'Italy'); Even with both structures optimised indexed, the latter is likely to be much faster. No join is required. Given the country containing a city is a pretty static relationship, I suggest predefining to optimise query performance makes sense. If you store the heirarchies as predefined levels then a heirarchical search using the recursive with capability- see: http://www.postgresql.org/docs/8.4/static/queries-with.html is perhaps possible, to invoke searches up down the tree. So use Postgis to determine the parent id using a spatial function, then store this as an indexed id. HTH, Brent Wood I'd say there are several approaches you could take to build a viable database, the optimal one is defined by your use case: the sorts of queries you want to apply. --- On Tue, 4/10/12, Michal Kubenka mkube...@gmail.com wrote: From: Michal Kubenka mkube...@gmail.com Subject: Re: [postgis-users] How to design a database for continents, countries, regions, cities and POIs? To: pcr...@pcreso.com Cc: PostGIS Users Discussion postgis-users@postgis.refractions.net Date: Tuesday, April 10, 2012, 9:59 AM Actually what we need is some hierarchical base for relationship between countries, cities, regions, etc. Main goal of the application will be collecting data from many sources about specific cities, regions, countries and so on, and store it in database. Let's say we will have city Rome, we collect some info about this city into database from couple sources. And we need to know that Rome is in province Rome, sub-region Lazio in region Lazio, country Italy. So system should be flexible to allow create such relation from real world. That's why I would choose two tables: 1) `polygons` - which can store countries, regions, sub-regions, provinces etc. 2) `points` - which can store cities and POIs Thanks. Michal K. On Mon, Apr 9, 2012 at 8:11 PM, pcr...@pcreso.com wrote: Are you planning to store multiple versions of these polygons, for zoom layers? Generally you need a high res version (eg: coastline) when zoomed in (large scale) and a lower resolution version when zoomed out (you can't see don't need the detail. This may or may not have an impact on your eventual data model, but it is worth ensuring you take this into account during the data modeling process. You can have a model where each feature has multiple geometry columns associated with it in the one table, or an approach which has the geometries in separate tables, using ID's to link to the aspatial attributes. The former is a simpler, monolithic solution, the latter is more complex but allows more use of tablespaces underlying Postgres optimisation. You may also find you need to carry out joins (identify relationships between types of polygon, eg: cities within counties within states within countries, and this may perform better with a denormalised structure with separate tables for different categories of polygon. One example you might look at is the OSM data model. Not quite what you are describing, but a robust well tested model for global roads related spatial data, which does not use Postgis at all. http://booki.flossmanuals.net/openstreetmap/_draft/_v/1.0/the-osm-data-model/ --- On Mon, 4/9/12, mkubenka
Re: [postgis-users] Simplifying a multipolygon layer, keeping polygon connection
On Wed, Apr 11, 2012 at 01:40:01PM +0200, Nicolas Ribot wrote: Method 1: using one SQL query: select (st_dump(st_polygonize(distinct geom))).geom as geom from ( select (st_dump(st_simplifyPreserveTopology(st_linemerge(st_union(geom)), 1))).geom as geom from ( select st_exteriorRing((st_dumpRings(geom)).geom) as geom from expl_dep ) as foo ) as bar; -- takes 6.1 to 6.4 seconds ... Method 2, using topology wrapped in a function: (topo1 and topo2 are created outside the function): create or replace function testSimplTopo() returns text as $$ select ST_CreateTopoGeo('topo1',ST_Collect(geom)) from expl_dep; select ST_CreateTopoGeo('topo2', geom) from ( select ST_Collect(st_simplifyPreserveTopology(geom, 1)) as geom from topo1.edge_data ) as foo; $$ language sql; -- takes 12.0 to 12.3 seconds. Why are you populating 'topo1' topology at all ? You're doing the hard work twice, which matches the double time it takes compared with the direct call. --strk; ,--o-. | __/ |Delivering high quality PostGIS 2.0 ! | / 2.0 |http://strk.keybit.net - http://vizzuality.com `-o--' ___ postgis-users mailing list postgis-users@postgis.refractions.net http://postgis.refractions.net/mailman/listinfo/postgis-users
Re: [postgis-users] Simplifying a multipolygon layer, keeping polygon connection
On Wed, Apr 11, 2012 at 02:21:15PM +0200, Sandro Santilli wrote: On Wed, Apr 11, 2012 at 01:40:01PM +0200, Nicolas Ribot wrote: Method 1: using one SQL query: select (st_dump(st_polygonize(distinct geom))).geom as geom from ( select (st_dump(st_simplifyPreserveTopology(st_linemerge(st_union(geom)), 1))).geom as geom from ( select st_exteriorRing((st_dumpRings(geom)).geom) as geom from expl_dep ) as foo ) as bar; -- takes 6.1 to 6.4 seconds ... Method 2, using topology wrapped in a function: (topo1 and topo2 are created outside the function): create or replace function testSimplTopo() returns text as $$ select ST_CreateTopoGeo('topo1',ST_Collect(geom)) from expl_dep; select ST_CreateTopoGeo('topo2', geom) from ( select ST_Collect(st_simplifyPreserveTopology(geom, 1)) as geom from topo1.edge_data ) as foo; $$ language sql; -- takes 12.0 to 12.3 seconds. Why are you populating 'topo1' topology at all ? You're doing the hard work twice, which matches the double time it takes compared with the direct call. Ok, I see you're using 'topo1' as input for the simplification. Instead you may want to simplify the input directly: select ST_CreateTopoGeo('topo1', ST_Collect( ST_SimplifyPreserveTopology(geom, 1) )) from expl_dep; --strk; ,--o-. | __/ |Delivering high quality PostGIS 2.0 ! | / 2.0 |http://strk.keybit.net - http://vizzuality.com `-o--' ___ postgis-users mailing list postgis-users@postgis.refractions.net http://postgis.refractions.net/mailman/listinfo/postgis-users
Re: [postgis-users] How to design a database for continents, countries, regions, cities and POIs?
On Mon, Apr 09, 2012 at 11:59:14PM +0200, Michal Kubenka wrote: Actually what we need is some hierarchical base for relationship between countries, cities, regions, etc. PostGIS topology supports hirearchical modeling of layers, with each object in upper layer defined by items of the lower layer, down to primitives. --strk; ,--o-. | __/ |Delivering high quality PostGIS 2.0 ! | / 2.0 |http://strk.keybit.net - http://vizzuality.com `-o--' ___ postgis-users mailing list postgis-users@postgis.refractions.net http://postgis.refractions.net/mailman/listinfo/postgis-users
Re: [postgis-users] Simplifying a multipolygon layer, keeping polygon connection
On Wed, Apr 11, 2012 at 02:30:31PM +0200, Nicolas Ribot wrote: Ok, I see you're using 'topo1' as input for the simplification. Instead you may want to simplify the input directly: select ST_CreateTopoGeo('topo1', ST_Collect( ST_SimplifyPreserveTopology(geom, 1) )) from expl_dep; Yes, this was the reason: hoping the noded edges could be simplified without loosing the connectivity (which is the case :) ) In one pass, I fall back to the classical simplify with broken topology as shown on the picture. But that's due to the way you associate faces with original polygons. If each and every face is associated to exactly _one_ polygon you'd be fine. You could get PointOnSurface of each face and check each of the original polygons for containment (ST_Contains). Alternatively you could still construct the topology with full-precision edges and then simplify the edges in-place (possibly after a CopyTopology), using ST_ChangeEdgeGeom. That way the _original_ nodes will remain the same and you'd be only simplifying the edges connecting them. Beware in this case that simplification might still introduce topology errors but such errors would be cought by ST_ChangeEdgeGeom with throwing an exception. You may catch the exception to reduce tolerance for a given edge or bypass all checks by editing the edge_data table directly and cross your fingers :) ValidateTopology should then tell you what you broke, if you want to fix manually after the fact. --strk; ,--o-. | __/ |Delivering high quality PostGIS 2.0 ! | / 2.0 |http://strk.keybit.net - http://vizzuality.com `-o--' ___ postgis-users mailing list postgis-users@postgis.refractions.net http://postgis.refractions.net/mailman/listinfo/postgis-users
Re: [postgis-users] extraction of pixel from raster
You will have to work with ST_Reclass(). See the doc: http://postgis.refractions.net/documentation/manual-svn/RT_ST_Reclass.html and mostly the third example. Pierre -Original Message- From: postgis-users-boun...@postgis.refractions.net [mailto:postgis-users- boun...@postgis.refractions.net] On Behalf Of jyoti gajrani Sent: Wednesday, April 11, 2012 7:36 AM To: postgis-users@postgis.refractions.net Subject: [postgis-users] extraction of pixel from raster hello can anybody plz help me in solving my problem -- I hv imported a raster in POSTGIS2.0 in a table. now i want to extract all pixels with only particular color like blue or some specific from that. how can i do that? Jyoti ___ postgis-users mailing list postgis-users@postgis.refractions.net http://postgis.refractions.net/mailman/listinfo/postgis-users
Re: [postgis-users] NDVI Calculation from two bands within one Raster
Sorry to keep reviving this thread, espcially with a basic question. I am using: SELECT ST_MapAlgebraExpr(a.rast, 4, b.rast, 3, '([rast1] - [rast2]) /([rast1] + [rast2])::float') WHEN (rast1 + rast2 = 0) THEN = 999 ELSE ([rast1] - [rast2]) /([rast1] + [rast2])::float END' FROM ndvi a, ndvi b; Where the WHEN statement is being used to prevent the error division by zero. However this returns: syntax error at or near WHEN Can anyone advise on where I'm structuring this query wrongly. Regards, James - GIS Undergraduate -- View this message in context: http://postgis.17.n6.nabble.com/NDVI-Calculation-from-two-bands-within-one-Raster-tp4656995p4832397.html Sent from the PostGIS - User mailing list archive at Nabble.com. ___ postgis-users mailing list postgis-users@postgis.refractions.net http://postgis.refractions.net/mailman/listinfo/postgis-users
Re: [postgis-users] How to design a database for continents, countries, regions, cities and POIs?
Thank you. On Wed, Apr 11, 2012 at 2:29 PM, Sandro Santilli s...@keybit.net wrote: On Mon, Apr 09, 2012 at 11:59:14PM +0200, Michal Kubenka wrote: Actually what we need is some hierarchical base for relationship between countries, cities, regions, etc. PostGIS topology supports hirearchical modeling of layers, with each object in upper layer defined by items of the lower layer, down to primitives. --strk; ,--o-. | __/ |Delivering high quality PostGIS 2.0 ! | / 2.0 |http://strk.keybit.net - http://vizzuality.com `-o--' ___ 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] How to design a database for continents, countries, regions, cities and POIs?
Question out ignorance because I was not aware that this capability existed. I might have a geometry for a country and separate geometries for the next layer of administrative regions (e.g. states in the United States). But there are often minor differences in the polygons such that the polygon for the country is nearly but not exactly the union of the states. How is this handled? Thanks. -=- Jerry On Apr 11, 2012, at 11:42 AM, Michal Kubenka wrote: Thank you. On Wed, Apr 11, 2012 at 2:29 PM, Sandro Santilli s...@keybit.net wrote: On Mon, Apr 09, 2012 at 11:59:14PM +0200, Michal Kubenka wrote: Actually what we need is some hierarchical base for relationship between countries, cities, regions, etc. PostGIS topology supports hirearchical modeling of layers, with each object in upper layer defined by items of the lower layer, down to primitives. --strk; ,--o-. | __/ |Delivering high quality PostGIS 2.0 ! | / 2.0 |http://strk.keybit.net - http://vizzuality.com `-o--' ___ 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] NDVI Calculation from two bands within one Raster
Your CASE syntax is wrong: http://www.postgresql.org/docs/9.1/interactive/functions-conditional.html#FUNCTIONS-CASE Should look like: ST_MapAlgebraExpr(a.rast, 4, b.rast, 3, 'CASE WHEN ([rast1] + [rast2]) = 0 THEN 999 ELSE ([rast1] - [rast2]) /([rast1] + [rast2])::float END') FROM ndvi a, ndvi b; Pierre -Original Message- From: postgis-users-boun...@postgis.refractions.net [mailto:postgis-users- boun...@postgis.refractions.net] On Behalf Of JamesH Sent: Wednesday, April 11, 2012 10:56 AM To: postgis-users@postgis.refractions.net Subject: Re: [postgis-users] NDVI Calculation from two bands within one Raster Sorry to keep reviving this thread, espcially with a basic question. I am using: SELECT ST_MapAlgebraExpr(a.rast, 4, b.rast, 3, '([rast1] - [rast2]) /([rast1] + [rast2])::float') WHEN (rast1 + rast2 = 0) THEN = 999 ELSE ([rast1] - [rast2]) /([rast1] + [rast2])::float END' FROM ndvi a, ndvi b; Where the WHEN statement is being used to prevent the error division by zero. However this returns: syntax error at or near WHEN Can anyone advise on where I'm structuring this query wrongly. Regards, James - GIS Undergraduate -- View this message in context: http://postgis.17.n6.nabble.com/NDVI- Calculation-from-two-bands-within-one-Raster-tp4656995p4832397.html Sent from the PostGIS - User mailing list archive at Nabble.com. ___ 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] How to design a database for continents, countries, regions, cities and POIs?
On Wed, Apr 11, 2012 at 11:47:02AM -0400, Jerry Carter wrote: I might have a geometry for a country and separate geometries for the next layer of administrative regions (e.g. states in the United States). But there are often minor differences in the polygons such that the polygon for the country is nearly but not exactly the union of the states. How is this handled? In the topology model your higher layer will be defined by composition of lower layer items. So you would only say that United States is formed by all the countries. And for each country you would only list the counties they are formed of. And for each county the parcels. And for each parcel the primitive faces. And each face is defined by its edges. The only actual geometries involved in all of the above would be the edge geometries. Not any other geometry. So there's no difference because primitive elements are singletons. --strk; ,--o-. | __/ |Delivering high quality PostGIS 2.0 ! | / 2.0 |http://strk.keybit.net - http://vizzuality.com `-o--' ___ postgis-users mailing list postgis-users@postgis.refractions.net http://postgis.refractions.net/mailman/listinfo/postgis-users
Re: [postgis-users] How to design a database for continents, countries, regions, cities and POIs?
Thanks for the quick answer. Just read your writeup on the subject [1]. This looks like a _very nice_ addition and should be extremely helpful. I've got some experimenting to do! -=- Jerry [1] http://strk.keybit.net/projects/postgis/Paris2011_TopologyWithPostGIS_2_0.pdf In otherwords, you skip On Apr 11, 2012, at 11:51 AM, Sandro Santilli wrote: On Wed, Apr 11, 2012 at 11:47:02AM -0400, Jerry Carter wrote: I might have a geometry for a country and separate geometries for the next layer of administrative regions (e.g. states in the United States). But there are often minor differences in the polygons such that the polygon for the country is nearly but not exactly the union of the states. How is this handled? In the topology model your higher layer will be defined by composition of lower layer items. So you would only say that United States is formed by all the countries. And for each country you would only list the counties they are formed of. And for each county the parcels. And for each parcel the primitive faces. And each face is defined by its edges. The only actual geometries involved in all of the above would be the edge geometries. Not any other geometry. So there's no difference because primitive elements are singletons. --strk; ,--o-. | __/ |Delivering high quality PostGIS 2.0 ! | / 2.0 |http://strk.keybit.net - http://vizzuality.com `-o--' ___ 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] NDVI Calculation from two bands within one Raster
Thanks, that has worked and can visualise through OpenJump. The Raster it produces (below), has only values of zero or one in its attributes. I'm assuming in the calculation its not keeping the values as Floats but its converting it to integer? Must be rounding them to either zero or one. Is there a way to solve this i.e. keep the pixels as floating integers? http://postgis.17.n6.nabble.com/file/n4832776/ndvi.jpg Regards, James - GIS Undergraduate -- View this message in context: http://postgis.17.n6.nabble.com/NDVI-Calculation-from-two-bands-within-one-Raster-tp4656995p4832776.html Sent from the PostGIS - User mailing list archive at Nabble.com. ___ postgis-users mailing list postgis-users@postgis.refractions.net http://postgis.refractions.net/mailman/listinfo/postgis-users
Re: [postgis-users] NDVI Calculation from two bands within one Raster
Thanks, that has worked and can visualise through OpenJump. The Raster it produces (below), has only values of zero or one in its attributes. I'm assuming in the calculation its not keeping the values as Floats but its converting it to integer? Must be rounding them to either zero or one. Is there a way to solve this i.e. keep the pixels as floating integers? http://postgis.17.n6.nabble.com/file/n4832782/ndvi.jpg Kind Regards, James - GIS Undergraduate -- View this message in context: http://postgis.17.n6.nabble.com/NDVI-Calculation-from-two-bands-within-one-Raster-tp4656995p4832782.html Sent from the PostGIS - User mailing list archive at Nabble.com. ___ postgis-users mailing list postgis-users@postgis.refractions.net http://postgis.refractions.net/mailman/listinfo/postgis-users
Re: [postgis-users] NDVI Calculation from two bands within one Raster
Are you using ST_dumpAsPolygons() to display in OpenJump? Depending on the GDAL version you're using, ST_DumpAsPolygons() might be converting floats to integers. Try looking at individual pixel values with ST_Value() or us ST_PixelAsPolygons() Pierre -Original Message- From: postgis-users-boun...@postgis.refractions.net [mailto:postgis-users- boun...@postgis.refractions.net] On Behalf Of JamesH Sent: Wednesday, April 11, 2012 12:39 PM To: postgis-users@postgis.refractions.net Subject: Re: [postgis-users] NDVI Calculation from two bands within one Raster Thanks, that has worked and can visualise through OpenJump. The Raster it produces (below), has only values of zero or one in its attributes. I'm assuming in the calculation its not keeping the values as Floats but its converting it to integer? Must be rounding them to either zero or one. Is there a way to solve this i.e. keep the pixels as floating integers? http://postgis.17.n6.nabble.com/file/n4832782/ndvi.jpg Kind Regards, James - GIS Undergraduate -- View this message in context: http://postgis.17.n6.nabble.com/NDVI- Calculation-from-two-bands-within-one-Raster-tp4656995p4832782.html Sent from the PostGIS - User mailing list archive at Nabble.com. ___ 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] NDVI Calculation from two bands within one Raster
Are you using ST_dumpAsPolygons() to display in OpenJump? Try looking at individual pixel values with ST_Value() or us ST_PixelAsPolygons() I was using DumpAsPolygons. Used ST_PixelAsPolygons and this returned the same, either zero, one or NullData. If it helps my GDAL version is: GDAL=GDAL 1.9.0, released 2011/12/29 Try looking at individual pixel values with ST_Value() Examined and in the final output, St_ValueCount returns: (0, 32980) (1, 519) Looked closer at the calculation with ST_ValueCount and [rast1] - [rast2] = (0, 17458) [rast1]+[rast2] = (0, 16223) so at some point these must be dividing by zero? Could this be affecting the NDVI? Kind Regards, James - GIS Undergraduate -- View this message in context: http://postgis.17.n6.nabble.com/NDVI-Calculation-from-two-bands-within-one-Raster-tp4656995p4832963.html Sent from the PostGIS - User mailing list archive at Nabble.com. ___ postgis-users mailing list postgis-users@postgis.refractions.net http://postgis.refractions.net/mailman/listinfo/postgis-users
Re: [postgis-users] NDVI Calculation from two bands within one Raster
Add a pixeltype parameter to your ST_MapAlgebraExpr call. Pierre -Original Message- From: postgis-users-boun...@postgis.refractions.net [mailto:postgis-users- boun...@postgis.refractions.net] On Behalf Of JamesH Sent: Wednesday, April 11, 2012 1:44 PM To: postgis-users@postgis.refractions.net Subject: Re: [postgis-users] NDVI Calculation from two bands within one Raster Are you using ST_dumpAsPolygons() to display in OpenJump? Try looking at individual pixel values with ST_Value() or us ST_PixelAsPolygons() I was using DumpAsPolygons. Used ST_PixelAsPolygons and this returned the same, either zero, one or NullData. If it helps my GDAL version is: GDAL=GDAL 1.9.0, released 2011/12/29 Try looking at individual pixel values with ST_Value() Examined and in the final output, St_ValueCount returns: (0, 32980) (1, 519) Looked closer at the calculation with ST_ValueCount and [rast1] - [rast2] = (0, 17458) [rast1]+[rast2] = (0, 16223) so at some point these must be dividing by zero? Could this be affecting the NDVI? Kind Regards, James - GIS Undergraduate -- View this message in context: http://postgis.17.n6.nabble.com/NDVI- Calculation-from-two-bands-within-one-Raster-tp4656995p4832963.html Sent from the PostGIS - User mailing list archive at Nabble.com. ___ 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] Method to remove overlaps in a layer
So I have been creating a very simple example based on the example in the toTopoGeom() documentation page: http://www.postgis.org/documentation/manual-svn/toTopoGeom.html SELECT CreateTopology('topo1',31467, 10); CREATE TABLE ovlp.test02_topo (id integer); SELECT AddTopoGeometryColumn('topo1', 'ovlp', 'test02_topo', 'topo', 'POLYGON'); INSERT INTO ovlp.test02_topo (id, topo) SELECT id, topology.toTopoGeom(geom, 'topo1', 1) FROM ovlp.test02 And I get: ** Error ** ERROR: interpolate_point4d: invalid F (1) SQL state: XX000 Context: PL/pgSQL function st_modedgesplit line 94 at assignment PL/pgSQL function topogeo_addpoint line 91 at assignment PL/pgSQL function topogeo_addlinestring line 132 at assignment SQL statement SELECT array_cat(edges, array_agg(x)) FROM ( select topology.TopoGeo_addLinestring(atopology, rec.geom, tol) as x ) as foo PL/pgSQL function topogeo_addpolygon line 27 at assignment SQL statement INSERT INTO topo1.relation(topogeo_id, layer_id, element_type, element_id) SELECT 5, 1, 3, topogeo_addPolygon('topo1', '010320EB7A01000880FD384B41F4615441D84B4B41004050695441B55D4B415C635441BD594B41965B5441C8694B41625554410080E15C4B410080474D54410080ED404B416C5054410080FD384B41F4615441'::geometry, 10); PL/pgSQL function totopogeom line 129 at EXECUTE statement I join a picture of the original layer to be converted to a topologic layer. I'm still at beta 5... Should upgrading solve that particular problem? Pierre So what would be the normal/easiest steps to convert a messy polygon coverage into a clean topology? Is it documented somewhere? My guess: 1- SET search_path TO topology,public; You shouldn't need this, when you load topology.sql you should get topology already appended to the end of the search_path associated with the database. 2- SELECT CreateTopoGeom('test') Yep. 3- SELECT toTopoGeom(geom, 'test', 1) FROM mymessyone; You didn't create a layer, see AddTopoGeometryColumn. The third argument is a layer id, as returned by that function. or SELECT ST_CreateTopoGeo('test', geom) FROM mymessyone; This one only works starting with an empty topology so you'll need to pass it a full collection: SELECT ST_CreateTopoGeo('test', ST_Collect(geom)) FROM mymessyone; But I'd recommend using toTopoGeom instead, to keep the linking between attributes and geometries. What happens when a polygon to be added to the topology overlaps a polygon already in the topology? Two overlapping rectangles would produce a total of 3 faces. If you're returning TopoGeometry objects (your step 2) both of them will be defined by 2 faces, having one face in common. attachment: overlaps.jpeg___ postgis-users mailing list postgis-users@postgis.refractions.net http://postgis.refractions.net/mailman/listinfo/postgis-users
Re: [postgis-users] Method to remove overlaps in a layer
An even simpler example leads to another type of error: SELECT CreateTopology('topo3',0, 10); DROP TABLE IF EXISTS ovlp.test03; CREATE TABLE ovlp.test03 AS SELECT 1 id, ST_GeomFromText('POLYGON((0 1, 2 2, 2 0, 0 1))') geom UNION ALL SELECT 2 id, ST_GeomFromText('POLYGON((1 1, 3 2, 3 0, 1 1))') geom CREATE TABLE ovlp.test03_topo (id integer); SELECT AddTopoGeometryColumn('topo3', 'ovlp', 'test03_topo', 'topo', 'POLYGON'); -- 1 INSERT INTO ovlp.test03_topo (id, topo) SELECT id, topology.toTopoGeom(geom, 'topo3', 1) FROM ovlp.test03 gives: ** Error ** ERROR: SQL/MM Spatial exception - edge crosses node. SQL state: P0001 Context: PL/pgSQL function topogeo_addpoint line 96 at assignment PL/pgSQL function topogeo_addlinestring line 125 at assignment SQL statement SELECT array_cat(edges, array_agg(x)) FROM ( select topology.TopoGeo_addLinestring(atopology, rec.geom, tol) as x ) as foo PL/pgSQL function topogeo_addpolygon line 27 at assignment SQL statement INSERT INTO topo3.relation(topogeo_id, layer_id, element_type, element_id) SELECT 2, 1, 3, topogeo_addPolygon('topo3', '01030001000400F03FF03F084000400840F03FF03F'::geometry, 10); PL/pgSQL function totopogeom line 129 at EXECUTE statement Pierre -Original Message- From: postgis-users-boun...@postgis.refractions.net [mailto:postgis-users- boun...@postgis.refractions.net] On Behalf Of Pierre Racine Sent: Wednesday, April 11, 2012 2:40 PM To: PostGIS Users Discussion Subject: Re: [postgis-users] Method to remove overlaps in a layer So I have been creating a very simple example based on the example in the toTopoGeom() documentation page: http://www.postgis.org/documentation/manual-svn/toTopoGeom.html SELECT CreateTopology('topo1',31467, 10); CREATE TABLE ovlp.test02_topo (id integer); SELECT AddTopoGeometryColumn('topo1', 'ovlp', 'test02_topo', 'topo', 'POLYGON'); INSERT INTO ovlp.test02_topo (id, topo) SELECT id, topology.toTopoGeom(geom, 'topo1', 1) FROM ovlp.test02 And I get: ** Error ** ERROR: interpolate_point4d: invalid F (1) SQL state: XX000 Context: PL/pgSQL function st_modedgesplit line 94 at assignment PL/pgSQL function topogeo_addpoint line 91 at assignment PL/pgSQL function topogeo_addlinestring line 132 at assignment SQL statement SELECT array_cat(edges, array_agg(x)) FROM ( select topology.TopoGeo_addLinestring(atopology, rec.geom, tol) as x ) as foo PL/pgSQL function topogeo_addpolygon line 27 at assignment SQL statement INSERT INTO topo1.relation(topogeo_id, layer_id, element_type, element_id) SELECT 5, 1, 3, topogeo_addPolygon('topo1', '010320EB7A01000880FD384B41F46154 41D84B4B41004050695441B55D4B415C6354 41BD594B41965B5441C8694B41625554 410080E15C4B410080474D54410080ED404B416C5054 410080FD384B41F4615441'::geometry, 10); PL/pgSQL function totopogeom line 129 at EXECUTE statement I join a picture of the original layer to be converted to a topologic layer. I'm still at beta 5... Should upgrading solve that particular problem? Pierre So what would be the normal/easiest steps to convert a messy polygon coverage into a clean topology? Is it documented somewhere? My guess: 1- SET search_path TO topology,public; You shouldn't need this, when you load topology.sql you should get topology already appended to the end of the search_path associated with the database. 2- SELECT CreateTopoGeom('test') Yep. 3- SELECT toTopoGeom(geom, 'test', 1) FROM mymessyone; You didn't create a layer, see AddTopoGeometryColumn. The third argument is a layer id, as returned by that function. or SELECT ST_CreateTopoGeo('test', geom) FROM mymessyone; This one only works starting with an empty topology so you'll need to pass it a full collection: SELECT ST_CreateTopoGeo('test', ST_Collect(geom)) FROM mymessyone; But I'd recommend using toTopoGeom instead, to keep the linking between attributes and geometries. What happens when a polygon to be added to the topology overlaps a polygon already in the topology? Two overlapping rectangles would produce a total of 3 faces. If you're returning TopoGeometry objects (your step 2) both of them will be defined by 2 faces, having one face in common. ___ postgis-users mailing list postgis-users@postgis.refractions.net http://postgis.refractions.net/mailman/listinfo/postgis-users
Re: [postgis-users] NDVI Calculation from two bands within one Raster
That has worked brilliantly and is exactly what I wanted. I added '32BF' to the the end. Thank you Pierre James - GIS Undergraduate -- View this message in context: http://postgis.17.n6.nabble.com/NDVI-Calculation-from-two-bands-within-one-Raster-tp4656995p4833176.html Sent from the PostGIS - User mailing list archive at Nabble.com. ___ postgis-users mailing list postgis-users@postgis.refractions.net http://postgis.refractions.net/mailman/listinfo/postgis-users
Re: [postgis-users] Method to remove overlaps in a layer
On Wed, Apr 11, 2012 at 02:40:22PM -0400, Pierre Racine wrote: ERROR: interpolate_point4d: invalid F (1) It's the very first time I see such error. Invalid F ? What's that ? I'm still at beta 5... Should upgrading solve that particular problem? If you send the data I'll be able to tell. Note that trunk (to be 2.0.1) already has a couple more fixes to topology. --strk; ,--o-. | __/ |Delivering high quality PostGIS 2.0 ! | / 2.0 |http://strk.keybit.net - http://vizzuality.com `-o--' ___ postgis-users mailing list postgis-users@postgis.refractions.net http://postgis.refractions.net/mailman/listinfo/postgis-users
Re: [postgis-users] Method to remove overlaps in a layer
On Wed, Apr 11, 2012 at 02:53:07PM -0400, Pierre Racine wrote: An even simpler example leads to another type of error: SELECT CreateTopology('topo3',0, 10); DROP TABLE IF EXISTS ovlp.test03; CREATE TABLE ovlp.test03 AS SELECT 1 id, ST_GeomFromText('POLYGON((0 1, 2 2, 2 0, 0 1))') geom UNION ALL SELECT 2 id, ST_GeomFromText('POLYGON((1 1, 3 2, 3 0, 1 1))') geom CREATE TABLE ovlp.test03_topo (id integer); SELECT AddTopoGeometryColumn('topo3', 'ovlp', 'test03_topo', 'topo', 'POLYGON'); -- 1 INSERT INTO ovlp.test03_topo (id, topo) SELECT id, topology.toTopoGeom(geom, 'topo3', 1) FROM ovlp.test03 gives: ** Error ** ERROR: SQL/MM Spatial exception - edge crosses node. Ouch, you're lucky. It kills the backend here !! I'm on it. --strk; ,--o-. | __/ |Delivering high quality PostGIS 2.0 ! | / 2.0 |http://strk.keybit.net - http://vizzuality.com `-o--' ___ postgis-users mailing list postgis-users@postgis.refractions.net http://postgis.refractions.net/mailman/listinfo/postgis-users
[postgis-users] linked point and polygon tables
Hello, Suppose we have in PostGIS table with point geometry, ID column and some other attributes. Is it possible to create table with polygon/polyline geometry which is linked to table with point geometry ? I mean such a link where: - geometry of polygon/polyline is not created by set of coordinates but is created by (refering to) set of ID's from point table, - when some point object (from point geometry table) is removed it should automatically influence the shape of polygon/polyline. Regards, Piotr ___ postgis-users mailing list postgis-users@postgis.refractions.net http://postgis.refractions.net/mailman/listinfo/postgis-users