[postgis-users] ST_Boundary(geometry) for polyon to line conversion?

2012-04-11 Thread Robert Buckley
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

2012-04-11 Thread Nicolas Ribot
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?

2012-04-11 Thread Nicolas Ribot
 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

2012-04-11 Thread Nicolas Ribot
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

2012-04-11 Thread Sandro Santilli
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

2012-04-11 Thread Sandro Santilli
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

2012-04-11 Thread Jose Carlos Martinez

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

2012-04-11 Thread Sandro Santilli
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

2012-04-11 Thread Sandro Santilli
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

2012-04-11 Thread Nicolas Ribot
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

2012-04-11 Thread Sandro Santilli
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

2012-04-11 Thread Paul Ramsey
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

2012-04-11 Thread jyoti gajrani
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

2012-04-11 Thread Nicolas Ribot
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?

2012-04-11 Thread Morin , Marc-André
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

2012-04-11 Thread Sandro Santilli
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

2012-04-11 Thread Sandro Santilli
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?

2012-04-11 Thread Sandro Santilli
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

2012-04-11 Thread Sandro Santilli
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

2012-04-11 Thread Pierre Racine
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

2012-04-11 Thread JamesH
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?

2012-04-11 Thread Michal Kubenka
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?

2012-04-11 Thread Jerry Carter
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

2012-04-11 Thread Pierre Racine
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?

2012-04-11 Thread Sandro Santilli
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?

2012-04-11 Thread Jerry Carter
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

2012-04-11 Thread JamesH
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

2012-04-11 Thread JamesH
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

2012-04-11 Thread Pierre Racine
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

2012-04-11 Thread JamesH
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

2012-04-11 Thread Pierre Racine
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

2012-04-11 Thread Pierre Racine
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

2012-04-11 Thread Pierre Racine
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

2012-04-11 Thread JamesH
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

2012-04-11 Thread Sandro Santilli
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

2012-04-11 Thread Sandro Santilli
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

2012-04-11 Thread Piotr Pachół

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