Re: [postgis-users] Snapping linestrings with given tolerance
Hi Nicolas, Thanks a lot for the tips. I will have a look at the thread and method you suggested and keep you posted. G ___ postgis-users mailing list postgis-users@lists.osgeo.org http://lists.osgeo.org/cgi-bin/mailman/listinfo/postgis-users
[postgis-users] Subject: Re: raster2pgsql 2.1.1 on Windows 7 64bit - ERREUR: erreur de syntaxe
Regina, Remi: No space in filenames. GDAL version is 1.11dev. I don't have any problems importing shapefiles into PostGIS. Encoding was already UTF8 but I tried your suggestion (SET PGCLIENT...) and I still have errors with SQL files having incomplete last line. I noticed that my LC_COLLATE and LC_TYPE are both set 'French.Canadian.1252'. Maybe this is the source of the problem. I tried creating a new cluster with --no-locale to make some tests with raster2pgsql but I screwed up and messed with my current cluster in the process... Will do more tests when I'm back from Easter Holiday and I fixed my cluster problem. Thanks for your help. Best, Guillaume ___ postgis-users mailing list postgis-users@lists.osgeo.org http://lists.osgeo.org/cgi-bin/mailman/listinfo/postgis-users
[postgis-users] Subject: Re: raster2pgsql 2.1.1 on Windows 7 64bit - ERREUR: erreur de syntaxe
I tried the same command as in my previous post on my PC at home. It's Windows 7 64bit as well and the cluster has the same locale ('French.Canadian.1252'.) and encoding (UTF8). The command was successful so my work colleague and I, who experiences the same behaviour as the one described previously on this list, suspect that it may be due to some problems or conflict between PostGIS and the hardware (we have similar systems): is this a sensible hypothesis? Thanks for your help, Guillaume > Message: 4 > Date: Thu, 17 Apr 2014 17:03:58 -0400 > From: Guillaume Drolet > To: "postgis-users@lists.osgeo.org" > Subject: [postgis-users] Subject: Re: raster2pgsql 2.1.1 on Windows 7 > 64bit - ERREUR: erreur de syntaxe > Message-ID: > > Content-Type: text/plain; charset=UTF-8 > > Regina, Remi: > > No space in filenames. GDAL version is 1.11dev. I don't have any > problems importing shapefiles into PostGIS. > > Encoding was already UTF8 but I tried your suggestion (SET > PGCLIENT...) and I still have errors with SQL files having incomplete > last line. > > I noticed that my LC_COLLATE and LC_TYPE are both set > 'French.Canadian.1252'. Maybe this is the source of the problem. I > tried creating a new cluster with --no-locale to make some tests with > raster2pgsql but I screwed up and messed with my current cluster in > the process... > > Will do more tests when I'm back from Easter Holiday and I fixed my > cluster problem. > Thanks for your help. > > Best, > > Guillaume ___ postgis-users mailing list postgis-users@lists.osgeo.org http://lists.osgeo.org/cgi-bin/mailman/listinfo/postgis-users
[postgis-users] Subject: Re: raster2pgsql 2.1.1 on Windows 7 64bit - ERREUR: erreur de syntaxe
Hi Mateusz, Thanks for helping. What you suggested, i.e. outputting raster2pgsql to an SQL file instead of piping it to psql, I have done already. I also posted examples of where the lines are broken in previous post (see also Jean-Daniel Sylvain's older post on the same topic). It doesn't seem to break on non-ASCII characters. For example, this one breaks in the middle of the file name at line 168: ... INSERT INTO "mod09a1.a249.h12v04.005.2006268184041_sur_refl_b01" ("rast","filename") VALUES ('010100A544B7031A11713FA544B7031A1171BF86BE196C857752C095AFEE8CEE6E4740E61032001F0087004D4F44303941312E41323030303034392E6831327630342E3030352E323030363236383138343034315F7375725F7265666C5F6230312E74696600'::raster,'MOD09A1.A249.h12v04.005.2006268184041_sur_refl_b01.tif'); INSERT INTO "mod09a1.a249.h12v04.005.2006268184041_sur_refl_b01" ("rast","filename") VALUES ('010100A544B7031A11713FA544B7031A1171BF59D7C60F306A52C095AFEE8CEE6E4740E61032001F0087004D4F44303941312E41323030303034392E6831327630342E3030352E323030363236383138343034315F7375725F7265666C5F6230312E74696600'::raster,'MOD09A1.A249.h12v04.005.2006268184041_sur_refl_b01.tif'); INSERT INTO "mod09a1.a249.h12v04.005.2006268184041_sur_refl_b01" ("rast","filename") VALUES ('010100A544B7031A11713FA544B7031A1171BF2BF073B3DA5C52C095AFEE8CEE6E4740E61032001F0087004D4F44303941312E41323030303034392E6831327630342E3030352E323030363236383138343034315F7375725F7265666C5F6230312E74696600'::raster,'MOD09A1.A249.h12v04.005.2006268184041_sur_refl_b01.tif'); INSERT INTO "mod09a1.a249.h12v04.005.20062 Or, this one (same tif file; I only shortened the filename) breaks in the middle of a cell (or tile?) representation, still at line 168: ... INTO "mod09a1.a249.sur_refl_b01" ("rast","filename") VALUES ('010100A544B7031A11713FA544B7031A1171BF86BE196C857752C095AFEE8CEE6E4740E61032001F0087004D4F44303941312E41323030303034392E7375725F7265666C5F6230312E74696600'::raster,'MOD09A1.A249.sur_refl_b01.tif'); INSERT INTO "mod09a1.a249.sur_refl_b01" ("rast","filename") VALUES ('010100A544B7031A11713FA544B7031A1171BF59D7C60F306A52C095AFEE8CEE6E4740E61032001F0087004D4F44303941312E41323030303034392E7375725F7265666C5F6230312E74696600'::raster,'MOD09A1.A249.sur_refl_b01.tif'); INSERT INTO "mod09a1.a249.sur_refl_b01" ("rast","filename") VALUES ('010100A544B7031A11713FA544B7031A1171BF2BF073B3DA5C52C095AFEE8CEE6E4740E61032001F0087004D4F44303941312E41323030303034392E7375725F7265666C5F6230312E74696600'::raster,'MOD09A1.A249.sur_refl_b01.tif'); INSERT INTO "mod09a1.a249.sur_refl_b01" ("rast","filename") VALUES ('010100A544B7031A11713FA544B7031A117 This is VERY frustrating and I hope an answer will be found! Guillaume Message: 2 Date: Thu, 24 Apr 2014 21:43:43 +0200 From: Mateusz ?oskot To: PostGIS Users Discussion Subject: Re: [postgis-users] Subject: Re: raster2pgsql 2.1.1 on Windows 7 64bit - ERREUR: erreur de syntaxe Message-ID: Content-Type: text/plain; charset=UTF-8 On 24 April 2014 21:19, Guillaume Drolet wrote: > I tried the same command as in my previous post on my PC at home. It's > Windows 7 64bit as well and the cluster has the same locale > ('French.Canadian.1252'.) and encoding (UTF8). > > The command was successful so my work colleague and I, who experiences > the same behaviour as the one described previously on this list, > suspect that it may be due to some problems or conflict between > PostGIS and the hardware (we have similar systems): is this a sensible > hypothesis? I doubt it. I use raster2pgsql frequently on Windows, since PostGIS pre-2.0 and I have never experienced any such issues. raster2pgsql is a standalone command line utility that generates SQL commands and outputs them to stdout. raster2pgsql has no run-time dependency on PostgreSQL or PostGIS. First thing to check is: does raster2pgsql generate valid output, if redirected to file, does the SQL commands in output file look well-formed and complete? If not, which SQL lines are broken? Do the broken lines correspond to any character data in non-ASCII encodings? For example, check if you run raster2pgsql with -F switch to add column with the name of your raster file and the raster file name contains characters with French accents, etc. I'd suggest, don't use psql.exe or load any data to PostgreSQL until you confirm raster2pgsql generates text output with complete well-formed and valid SQL commands first. Best regards, -- Mateusz ?oskot, http://mateusz.loskot.net ___ postgis-users mailing list postgis-users@lists.osgeo.org http://lists.osgeo.org/cgi-bin/mailman/listinfo/postgis-users
[postgis-users] Subject: Re: raster2pgsql 2.1.1 on Windows 7 64bit - ERREUR: erreur de syntaxe
Hi Regina and others: We found the problem with raster2pgsql! It was my version of GDAL (1.11dev) that my GDAL environment variable was pointing to that created the bug. I simply moved my GDAL folder to some other location (without changing GDAL environment variable) and ran raster2pgsql again so that PostGIS could use its own GDAL instead. This resulted in a complete (i.e. no broken line at the end) output SQL file from raster2pgsql. Since I want to keep using my GDAL version (1.11dev, with all my add-ons compiled and the version used by other software), this raises the question: how to force PostGIS to not use the GDAL installation referenced by the GDAL environment variable but instead use its own GDAL? Thanks! G 2014-04-30 15:00 GMT-04:00 : > Send postgis-users mailing list submissions to > postgis-users@lists.osgeo.org > > To subscribe or unsubscribe via the World Wide Web, visit > http://lists.osgeo.org/cgi-bin/mailman/listinfo/postgis-users > or, via email, send a message with subject or body 'help' to > postgis-users-requ...@lists.osgeo.org > > You can reach the person managing the list at > postgis-users-ow...@lists.osgeo.org > > When replying, please edit your Subject line so it is more specific > than "Re: Contents of postgis-users digest..." > > > Today's Topics: > >1. Re: Subject: Re: raster2pgsql 2.1.1 on Windows 7 64bit - > ERREUR: erreur de syntaxe (Paragon Corporation) > > > -- > > Message: 1 > Date: Tue, 29 Apr 2014 18:51:16 -0400 > From: "Paragon Corporation" > To: "'PostGIS Users Discussion'" > Subject: Re: [postgis-users] Subject: Re: raster2pgsql 2.1.1 on > Windows 7 64bit - ERREUR: erreur de syntaxe > Message-ID: <40C3ECB89B9F4C388A6CE0F390D877C5@O> > Content-Type: text/plain; charset="us-ascii" > > Guillaume, > > Did you say this works fine on your home pc, but not your work one? > > I'm wondering if maybe its some malware or virus checker shutting down the > execution of raster2pgsql for some reason. > > Have you looked at the windows event viewer to see if it is showing any > messages? I recall one time a behind the scenes VMWare virus checker kept > on deleting some of my Msys files because it thought it was hacker ware. No > duh. Had to tell sysadmins to exclude the folder from their checking. > > Hope that helps, > Regina > > -Original Message- > From: postgis-users-boun...@lists.osgeo.org > [mailto:postgis-users-boun...@lists.osgeo.org] On Behalf Of Guillaume Drolet > Sent: Monday, April 28, 2014 2:30 PM > To: postgis-users@lists.osgeo.org > Subject: [postgis-users] Subject: Re: raster2pgsql 2.1.1 on Windows 7 64bit > - ERREUR: erreur de syntaxe > > Hi Mateusz, > > Thanks for helping. What you suggested, i.e. outputting raster2pgsql to an > SQL file instead of piping it to psql, I have done already. I also posted > examples of where the lines are broken in previous post (see also > Jean-Daniel Sylvain's older post on the same topic). > > It doesn't seem to break on non-ASCII characters. For example, this one > breaks in the middle of the file name at line 168: > > ... > INSERT INTO "mod09a1.a249.h12v04.005.2006268184041_sur_refl_b01" > ("rast","filename") VALUES > ('010100A544B7031A11713FA544B7031A1171BF86BE196C857752C095AFEE8CEE6E4740 > E61032001F0087004D4F44303941312E > 41323030303034392E6831327630342E3030352E323030363236383138343034315F7375725F > 7265666C5F6230312E74696600'::raster,'MOD09A1.A249.h12v04.005.20062681840 > 41_sur_refl_b01.tif'); > INSERT INTO "mod09a1.a249.h12v04.005.2006268184041_sur_refl_b01" > ("rast","filename") VALUES > ('010100A544B7031A11713FA544B7031A1171BF59D7C60F306A52C095AFEE8CEE6E4740 > E61032001F0087004D4F44303941312E > 41323030303034392E6831327630342E3030352E323030363236383138343034315F7375725F > 7265666C5F6230312E74696600'::raster,'MOD09A1.A249.h12v04.005.20062681840 > 41_sur_refl_b01.tif'); > INSERT INTO "mod09a1.a249.h12v04.005.2006268184041_sur_refl_b01" > ("rast","filename") VALUES > ('010100A544B7031A11713FA544B7031A1171BF2BF073B3DA5C52C095AFEE8CEE6E4740 > E61032001F0087004D4F44303941312E > 41323030303034392E6831327630342E3030352E323030363236383138343034315F7375725F > 7265666C5F6230312E74696600'::raster,'
[postgis-users] raster2pgsql 2.1.1 on Windows 7 64bit - ERREUR: erreur de syntaxe
Regina, When I do SELECT postgis_full_version(); It says (from what I remember. I won't have access to my database until next week) : GDAL 1.11dev (of course if my 1.11dev installation is at the place described by the GDAL variable). If my GDAL installation is moved as I did in my previous post, the GDAL version shown by SELECT postgis_full_version(); is 1.10. So it looks like a problem with raster2pgsql and rather than with PostGIS as the loading of the SQL file from raster2pgsql runs fine. I don't know who maintains this program but maybe it'd be worth checking it against versions of GDAL after 1.10. Regards, Guillaume > > Message: 2 > Date: Wed, 30 Apr 2014 21:24:06 -0400 > From: "Paragon Corporation" > To: "'PostGIS Users Discussion'" > Subject: Re: [postgis-users] Subject: Re: raster2pgsql 2.1.1 on > Windows 7 64bit - ERREUR: erreur de syntaxe > Message-ID: > Content-Type: text/plain; charset="us-ascii" > > Guillaume, > > It was my understanding that PostGIS only uses the libgdal in the same > folder as raster2pgsql since that would be the first in its path. However it > uses the environment variable for the proj files and so forth. So I'm still > a bit surprised this is an issue unless those files are vastly different > from 1.10. > > When you do > > SELECT postgis_full_version(); > > What does it give for the GDAL version. > > > If its just a during load problem, I think you can just simply override the > environment variable in your loader batch script > > SET GDAL_DATA="C:\Program Files\PostgreSQL\9.3\gdal-data" > > Of course this won't fix your running install. In PostGIS 2.2 we introduced > a PostgreSQL GUC variable to handle this: > > http://postgis.net/docs/manual-dev/postgis_gdal_datapath.html > > But unfortunately that won't help you for 2.1 and not sure how much that > helps with raster2pgsql anyway since I don't think raster2pgsql makes > connection to Postgres to be able to even see that variable. > > What I could do is create a build that uses the newer GDAL to see if that > improves your situation (you can be my test guinea pig :) ). IN 2.2 I plan > to ship with GDAL 1.11 and probably upgrade next 2.1 to GDAL 1.11 as well > since unfortunately they have to share the same GDAL on same PostgreSQL > install and some people I expect will be running both. > > > Opinions? > Regina > http://www.postgis.us > > > -Original Message- > From: postgis-users-boun...@lists.osgeo.org > [mailto:postgis-users-boun...@lists.osgeo.org] On Behalf Of Guillaume Drolet > Sent: Wednesday, April 30, 2014 4:39 PM > To: postgis-users@lists.osgeo.org > Subject: [postgis-users] Subject: Re: raster2pgsql 2.1.1 on Windows 7 64bit > - ERREUR: erreur de syntaxe > > Hi Regina and others: > > We found the problem with raster2pgsql! It was my version of GDAL > (1.11dev) that my GDAL environment variable was pointing to that created the > bug. I simply moved my GDAL folder to some other location (without changing > GDAL environment variable) and ran raster2pgsql again so that PostGIS could > use its own GDAL instead. This resulted in a complete (i.e. no broken line > at the end) output SQL file from raster2pgsql. > > Since I want to keep using my GDAL version (1.11dev, with all my add-ons > compiled and the version used by other software), this raises the question: > how to force PostGIS to not use the GDAL installation referenced by the GDAL > environment variable but instead use its own GDAL? > > Thanks! > > G > > 2014-04-30 15:00 GMT-04:00 : >> Send postgis-users mailing list submissions to >> postgis-users@lists.osgeo.org >> >> To subscribe or unsubscribe via the World Wide Web, visit >> http://lists.osgeo.org/cgi-bin/mailman/listinfo/postgis-users >> or, via email, send a message with subject or body 'help' to >> postgis-users-requ...@lists.osgeo.org >> >> You can reach the person managing the list at >> postgis-users-ow...@lists.osgeo.org >> >> When replying, please edit your Subject line so it is more specific >> than "Re: Contents of postgis-users digest..." >> >> >> Today's Topics: >> >>1. Re: Subject: Re: raster2pgsql 2.1.1 on Windows 7 64bit - >> ERREUR: erreur de syntaxe (Paragon Corporation) >> >> >> -- >> >> Message: 1 >> Date: Tue, 29 Apr 2014 18:51:16 -0400 >> From: "Paragon Corporation" >> To: "'PostGIS Users Discussion'&quo
[postgis-users] Topology: cannot delete slivers (or gaps)
Hi! I'm am turning to you all hoping you'll have some tricks to help with this problem: I created a topology for 2290 polygons and populated a topogeom column using these commands: SELECT topology.CreateTopology('de_20k_topo', find_srid('syshiera', 'de_20k', 'geom_32198'), 1e-4); SELECT topology.AddTopoGeometryColumn('de_20k_topo', 'syshiera', 'de_20k', 'topogeom', 'MULTIPOLYGON'); DO $$DECLARE r record; BEGIN FOR r IN SELECT * FROM syshiera.de_20k LOOP BEGIN UPDATE syshiera.de_20k SET topogeom = toTopoGeom(geom_32198, 'de_20k_topo', 1, 1e-4) WHERE gid = r.gid; EXCEPTION WHEN OTHERS THEN RAISE WARNING 'Loading of record % failed: %', r.gid, SQLERRM; END; END LOOP; END$$; -- For testing: SELECT topology.CopyTopology('de_20k_topo', 'de_20k_topo_tests'); I need to fix 1) extra nodes along edges that are not at an intersection and 2) sliver polygons. For cases in 1, I was successful using 'ST_NewEdgeHeal' and a plPgsql script. For cases in 2, I can remove some slivers with 'ST_RemEdgeModFace' but some of them I just can't delete: I get either: ERROR: TopoGeom 1123 in layer 1 (de_20k_topo_tests.LAYER1.) cannot be represented healing faces 1514 and 1502 Or, when trying with remove a node with 'topology.ST_NewEdgeHeal': ERROR: SQL/MM Spatial exception - other edges connected (3912) In this particular example, two edges are connected two the same two nodes, similar to what was discussed here: http://postgis.17.x6.nabble.com/PostGIS-1955-Exception-when-2-edges-passed-to-ST-ModEdgeHeal-that-are-attached-to-the-same-2-nodes-td4999383.html But according to that discussion thread, this was fixed in an earlier version (I'm using topology version 2.1.3). Maybe my problem is different but I'll need your help to get it sorted. Here's an example showing the small face (1502) I want remove: https://dl.dropboxusercontent.com/u/5196336/example1.bmp Thanks for your help! ___ postgis-users mailing list postgis-users@lists.osgeo.org http://lists.osgeo.org/cgi-bin/mailman/listinfo/postgis-users
Re: [postgis-users] Topology: cannot delete slivers (or gaps)
Hi Sandro, Thanks for your help! topology.ValidateTopology gives this: error id1 id2 character varying integer integer -- --- face has no rings 1831 If I calculate the area of face 1831 (st_area(st_getfacegeometry('de_20k_topo', face_id))), I get a NULL value. I also find I have some 650 of these slivers with an area below 75 m2 while most of my correct polygons have an area above 300 000 000 m2. Re the edges passed to ST_NewEdgeHeal, in the example I sent it was edges 3908 and 3909. If you think about something else that could be problematic, let me know. Meanwhile I'll dig into the approach you suggested in your reply (i.e. fixing TopoGeometry objects). Thanks! PS. Forgot to mention that in the first topology I had created, I used a tolerance of 2.0 meters instead of 1e-4 meters. While that yielded much less slivers, I was getting way more of those nodes not located at an intersection. What is best? -- View this message in context: http://postgis.17.x6.nabble.com/Topology-cannot-delete-slivers-or-gaps-tp5007250p5007257.html Sent from the PostGIS - User mailing list archive at Nabble.com. ___ postgis-users mailing list postgis-users@lists.osgeo.org http://lists.osgeo.org/cgi-bin/mailman/listinfo/postgis-users
Re: [postgis-users] Topology: cannot delete slivers (or gaps)
Sandro Santilli wrote >> I suggest you take a look at edges having face 1831 on either >> the right or left edge. Do any edge exist ? > > edge 5187 has face 1831 as its right face > >> The correct procedure to drop those slivers would be to assign >> those small areas to the (closest|biggest|leftmost|youpickit) >> TopoGeometry object, to one and only one, and then, if you really >> need it, drop the edge internal to the TopoGeometry. > >> I made a QGIS plugin which may help you do some of these things >> manually, which is useful in some cases, mainly to understand the >> issues. > > The plugin you made is a great tool (thanks). I will first try to get the > number of slivers down to > a manageable number, which I had with setting the tolerance at 2.0 m, and > then apply the > TopoGeometry approach above. > >> So the error was right about there being another edge attached to the >> node: > >> ERROR: SQL/MM Spatial exception - other edges connected (3912) > >> I'm guessing 3912 and 3908 are the two parallel edges in the picture >> you attached: https://dl.dropboxusercontent.com/u/5196336/example1.bmp > > You're right, they are the two parallel edges. And it's the same for most > of > the some 650 slivers I have. > >> Again, QGIS may help you getting more info out of the visual, enable >> the DBManager core plugin, select your topology schema and select >> Schema->TopoViewer from the menu. > > I didn't know about the TopoViewer tool. I was importing the topology > tables (node, edge, face) > like any other PG geometry tables. Will explore with this new tool. > > >> It's probably easier to clean up nodes for an areal topology than to >> clean up faces. With my QGIS plugin (PostGIS Topology Editor) you can >> select all nodes and click on a button which will drop all the ones >> that are safe to drop. > > You're right! I'll create my topology again but with the 2.0 m tolerance, > then run > my function for getting rid of nodes that are not at intersections and > finally, > I will use the approach you suggested of working with the TopoGeometry > objects and the > relation table to eliminate the remaining slivers. > > Will keep you posted on how it goes. Again, I really appreciate your help. > Thanks. > > Guillaume. -- View this message in context: http://postgis.17.x6.nabble.com/Topology-cannot-delete-slivers-or-gaps-tp5007250p5007266.html Sent from the PostGIS - User mailing list archive at Nabble.com. ___ postgis-users mailing list postgis-users@lists.osgeo.org http://lists.osgeo.org/cgi-bin/mailman/listinfo/postgis-users
Re: [postgis-users] Topology: cannot delete slivers (or gaps)
Hi Sandro and others, First, thanks Jon for the hint on pprepair: I ran pprepair on my shapefile and I'm now in the process of importing it into PostGIs for topology creation. Then I'll be able to tell if pprepair really fixed everything. But because I'm stubborn and never give up, I still want to succeed on fixing TopoGeometry objects in PostGIS (bear with me, what follows is a bit long): Sandro, you said: > The fix here would be to find faces that take part in the composition > of multiple TopoGeometry objects and those NOT taking part in the > composition of any TopoGeometry object. I've seen GRASS GUIs referring > to those 2 classes as "polygon0" and "polygon2", enough that I've been > thinking it could be useful to have them shown as layers in the QGIS > viewer for PostGIS Topology. You should be able to extract those faces > with a query to the relation and the layers tables. > So basically you first have to fix the TopoGeometry objects and then > can safely drop the no-more-needed edges. First, I looked at faces NOT taking part in the composition of any TopoGeometry object: /SELECT face_id FROM de_20k_topo.face WHERE face_id NOT IN (SELECT element_id FROM de_20k_topo.relation);/ face_id 0 2418 2650 2663 2766 3077 3118 3146 3267 3270 None of these 10 faces are involved with my three remaining slivers... And when I look at them, they ARE related to edges, as a next_face or a right_face. Is there any useful information in that for my purpose? Second, I tried to work along the lines suggested by Sandro, using a new topology (tolerance 2.0 m). As in one of my previous attempts, I successfully fixed 20 of the 23 slivers using ST_RemEdgeModFace. For the remaining three slivers, however, I tried modyfying TopoGeometry objects but with no success. I must admit I'm a bit lost and not too sure if I really understand what TopoGeometry objects are: For example, let's look at problematic face_id: 2227 (https://dl.dropboxusercontent.com/u/5196336/example_2227.bmp). In this case, the light grey line is where a single edge should be. Blue edges 6381 and 6341 should not be there. The light grey line should start from the node at the bottom right and connect at its other end somewhere in the middle of edge 6340: I thought, let's first remove edge 6346, which should not exist in the first place, and then split edge 6340 with a node and connect it to node 4753 on the lower right corner of the image: SELECT ST_RemEdgeModFace('de_20k_topo', 6346); As in my previous post, this yields: ERROR: TopoGeom 1546 in layer 1 (syshiera.de_20k.topogeom) cannot be represented healing faces 2225 and 2227 CONTEXT: instruction SQL « SELECT topology._ST_RemEdgeCheck(toponame, topoid, e1id, e1rec.left_face, e1rec.right_face) » function PL/pgsql st_remedgemodface(character varying,integer), line 43 at PERFORM So I looked at TopoGeom 1546: /SELECT * FROM de_20k_topo.relation WHERE topogeo_id = 1546;/ topogeo_id;layer_id;element_id;element_type 1546;1;2225;3 1546;1;2233;3 Edge 2227 is missing from TopoGeom 1546: Where is it? /SELECT * FROM de_20k_topo.relation WHERE element_id = 2227;/ topogeo_id;layer_id;element_id;element_type 1551;1;2227;3 There it is, in TopoGeom 1551. What else is there: /SELECT * FROM de_20k_topo.relation WHERE topogeo_id = 1551;/ topogeo_id;layer_id;element_id;element_type 1551;1;2225;3 1551;1;2234;3 1551;1;2227;3 So edges 2225 and 2227 are associated with TopoGeom 1551. Will change that: / UPDATE de_20k_topo.relation SET topogeo_id = 1546 WHERE topogeo_id = 1551 AND element_id = 2227;/ Query returned successfully: one row affected, 481 ms execution time. And then again: /SELECT ST_RemEdgeModFace('de_20k_topo', 6346); / ERROR: TopoGeom 1551 in layer 1 (syshiera.de_20k.topogeom) cannot be represented healing faces 2225 and 2227 CONTEXT: instruction SQL « SELECT topology._ST_RemEdgeCheck(toponame, topoid, e1id, e1rec.left_face, e1rec.right_face) » function PL/pgsql st_remedgemodface(character varying,integer), line 43 at PERFORM Obviously, as mentionned above, I don't quite understand how TopoGeometry objects work and the link between TopoGeoms, geometries and primitives. I need a bit of education here. Thanks a lot for your patience and your help. PS. Remi, thanks for suggesting GRASS but I've never used it and the learning curve seems pretty steep, no? -- View this message in context: http://postgis.17.x6.nabble.com/Topology-cannot-delete-slivers-or-gaps-tp5007250p5007273.html Sent from the PostGIS - User mailing list archive at Nabble.com. ___ postgis-users mailing list postgis-users@lists.osgeo.org http://lists.osgeo.org/cgi-bin/mailman/listinfo/postgis-users
Re: [postgis-users] Topology: cannot delete slivers (or gaps)
Sandro Santilli wrote >> Except face 0, which is the "universe" face, all the others >> should be "gaps" in your topology. You can zoom to the face's >> bounding box ("mbr" field) to take a closer look (do you have qgis >> set up by now?). > > Ok. I had already looked at some of them but as face-derived geometries > (using ST_GetFaceGeometry) and not as real faces. As geometries, to me > they look like polygons and I didn't realize they were in fact holes. See > this example in QGIS for face_id 2418 turned into a polygon: > https://dl.dropboxusercontent.com/u/5196336/geomFromFace_id_2418.bmp. Now, > look at its mbr representation: > https://dl.dropboxusercontent.com/u/5196336/mbrFace_id_2418.bmp. How do > you tell it's a hole and not just a regular face? > >> What you could do programmatically is for each face to fetch the list >> of edges bounding it and pick the face on the other side as a candidate; >> then for each candidate face you find TopoGeometry objects containing >> it in their definition and pick one to assign that orphaned face >> (you could use smaller TopoGeometry, or larger, or whatever). > > Let's try this or face_id 2418 (which is, btw, a meaningful unit and not > an error...): / > SELECT f.face_id, e.edge_id, e.left_face, e.right_face > FROM de_20k_topo.face f JOIN > de_20k_topo.edge_data e ON (f.face_id = e.left_face OR >f.face_id = e.right_face) > WHERE face_id = 2418; / > * > face_id, edge_id, left_face, right_face * > 2418;2663;2418;950 > 2418;6986;2418;340 > 2418;2672;2418;951 > 2418;2665;2418;344 > 2418;6989;2418;2417 > > I pick the first edge (2663) and take it's right_face, 950, as a > candidate: / > SELECT topogeo_id > FROM de_20k_topo.relation > WHERE element_id = 950; / > * > topogeo_id * > 814 > > Only TopoGeom 814 has face_id 950 as one of its elements. Will assign > face_id 2418 > to this TopoGeom: / > INSERT INTO de_20k_topo.relation > VALUES (814, 1, 2418, 3); / > > Then what? Is that it? Does that mean that that big hole 2418 is now > unioned with face_id 950 (in which case it isn't what I want as it's a > meaningful ecological unit)? > > >> It might be better first to deal with the faces that are associated >> with multiple TopoGeometry, just to avoid assigning an orphaned face >> to a TopoGeometry based on it's containing a face which will later >> be removed from the definition because duplicated. > > Good point. > >> Eh, that part is a still a bit weak in documantation, maybe. >> A TopoGeometry is a Geometry defined by its components. > >> In your case (areal TopoGeometry objects) it is a (multi)Polygon >> defined by a set of faces. Note that QGIS is also able to >> show you the TopoGeometry objects, and lets you edit them "spatially" >> (less robust than editing them by adding/removing components but with >> some care easier to use). > >>> In this case, the light grey line is where a single edge should be. Blue >>> edges 6381 >>> and 6341 should not be there. The light grey line should start from the >>> node >>> at the bottom right and connect at its other end somewhere in the middle >>> of >>> edge 6340: > >> Uhm, edge identifiers are missing from that dump. >> I take it 6381 is the lower blue and 6341 the upper blue. > > You're correct here. > >> I suggest you load the TopoGeometry layer with qgis and see for yourself. > > When you say load the TopoGeometry layer, I assume you mean the topogeom > column I populated with: / > UPDATE syshiera.de_20k SET topogeom = toTopoGeom(geom_32198, > 'de_20k_topo', 1, 2.0) / > > Right? > > When I load it in QGIS using the PostGIS table loader, I can visualize it > and see there are holes. E.g. > https://dl.dropboxusercontent.com/u/5196336/holeInTopogeometry.bmp > >> Unfortunately (that'd be useful to add in qgis) it won't be easy to >> find TopoGeometry 1546 as extracting the id of a TopoGeometry requires >> running an expression on the database: id(topogeom). > > I can either / > SELECT * > FROM syshiera.de_20k > WHERE id(topogeom) = 1546; / > > or, as I mentioned before, I can superimpose my 'faces' layer, which I > derived from de_20k_topo.face > using ST_GetFaceGeometry, to 'topogeom' to get the topogeo_id. No? > >> You could create a view selecting just that topoGeometry, or a subset >> of them, or you could just load them all, hoping it won't be too >> expensive >> (another spot where qgis support for postgis topology could be improved). > > I agree, just loading the TopoViewer freezes QGIS for a while and I don't > lack ressources on my machine. > >> Chances are it would be enough to assign face 2227 to the TopoGeometry >> with id 1546 (likely the one not-assigned). But as face 2227 was not >> in your initial set (unassigned faces) I guess there's another >> TopoGeometry >> which does have 2227 in his definition, but not the one on the right. >> In that case you have to pick which of the two contending TopoGeometry >> objects should get face 2227 assign
Re: [postgis-users] Topology: cannot delete slivers (or gaps)
Hi John, I tried pprepair and imported the fixed shapefile into PostGIS. Although pprepair said all polygons were fixed, there are still 22 slivers in my topology. Maybe due to the value I used for tolerance? Finding the right tolerance for building a topology isn't an easy task... -- View this message in context: http://postgis.17.x6.nabble.com/Topology-cannot-delete-slivers-or-gaps-tp5007250p5007279.html Sent from the PostGIS - User mailing list archive at Nabble.com. ___ postgis-users mailing list postgis-users@lists.osgeo.org http://lists.osgeo.org/cgi-bin/mailman/listinfo/postgis-users
Re: [postgis-users] Topology: cannot delete slivers (or gaps)
Sandro, Using your help, I filled most holes but three I couldn't fill using the "UPDATE ... SET topogeom ..." approach. Will explore other approaches later. Re the presence of holes, I may have deleted faces when running my function for removing extra nodes. This needs checking. Advancing in my work, I also got to drop edge 6341 and I end up with this: https://dl.dropboxusercontent.com/u/5196336/edge6340.bmp So the next step is to split edge 6340 about its middle, with the goal of then connecting the newly created node with node 4753, thus creating an edge where there should be one (i.e. on the light grey line). However, I get the error: *point not on edge* when running ST_ModEdgeSplit. First, I tried to find the new node location using the coordinate capture tool in QGIS but because I got that error I thought maybe I really need to find a point that's on the edge so I used this approach: /WITH edge AS (SELECT geom line FROM de_20k_topo.edge_data WHERE edge_id = 6340 ), point AS (SELECT e.line, ST_ClosestPoint(e.line, ST_GeomFromText('POINT(-208737.346 811681.163)', 32198)) point FROM edge e ) SELECT ST_ModEdgeSplit('de_20k_topo', 6340, point) FROM point; / I still get the same error: point not on node. How is one supposed to find a point that IS on an edge? -- View this message in context: http://postgis.17.x6.nabble.com/Topology-cannot-delete-slivers-or-gaps-tp5007250p5007281.html Sent from the PostGIS - User mailing list archive at Nabble.com. ___ postgis-users mailing list postgis-users@lists.osgeo.org http://lists.osgeo.org/cgi-bin/mailman/listinfo/postgis-users
Re: [postgis-users] Topology: cannot delete slivers (or gaps)
Sandro: TopoGeo_addPoint sounds like the appropriate choice. I hadn't seen this alternative topology function. Will try that, thanks. Remi: Thanks for your advice on when to use what tools. I should give GRASS a shot at topology cleaning. Since topology-related work isn't something I do everyday and PostGIS is a tool I use frequently, I naturally went in that direction. -- View this message in context: http://postgis.17.x6.nabble.com/Topology-cannot-delete-slivers-or-gaps-tp5007250p5007284.html Sent from the PostGIS - User mailing list archive at Nabble.com. ___ postgis-users mailing list postgis-users@lists.osgeo.org http://lists.osgeo.org/cgi-bin/mailman/listinfo/postgis-users
Re: [postgis-users] Topology: cannot delete slivers (or gaps)
I finally finished cleaning my topology, using a mix of SQL commands, functions, and also using Sandro's QGIS topology editor plugin. In parallel, I also tried going the GRASS route, importing my original shapefile in GRASS and applying the cleaning function (v.clean ...) as suggested by Remi. I injected the result into PostGIS and created a topology. Most of the 23 slivers were there in the topology so one of two things: 1) GRASS didn't fix them (maybe I didn't choose an appropriate tolerance) or 2) they were created by PostGIS when building the topology. Anyway, I got what I finally got what I need, thanks a lot to all of you who helped me. I'm gonna ask you one last thing if I may: I want to replace the original geom column (i.e. that used to build the topogeom) with a topologically correct one, and keep all the associated ecological attributes. Is this the right way to do it: UPDATE syshiera.de_20k SET geom = topogeom::geometry; ? Thanks. -- View this message in context: http://postgis.17.x6.nabble.com/Topology-cannot-delete-slivers-or-gaps-tp5007250p5007286.html Sent from the PostGIS - User mailing list archive at Nabble.com. ___ postgis-users mailing list postgis-users@lists.osgeo.org http://lists.osgeo.org/cgi-bin/mailman/listinfo/postgis-users
[postgis-users] Dynamic parameters to ST_MapAlgebra
Hi, I'm trying to find the best way to dynamically apply scaling factors to a raster. I've tried different approaches using ST_MapAlgebra, both callback and expression versions. The fastest method so far is this one but it doesn't allow for dynamic parameters: SELECT ST_MapAlgebra(rast, 1, '32BF', '([rast] * 0.01) + 0') AS rast, FROM evi; Compared to using the following function, the query above is ten times faster: CREATE OR REPLACE FUNCTION public.scale_sds_plv8(value double precision[][][], pos integer[][], VARIADIC userargs text[]) RETURNS double precision AS $$ var g = Number(userargs[0]); var o = Number(userargs[1]); return (value * g + o); $$ LANGUAGE plv8 IMMUTABLE; WITH scaling_params AS ( SELECT gain, off_set FROM meta WHERE product = 'MCD12Q2' AND sds = 'EVI' ) SELECT ST_MapAlgebra(rast, 1, 'scale_sds_plv8(double precision[], integer[], text[])'::regprocedure, '32BF', 'FIRST', NULL::raster, 0, 0, VARIADIC ARRAY[gain, off_set]::text[]) AS rast FROM evi, scaling_params I want to take advantage of the speed of the first method above but be able to pass parameters dynamically instead of hard-coding them in the query. To this aim I tried this approach: CREATE OR REPLACE FUNCTION scale_sds_plpgsql(rast raster, gain float8, offs float8 DEFAULT 0) RETURNS raster AS $$ DECLARE sql text; BEGIN sql := 'SELECT ST_MapAlgebra(' || rast || ', 1, ''32BF'', ''([rast] * ' || gain || ') + ' || offs || ''')'; EXECUTE sql; END; $$ LANGUAGE plpgsql IMMUTABLE; SELECT scale_sds_plpgsql(rast, 0.01, 0.0) FROM evi; Running this last query, I get a synthax error which I don't understand: ERREUR: erreur de syntaxe sur ou près de « F4000407FC05E10B2668C203F41A0BD86AA1F7C2 CONTEXT: fonction PL/pgsql scale_rast(raster,double precision,double precision), ligne 12 à instruction EXECUTE ** Error ** I'm sure there are experienced folks on this list who will have some clues to solve this or who will come up with a better approach. Thanks a lot for your time. G -- View this message in context: http://postgis.17.x6.nabble.com/Dynamic-parameters-to-ST-MapAlgebra-tp5007468.html Sent from the PostGIS - User mailing list archive at Nabble.com. ___ postgis-users mailing list postgis-users@lists.osgeo.org http://lists.osgeo.org/cgi-bin/mailman/listinfo/postgis-users
Re: [postgis-users] Dynamic parameters to ST_MapAlgebra
Hugues: thank a lot for your fix. The function now runs without error... but it's been running for almost an hour now :( Compared to the 30-ish seconds it takes using the hard-coded query, and the some 120 seconds using the PLV8 version, I cannot think of using the PLPGSQL version in its current formulation. Does anyone else has a possibly faster suggestion, more in the ballpark of the two other methods? Thanks, G -- View this message in context: http://postgis.17.x6.nabble.com/Dynamic-parameters-to-ST-MapAlgebra-tp5007468p5007470.html Sent from the PostGIS - User mailing list archive at Nabble.com. ___ postgis-users mailing list postgis-users@lists.osgeo.org http://lists.osgeo.org/cgi-bin/mailman/listinfo/postgis-users
Re: [postgis-users] Dynamic parameters to ST_MapAlgebra
Hi Regina, " Paragon Corporation-2 wrote > Guillaume, > > If the expression approach is faster, is there some reason you aren't just > doing > > > WITH scaling_params AS ( > SELECT gain::text As gain, off_set::text as off_set > FROM meta > WHERE product = 'MCD12Q2' AND sds = 'EVI' > ) > SELECT ST_MapAlgebra(rast, 1, '32BF', '([rast] * ' || gain || ') + ' > || > off_set ) AS rast > FROM evi, scaling_params; > > Or is that too hard-coded for you? No it isn't. My SQL skills are just not enough sharp yet. I can't say I have loads of experience with this language and some of the possibilities still slip my mind. Thanks for the enlightening me! > Now regarding the one that Hugues gave you, you have to run it differently > from how you were calling your others, otherwise it would take forever > since > for every raster row you have it would apply map algebra to all the rows > in > your table. So how did you call his? Called it pretty much the same way as I showed in my first post: SELECT scale_sds_plpgsql(rast, 0.01, 0.0) FROM evi WHERE date = mydate; > If you wanted to use the syntax you have I would rewrite your plpgsql as > follows (and not use dynamic sql) - you don't need dynamic sql to do this > just call ST_MapAlgebra directly. When I referred to 'dynamic', I meant with passing parameters resulting from a query. Maybe 'dynamic' isn't the right word? Thanks for all your suggestions, I appreciate it. Will try them and see how they improve things. Best wishes, G -- View this message in context: http://postgis.17.x6.nabble.com/Dynamic-parameters-to-ST-MapAlgebra-tp5007468p5007475.html Sent from the PostGIS - User mailing list archive at Nabble.com. ___ postgis-users mailing list postgis-users@lists.osgeo.org http://lists.osgeo.org/cgi-bin/mailman/listinfo/postgis-users