Re: [postgis-users] Getting TopologyExections when trying to node linestrings to create an overlay
On Mon, Feb 9, 2015 at 6:25 AM, Sandro Santilli s...@keybit.net wrote: First of all I confirm it still happens with GEOS=3.5.0dev-CAPI-1.9.0 r4038. Second, I took a look at a random set (geom_set_id=1) and I found it pretty big. That's to say you could probably further reduce the dataset for the ticket. That set contains 109 polygons, I can get the error by attempting to union the boundaries of the first 40 in that set, and I'm sure you can further reduce the input. Thanks for taking a look. I'll work on doing that when I can find the time, but I don't expect that to be a fast process at all. Even just checking for the pairwise case took a decent amount of time to develop the query and took overnight to finish running (even with optimizations like a bounding box intersection test). I don't really have any good heuristic that could narrow down the possibilies for reproducing, so I don't see much option other than to brute force it possibly with some kind of filter. That's why I didn't put more effort into shrinking the input set to begin with. Is a PostGIS database dump an okay format to provide the shapes, or would you prefer something else? I suppose I could dump groups into shapefiles or something like that if it's more convenient. On Mon, Feb 9, 2015 at 7:00 AM, Rémi Cura remi.c...@gmail.com wrote: this is a precision related issue, coordinates are way too big and should be translated. I understand what you mean by that (as in floating point problems due to size), but these coordinates are pretty typical. This is a standard UTM projection, zone 15N in middle America. It's even predefined in PostGIS' list of spatial references. I'm told ESRI had this kind of problem years ago, but they dealt with it as far as I know. While I would choose PostGIS over ESRI any day, this could be viewed by my coworkers as a good argument against using PostGIS; it represents a serious reliability concern since it applies to a broad range of functions and 2D projections. Basically, if you use any projection with coordinates of this size (of which there are a good number), there seems to be no telling when any function will just blow up in your face at random. On Mon, Feb 9, 2015 at 7:04 AM, Rémi Cura remi.c...@gmail.com wrote: Hey, I executed your data, the following command solve the problem (with very recent GEOS for me) (POSTGIS=2.2.0dev r12846 GEOS=3.5.0dev-CAPI-1.9.0 r0 PROJ=Rel. 4.8.0, 6 March 2012 GDAL=GDAL 2.0.0dev, released 2014/04/16 LIBXML=2.8.0 RASTER) [Snip] The change compared to your approach : convert input to table of simple polygons, (no array, no multi). Then translate to improve precision in geos computing Then the union. I don't really understand what you are trying to do, but ist_union seems dangerous and quit ineffective for that . I'm looking at this function call in your code: ST_Union( ST_MakePolygon(ST_ExteriorRing(geom)) ). That call seems to remove all holes and then create the union of all the covered areas (a single multipolygon that covers all areas covered by the originals). That is not what I'm trying to do; I already have an outer boundary that I could use for that purpose. Here's what I'm trying to do. I'm trying to create an overlay (as I said in the first sentence of my original e-mail but didn't elaborate on), in the sense that the term is used by these two articles: http://boundlessgeo.com/2014/10/postgis-training-creating-overlays/ or http://trac.osgeo.org/postgis/wiki/UsersWikiExamplesOverlayTables. What they do is they create a sort of Venn diagram, if you will; they take all the polygons and create a new polygon for each area with a different set of overlappying polygons. Both of them use ST_Union in the same way I am: they take a set of linestrings and union them to create a fully noded linestring, and then they use this linework to create a bunch of new polygons. I'm leveraging ST_Boundary instead of ST_ExteriorRing in my code because I need to preserve holes, but that shouldn't change the results as far as I know. In what way is ST_Union dangerous and quite ineffective for the purpose of noding lingstrings? If that's true, I'm apparently not the only one who has some misconceptions, since both these articles use it the same way. I should also note that I don't show the whole process here; I'm only showing the part where I'm noding the linestrings because that's where this error occurs. Once I have the polygons, I also go back and relate them to the attributed rows in the original tables (of which there are several). That said, I'm trying out the translate and cutting out multi, but I am still using an array in my actual code for a reason. Namely, I need to be able to do this with *different* sets of geometries. These geometries come as the result of selecting polygons from 3 or 4 different tables, each having completely disparate sets of attributes. So basically, I need to be able to use an arbitrary query to get the group of polygons
Re: [postgis-users] operator is not unique: text || geometry
(better to stay on list )
I meant
---
RAISE EXCEPTION '%',_q ;
---
You must understand that plpgsql function fabricate on the fly SQL
statement (meaning, at execution time).
That means that without actually executing the function, there is no way to
know exactly what it does.
Now i I __*can't*__ execute your function, not having your table
Now at execution,
it will stop you function there, and print the UPDATE query that should
have been executed.
Then you can analyse the UPDATE query that have been printed, and test it
to see why it doesn't work and how you could make it work (how which I have
no idea without the query).
You should see something like (I put xxx because I don' have the value.)
-
sql NOTICE :
UPDATE fgcm. SET (x,x,x,x)::topogeometry
= topology.toTopoGeom(ST_Transform(x::geometry,32648),'', 1, 1.0)
WHERE objectid = '
-
Maybe you need to replace the
---
= topology.toTopoGeom(ST_Transform($1::geometry,32648), %I, 1, 1.0)
---
with
-
= topology.toTopoGeom(ST_Transform($1::geometry,32648), %L, 1, 1.0)
-
I can't know.
Cheers,
Rémi-C
2015-02-18 20:43 GMT+01:00 Miller, Stephan smill...@harris.com:
Remi –
I didn’t understand.
Adding RAISE EXCEPTION '%',-q ; before the EXECUTE generates a syntax
error. Did you mean perhaps
RAISE EXCEPTION '%',_q ;
Instead of
EXECUTE _q USING r.shape, cleantopo;
How do I specify the r.shape and cleantopo parameters?
Sorry to be so dense.
Thanks,
Steve
*From:* Rémi Cura [mailto:remi.c...@gmail.com]
*Sent:* Wednesday, February 18, 2015 1:59 PM
*To:* Miller, Stephan
*Subject:* Re: [postgis-users] operator is not unique: text || geometry
As I wrote before,
simply print the update query (and don't execute it)
You can do this by adding a RAISE EXCEPTION '%',-q ; before the EXECUTE
then test it !
Cheers,
Rémi-C
2015-02-18 19:57 GMT+01:00 Miller, Stephan smill...@harris.com:
Remi –
I forced the transform as you suggested using SetSRID. Now I am failing
the UPDATE query somehow.
_
_q := format('SELECT objectid, f_code, shape, topo_shape
FROM fgcm.%I',updatedtablename);
FOR r IN
EXECUTE _q
LOOP
BEGIN
RAISE NOTICE 'Loading % attempt with shape
= % and topo_shape = %' , r.objectid, r.shape, r.topo_shape;
RAISE NOTICE 'Table % Shape %',
updatedtablename, r.topo_shape;
_q :=
format('UPDATE fgcm.%I SET
%I = topology.toTopoGeom(ST_Transform(ST_SetSRID($1, 4326),32648), $2, 1,
1.0)
WHERE objectid = r.objectid' ,updatedtablename, r.topo_shape);
EXECUTE _q USING r.shape, cleantopo;
raise NOTICE 'After % Shape
%',updatedtablename,r.topo_shape;
RAISE NOTICE 'Object % after conversion from shape = %
to topo_shape = %', r.objectid, (ST_AsText(r.shape)),
(ST_AsText(r.topo_shape));
EXCEPTION
WHEN OTHERS THEN
RAISE WARNING 'Loading of record % failed: % %', r.objectid,
SQLSTATE, SQLERRM;
END;
END LOOP;
RETURN;
END
$BODY$
LANGUAGE plpgsql VOLATILE
COST 100
ROWS 2000;
SELECT * FROM fgcm.hc_check_gaps_in_linear_topology('vnroadsclipped',
'VNclippedroadscleantopo');
__
The results for the first feature is shown below.
NOTICE: Loading 1 attempt with shape =
0102E0E6100200380952E7B97B5A40F074DD1774CD3440006AE8C0F87FE825AB94B17B5A40F013885085CD3440006AE8C0F87F
and topo_shape = NULL
NOTICE: Table updatedvnroadsclipped Shape NULL
WARNING: Loading of record 1 failed: 22004 null values cannot be
formatted as an SQL identifier
The absence of the two RAISE NOTICE prints means the UPDATE is failing
somehow. Any suggestions?
Thanks,
Steve
*From:* Rémi Cura [mailto:remi.c...@gmail.com]
*Sent:* Wednesday, February 18, 2015 4:47 AM
*To:* Miller, Stephan
*Subject:* Re: [postgis-users] operator is not unique: text || geometry
Good,
maybe the srid of $1::geometry is not what it should be, you could try
to force it (ST_SetSRID($1::geometry,your_srid)
ST_Transform($1::geometry, 32468) ---
ST_Transform(ST_SetSRID($1::geometry,your_srid),
32468)
Cheers,
Rémi-C
2015-02-17 20:52 GMT+01:00 Miller, Stephan smill...@harris.com:
Remi –
I have it working with one exception: my embedded
ST_Transform($1::geometry, 32468) has stopped working. It is not
transforming lat/lon to a local
Re: [postgis-users] Getting TopologyExections when trying to node linestrings to create an overlay
Hey, last chance for you ^^ I successfully was able to remove errors in 2 ways : - first your geometry array contains a lot of duplicated geometry. ( rest of this is based on geom_set_id= 1) I suspect you forgot a where in an inner join, or something like this. Removing the duplicate in your data seems to solve the problem (removing duplicate : DISTINCT ON (geom)) I generated each pair of geometry in this set, none gave the error individually (which lessen the chance of a legit GEOS problem) - second, without removing the duplicates, a combination of tanslate and snaptogrid(0.5) was sufficient to also remove the error. I would recommend that you analyse a little bit your data. For instance, simply deuplicating on geom make your data set going from --15317 geometry in arrays to --9934 Snapping to a 0.5 grid before deduplicating further reduce the data set to 9617 (which might indicates that somewhere in your workflow you have a precision related issue). I never used ARCGIS, but I can bet you that ARCGIS is not precision-safe. The only product that is truly safe is CGAL, which comes with other type of constraints. You could also probably use GRASS safely if cleaning the data with v.clean on import. Anyway, there is no safe tool, only safe way to use it. You could very easily create another SRS that is a translation of your original srs to increase precision. So you simply change your workflow to ST_Transform before computing, then ST_Transform after computing. Lastly, I don't see the interest of using an array of geom here because you don't need to be able to do my_array[N] (no need to access). So you could simply use a geometry collection, so your input is a geom, and not a geom[]. Cheers, Rémi-C 2015-02-18 10:39 GMT+01:00 BladeOfLight16 bladeofligh...@gmail.com: On Mon, Feb 9, 2015 at 6:25 AM, Sandro Santilli s...@keybit.net wrote: First of all I confirm it still happens with GEOS=3.5.0dev-CAPI-1.9.0 r4038. Second, I took a look at a random set (geom_set_id=1) and I found it pretty big. That's to say you could probably further reduce the dataset for the ticket. That set contains 109 polygons, I can get the error by attempting to union the boundaries of the first 40 in that set, and I'm sure you can further reduce the input. Thanks for taking a look. I'll work on doing that when I can find the time, but I don't expect that to be a fast process at all. Even just checking for the pairwise case took a decent amount of time to develop the query and took overnight to finish running (even with optimizations like a bounding box intersection test). I don't really have any good heuristic that could narrow down the possibilies for reproducing, so I don't see much option other than to brute force it possibly with some kind of filter. That's why I didn't put more effort into shrinking the input set to begin with. Is a PostGIS database dump an okay format to provide the shapes, or would you prefer something else? I suppose I could dump groups into shapefiles or something like that if it's more convenient. On Mon, Feb 9, 2015 at 7:00 AM, Rémi Cura remi.c...@gmail.com wrote: this is a precision related issue, coordinates are way too big and should be translated. I understand what you mean by that (as in floating point problems due to size), but these coordinates are pretty typical. This is a standard UTM projection, zone 15N in middle America. It's even predefined in PostGIS' list of spatial references. I'm told ESRI had this kind of problem years ago, but they dealt with it as far as I know. While I would choose PostGIS over ESRI any day, this could be viewed by my coworkers as a good argument against using PostGIS; it represents a serious reliability concern since it applies to a broad range of functions and 2D projections. Basically, if you use any projection with coordinates of this size (of which there are a good number), there seems to be no telling when any function will just blow up in your face at random. On Mon, Feb 9, 2015 at 7:04 AM, Rémi Cura remi.c...@gmail.com wrote: Hey, I executed your data, the following command solve the problem (with very recent GEOS for me) (POSTGIS=2.2.0dev r12846 GEOS=3.5.0dev-CAPI-1.9.0 r0 PROJ=Rel. 4.8.0, 6 March 2012 GDAL=GDAL 2.0.0dev, released 2014/04/16 LIBXML=2.8.0 RASTER) [Snip] The change compared to your approach : convert input to table of simple polygons, (no array, no multi). Then translate to improve precision in geos computing Then the union. I don't really understand what you are trying to do, but ist_union seems dangerous and quit ineffective for that . I'm looking at this function call in your code: ST_Union( ST_MakePolygon(ST_ExteriorRing(geom)) ). That call seems to remove all holes and then create the union of all the covered areas (a single multipolygon that covers all areas covered by the originals). That is not what I'm trying to do; I already have an outer
Re: [postgis-users] Getting TopologyExections when trying to node linestrings to create an overlay
Computing finished successfuly on your data. The change I did is : ST_Polygonize(ST_Node(ST_Collect(geom)) ) It is awfully slow tough (1 hour for all you dataset I think) Cheers, Rémi-C 2015-02-18 13:07 GMT+01:00 Rémi Cura remi.c...@gmail.com: Hey, last chance for you ^^ I successfully was able to remove errors in 2 ways : - first your geometry array contains a lot of duplicated geometry. ( rest of this is based on geom_set_id= 1) I suspect you forgot a where in an inner join, or something like this. Removing the duplicate in your data seems to solve the problem (removing duplicate : DISTINCT ON (geom)) I generated each pair of geometry in this set, none gave the error individually (which lessen the chance of a legit GEOS problem) - second, without removing the duplicates, a combination of tanslate and snaptogrid(0.5) was sufficient to also remove the error. I would recommend that you analyse a little bit your data. For instance, simply deuplicating on geom make your data set going from --15317 geometry in arrays to --9934 Snapping to a 0.5 grid before deduplicating further reduce the data set to 9617 (which might indicates that somewhere in your workflow you have a precision related issue). I never used ARCGIS, but I can bet you that ARCGIS is not precision-safe. The only product that is truly safe is CGAL, which comes with other type of constraints. You could also probably use GRASS safely if cleaning the data with v.clean on import. Anyway, there is no safe tool, only safe way to use it. You could very easily create another SRS that is a translation of your original srs to increase precision. So you simply change your workflow to ST_Transform before computing, then ST_Transform after computing. Lastly, I don't see the interest of using an array of geom here because you don't need to be able to do my_array[N] (no need to access). So you could simply use a geometry collection, so your input is a geom, and not a geom[]. Cheers, Rémi-C 2015-02-18 10:39 GMT+01:00 BladeOfLight16 bladeofligh...@gmail.com: On Mon, Feb 9, 2015 at 6:25 AM, Sandro Santilli s...@keybit.net wrote: First of all I confirm it still happens with GEOS=3.5.0dev-CAPI-1.9.0 r4038. Second, I took a look at a random set (geom_set_id=1) and I found it pretty big. That's to say you could probably further reduce the dataset for the ticket. That set contains 109 polygons, I can get the error by attempting to union the boundaries of the first 40 in that set, and I'm sure you can further reduce the input. Thanks for taking a look. I'll work on doing that when I can find the time, but I don't expect that to be a fast process at all. Even just checking for the pairwise case took a decent amount of time to develop the query and took overnight to finish running (even with optimizations like a bounding box intersection test). I don't really have any good heuristic that could narrow down the possibilies for reproducing, so I don't see much option other than to brute force it possibly with some kind of filter. That's why I didn't put more effort into shrinking the input set to begin with. Is a PostGIS database dump an okay format to provide the shapes, or would you prefer something else? I suppose I could dump groups into shapefiles or something like that if it's more convenient. On Mon, Feb 9, 2015 at 7:00 AM, Rémi Cura remi.c...@gmail.com wrote: this is a precision related issue, coordinates are way too big and should be translated. I understand what you mean by that (as in floating point problems due to size), but these coordinates are pretty typical. This is a standard UTM projection, zone 15N in middle America. It's even predefined in PostGIS' list of spatial references. I'm told ESRI had this kind of problem years ago, but they dealt with it as far as I know. While I would choose PostGIS over ESRI any day, this could be viewed by my coworkers as a good argument against using PostGIS; it represents a serious reliability concern since it applies to a broad range of functions and 2D projections. Basically, if you use any projection with coordinates of this size (of which there are a good number), there seems to be no telling when any function will just blow up in your face at random. On Mon, Feb 9, 2015 at 7:04 AM, Rémi Cura remi.c...@gmail.com wrote: Hey, I executed your data, the following command solve the problem (with very recent GEOS for me) (POSTGIS=2.2.0dev r12846 GEOS=3.5.0dev-CAPI-1.9.0 r0 PROJ=Rel. 4.8.0, 6 March 2012 GDAL=GDAL 2.0.0dev, released 2014/04/16 LIBXML=2.8.0 RASTER) [Snip] The change compared to your approach : convert input to table of simple polygons, (no array, no multi). Then translate to improve precision in geos computing Then the union. I don't really understand what you are trying to do, but ist_union seems dangerous and quit ineffective for that . I'm looking at this function call in
[postgis-users] [OFF-TOPIC] SQL question
I have an interesting join problem that I'm not sure how to tackle in sql. I have two tables that I need to join (well really match) like this: create table a ( gid serial, item text, data text ); create table b ( gid serial, item text, moredata text ); gid can be used to force correct ordering of the records. item can be one or more rows with the same value in it (for example there might be 10 item='foo' records) data and more data is stuff related to the given record. What I know is that: 1. if table a has N rows for a given item then table b will have N rows 2. that the N rows in table a and table b are in the same order by gid 3. it is not safe to assume that a.gid=b.gid will link the correct records What I need is to match/join is: a.item.row[1] to b.item.row[1] a.item.row[2] to b.item.row[2] ... a.item.row[N] to b.item.row[N] where a.item=b.item Any thoughts on how to solve this with SQL? select aa.item, aa.cnt, bb.cnt from (select item, count(*) as cnt from a group by item) aa left outer join (select item, count(*) as cnt from b group by item) b on aa.item=bb.item order by aa.item; So likewise, we should be able to do something like: select aa.item, aa.data, bb.moredata from (select item, data from a order by item, gid) aa left outer join (select item, moredata from b order by item, gid) b on aa.item=bb.item and a.itemrow=b.itemrow order by aa.item; This needs some way of assigning itemrow numbers that can be matched. I think this is an application of over() but I'm not quite sure how to apply it. Thoughts would be appreciated :) Thanks, -Steve ___ postgis-users mailing list postgis-users@lists.osgeo.org http://lists.osgeo.org/cgi-bin/mailman/listinfo/postgis-users
