Re: [postgis-users] Timezone for a given lat/long
I use select tz.tzid from location join tz_world tz on st_within(location."locationPoint", tz.geom) with the tz_world shapefile that I'm 99% sure I got from http://efele.net/maps/tz/world/ -- John Abraham j...@hbaspecto.com > On Mar 22, 2017, at 3:23 PM, suraj birla wrote: > > How to find the timezone for a given lat/long using postgis? > > We have to requirement to fnd the local timezone for millions of records > which has UTC timestamp and lat/long info. > > Google provide API but it's single request per call.. Not feasible for milion > records. > > Any suggestion? > > Thanks > Suraj > ___ > postgis-users mailing list > postgis-users@lists.osgeo.org > https://lists.osgeo.org/mailman/listinfo/postgis-users ___ postgis-users mailing list postgis-users@lists.osgeo.org https://lists.osgeo.org/mailman/listinfo/postgis-users
[postgis-users] A little freelance job to write a postgis function to cut a compass corner off a polygon
I need some functionality that I think is best put into a function, probably in pgsql. This has been on the todo list for years, but since I can't seem to get to it, maybe one of you would like to help out. I don't think it should take long for someone who has the skills. I'd really appreciate the help. You obviously need to have some experience writing postgis functions. The function would divide up a polygon by cutting off one of the compass corners using east-west and north-south lines. I think you'll have to use a bisection method to find the correct position of one of the lines, as the area (rather than the extent) of the corner-section to be cut off will be specified as a function parameter. I wrote down the planned algorithm and threw it on upwork here https://www.upwork.com/jobs/~01b6948ce256e51f85 <https://www.upwork.com/jobs/~01b6948ce256e51f85> (I'm not sure whether upwork is the best platform to ask for postgis freelance work, let me know if you have other ideas.) Hope someone can help, -- John Abraham j...@hbaspecto.com 403-232-1060 ___ postgis-users mailing list postgis-users@lists.osgeo.org http://lists.osgeo.org/mailman/listinfo/postgis-users
Re: [postgis-users] ST_Intersection very slow
Wow guys. 131 seconds is a lot faster than 10 days. Remi, if you can work out a fast intersection method that uses points and lines as a preprocessor to dealing with the polygons, I think it would be a great addition to PostGIS. ST_Intersection in PostGIS is often quite a bit slower than the implementation in Geomedia, Mapinfo, and (I hear) ArcGIS, so any functionality that results in speed improvements would be great. Mark, I can't wait to figure out why your system was fast! I was following your (preliminary) tutorial and gridding the data was progressing very slowly. I have a provincial boundary file but there seems to be much ambiguity in GIS representations of the provincial boundary, so I won't send you the one I have. I can try to assemble one from other sources. -- John Abraham j...@hbaspecto.com 403-232-1060 On Feb 25, 2015, at 6:28 AM, Mark Wynter wrote: > Hi John > > I’ve just crunched your whole dataset. The process takes 131 seconds for the > vector tiling (using a 16 CPU machine). Plus another 170 seconds for data > prep at the start including making the poly's valid. > For a 2 CPU machine, it will take circa 15 minutes, or double that using a > single CPU. > > Only one small issue outstanding - and that relates to clipping the regular > grid prior to tiling. For clipping I used the union of the > abmiw2wlcv_48tiles as supplied with the data - the problem is the > abmiw2wlcv_48tiles are rough and ready, which produces voids. The voids using > my method unfortunately get the same feature class as lc_class = 32. You’ll > see this clearly on second screenshot. > The way around this is to clip the regular grid using a high-res shapefile of > the Alberta state boundary prior to the tile crunching the > lancover_polygons_2010 table. This is easily done - I just didn’t have the > state boundary data. > > I need to get some sleep. John, Remi, I will share with you the code > tomorrow. For others, I’ll be posting a tutorial that steps through the > methods touched upon in this thread.. > > John, the only difference between my tutorial and its application to your > land cover data was a few tweaks to the data prep stage. Otherwise the same > code pattern (no modification at all needed to the worker_function). It was > great to test the code with your data. > > Speak tomorrow. > Mark > > > > > > > > > > > > ___ postgis-users mailing list postgis-users@lists.osgeo.org http://lists.osgeo.org/cgi-bin/mailman/listinfo/postgis-users
Re: [postgis-users] ST_Intersection very slow
Following advice from Mark and Remi, I am now tiling the troublesome layer. It's not going quickly. I'm trying an intersection of the layer with only 5 of my 10kmX10km grids and it has been running for 30 min. I have 8680 grids, which is more features than the layer I was originally working with (2709 polygons). I was hoping that intersection with squares would go much faster than an intersection with arbitrary polygons, but I can confirm that it's definitely not orders of magnitude faster. So, I think my idea of simplifying things recursively with a quad grid is worth a shot. Mark, can you point me to your Quadgrid function? Mark suggested counting points. Counting points (and rounding to the nearest 1) gives: select round(st_npoints(geom)/100,-2) as point100, count(*) from lancover_polygons_snap group by round(st_npoints(geom)/100,-2) order by point100 desc; 26300;1 1500;1 1100;1 800;1 600;2 500;1 400;5 300;10 200;22 100;141 0;997846 The one with about 263 points represents the developed land in the region (towns and cities including roads that connect them). Mark also suggested counting interior rings. select round(ST_NumInteriorRings(st_geometryN(geom,1)),-1) as rings10, avg(st_numpoints(geom)) as avgpts, count(*) as numfeatures from lancover_polygons_snap group by round(ST_NumInteriorRings(st_geometryN(geom,1)),-1) order by rings10 desc; 43140;2629905;1 4230;147121;1 3800;110083;1 2480;76515;1 1760;50261;1 1720;55576;1 1680;64588;1 1350;41826;1 1130;38059;1 970;36425;1 950;37686;1 860;32952;1 850;34657;1 820;35285;2 810;24769;1 800;28113;1 760;30726;1 740;26296;1 720;24788;2 670;18874;1 660;22428;1 640;26258;1 630;27213;2 620;22056;2 610;15771;1 560;19942;2 540;23988;1 500;15763;2 480;16598;1 470;19865;2 410;15466;1 400;16023;1 380;13536;4 370;14646;3 360;14122;1 340;12910;1 330;13273;1 320;15222;2 300;13421;1 290;9072;3 280;10025;4 270;12505;4 260;11002;5 250;9843;3 240;9980;6 230;11456;3 220;9049;3 210;9194;6 200;8586;5 190;7089;5 180;6796;8 170;6913;8 160;6955;6 150;5911;10 140;5769;15 130;5377;13 120;5301;17 110;5246;17 100;4422;28 90;4332;28 80;3664;35 70;3292;43 60;2860;66 50;2236;74 40;1882;141 30;1493;215 20;993;618 10;445;3325 0;20;993265 -- John Abraham j...@hbaspecto.com 403-232-1060 "POSTGIS="2.0.3 r11128" GEOS="3.3.8-CAPI-1.7.8" PROJ="Rel. 4.8.0, 6 March 2012" GDAL="GDAL 1.9.2, released 2012/10/08" LIBXML="2.7.8" LIBJSON="UNKNOWN" (core procs from "2.0.3 r11132" need upgrade) TOPOLOGY (topology procs from "2.0.3 r11132" need upgrade) RASTER (raster procs from "2.0.3 r11132" need upgrade)" On Feb 24, 2015, at 3:31 PM, Rémi Cura wrote: > Hey, > I highjack this interesting discussion, sorry =) > > If you have one single massive Polygon (no multi). > You can generate a point table (keeping path of course) with ST_DumpPoints. > > Then you construct an index on this table. > Then you perform your Intersection with indexed points, which will be very > fast. > Then you reconstruct parts of polygon. > It will require some advanced to very advanced SQL (like windows function) to > make it work tough. > > I think it would be faster than other alternatives, because basicaly it > amounts to using the optimized R-Tree index on points, vs a custom and > non-otpimized Quad tree index on polygon. > Moreover when using a recusrive strategy, you end up computing the same think > a lot with geos. > > Cheers, > Rémi-C > > 2015-02-24 23:10 GMT+01:00 Mark Wynter : > > > Thanks Mark. I will indeed do the st_dump, but I don't think it will help > > much (very few multis). I think the tiling will help a lot. What I > > wonder, though, is how long it will take to tile? Afterall, it's still an > > st_intersection operation to tile each geometry against each tile. > > I’ve almost finished writing the tutorial - where I address many of these > points.The variables that affect performance are: > * how you’ve written your ST_Intersection query > * multi vs. non-multi > * size and complexity of geoms > * no. available CPUs (for parallelisation) > * tile batch size - important!!! > > All strategies in combination may be necessary if your queries are taking > forever. > > For the demonstration dataset (a multi polygon representing whole of > Australia), my tutorial tiling query incorporates ST_Intersection and > ST_Difference subqueries to produce tiled features representing land and > water. > I achieved a 49x reduction in the query time to tile the whole of Australia, > starting with a single multi polygon. > > The more complex the query, the more significant this time saving is in > absolute terms. > > > Is there a quad tree type tiling algorithm i
Re: [postgis-users] ST_Intersection very slow.
Thanks Mark. I will indeed do the st_dump, but I don't think it will help much (very few multis). I think the tiling will help a lot. What I wonder, though, is how long it will take to tile? Afterall, it's still an st_intersection operation to tile each geometry against each tile. Is there a quad tree type tiling algorithm in a function? If I do 256 x 256 tiles, doing it all at once would be 65536 operations of st_intersection(complex_geom, square). With a quad tree I'll only have 4 operations of st_intersection(complex_geom, square), and then 16 of (a_little_less_complex_geom, square), and then 64 of (even_less_complex_geom, square) and so on, for 8 nests. The last one will still be 65536 operations, but the complex geoms should be a lot simpler by the time I get down to that level. What do you think, is this worth trying? Or is intersecting with a square fairly simple regardless of how complex the other geometry is? (Your reply initially fell into my junk mail, so I missed it, and also haven't tried these things yet, partly because I missed your email.) -- John Abraham j...@hbaspecto.com 403-232-1060 On Feb 19, 2015, at 6:52 PM, Mark Wynter wrote: > >> So I've was running this query for 866000 s (10 days) before I decided to >> kill it: > > >> One potential thing I've realized is that a few of the geometries in tazjan2 >> are multipolygons, not single polygons. But it's only a few. There are a >> few very large and complex polygons in lancover_polygons_snap, but again, >> most of the 998031 are simple small polygons, about half would be >> ST_CoveredBy the polygons in tazjan2 and most of the rest would only overlap >> two or three of the polygons in tazjan2. > > > A common offender - the multipolygons. > > You can get an immediate performance improvement by > ST_Dump(wkb_geometry).geom into a new table, with a new serial id, but > retaining the original poly_id so you can don’t lose the multipolygon > provenance of each polygon. > Complex operations on the multipolygon table, like ST_Intersection(), are > extremely slow - Because the query has to access and compute the operation on > each of the geometry elements contained within the multi. > > Another massive performance gain can be achieved by “Tiling” your dumped > polygon table - which has the effect of breaking your geometries into smaller > features, which you can then “union” back together again, after running > ST_Intersection on your vector tiles. Don’t try to tile your multi polygon > table - you’ll still be waiting days. > > Take the coastline of Australia - a single multi polygon > > SELECT Count(ogc_fid) as num_rows, SUM(ST_NumGeometries(wkb_geometry)) as > num_geoms, SUM(ST_NPoints(wkb_geometry)) as num_points FROM abs_aus11; > num_rows | num_geoms | num_points > ---+-- > 1 | 6718 |1622598 > > If you need to query the intersection of a polygon with land, and the water, > its near impossible. Many, many days. > > Hence we need to take a different approach. > > The first thing you can do is to dump the multi polygon table. > Then create a regular vector grid. > Then tile the both your Polygon tables using the vector grid. > Put indexes on the resultant tiled tables. > Then run your ST_Intersects. > The result is many more geometries, and points, but now we get the benefit of > geometry indexes. > > As a result of tiling Australia coastline, we now get > > SELECT Count(tid) as num_rows, SUM(ST_NumGeometries(wkb_geometry)) as > num_geoms, SUM(ST_NPoints(wkb_geometry)) as num_points FROM abs_aus11_tiled; > num_rows | num_geoms | num_points > ---+-- > 17222 | 29020 |6513770 > > Each geometry also has a tile_id, and the original poly_id. The tile_id’s > are really useful for subsetting your queries, and for using tile-id’s as an > additional join condition. At any time you can rapidly rebuild your original > geometries doing ST_Union() GROUP BY poly_id. > > Using the approach of dumping and tiling, queries that once took days now > takes minutes at most. And as Remi mentioned, you can parallelise the > process. The concept of vector tiling in PostGIS is analogous to Map Reduce > and assigning Key Value Pairs. Its about chunking your data, doing lots of > fast operations on the smaller bits, and then consolidating your outputs. > You don’t necessarily have to write a SQL function to do the SQL > parallelisation. My preference is Bash and GNU Parallel because you can > write vanilla SQL and execute via psql in parallel. For me its now about > the speed I can write the SQL. > > So we’re
Re: [postgis-users] ST_Intersection very slow.
Thanks for the hint, Rémi. I'll give it a try, to see if the st_intersection is slow (as compared to just st_intersects). Here's my EXPLAIN ANALYZE for my 5 zone test: "Nested Loop (cost=1.56..80.91 rows=8 width=560) (actual time=3929258.738..19756123.479 rows=29 loops=1)" " Join Filter: (_st_intersects(p.geom, n.the_geom) AND ((NOT (p.geom && n.the_geom)) OR (NOT _st_touches(p.geom, n.the_geom" " CTE sometaz" "-> Limit (cost=0.00..1.56 rows=5 width=9537) (actual time=0.433..0.576 rows=5 loops=1)" " -> Seq Scan on tazjan2 (cost=0.00..843.09 rows=2709 width=9537) (actual time=0.415..0.511 rows=5 loops=1)" " -> CTE Scan on sometaz n (cost=0.00..0.10 rows=5 width=80) (actual time=0.460..0.724 rows=5 loops=1)" " -> Index Scan using lancover_polygons_snap_geom_idx on lancover_polygons_snap p (cost=0.00..12.47 rows=5 width=480) (actual time=0.242..65.640 rows=16 loops=5)" "Index Cond: (geom && n.the_geom)" "Total runtime: 19756143.935 ms" Here it is explained in layspeak: http://explain.depesz.com/s/dbVi So, am I reading this right, that it took 5.5 hours to st_intersects 5 (multi)polygons with 29 (multi)polygons, and find the 8 that intersected? That the index scan was fast (328ms), and returned only a few (16, or maybe 29) of the 998031 polygons, but that the actual intersection afterwards, between just a few polygons, was so very slow? (The query was this: with sometaz as (select * from tazjan2 limit 5) SELECT p.lc_class, n.taz , CASE WHEN ST_CoveredBy(p.geom, n.the_geom) THEN p.geom ELSE ST_Multi( ST_Intersection(p.geom,n.geom) ) END AS geom FROM lancover_polygons_snap AS p INNER JOIN sometaz AS n ON (ST_Intersects(p.geom, n.the_geom) AND NOT ST_Touches(p.geom, n.the_geom) ); ) -- John Abraham j...@hbaspecto.com 403-232-1060 On Feb 19, 2015, at 11:56 AM, Rémi Cura wrote: > Hey, > you could try to not use CASE (so separate the spatial join from the > processing, which is easy to parallelize (assuming you have more than one > core )). > > First generate the table with > > CREATE TABLE psatial_mapping_between_lancover_and_taz AS > SELECT row_number() over() as row_id, p.your_primary_key, n.your_primary_key > FROM lancover_polygons_snap AS p >INNER JOIN tazjan2 AS n > ON (ST_Intersects(p.geom, n.the_geom) > AND NOT ST_Touches(p.geom, n.the_geom) ); > > > Then when you have your table, you can do your processing (put into a SQL > function , it is cleaner), > and you parallelize on the "row_id". > > Basically you get min(row_id), ax(row_id), you spearate it into N parts, then > execute the parts with your K processes. > (both bash and python have utilities for this). > > Of course if it is actually the spatial join that is slow, you can also > parallelise this. > (by cutting lancoverinto N pieces for instance ) > Cheers, > Rémi-C > > 2015-02-19 19:45 GMT+01:00 John Abraham : > So I've was running this query for 866000 s (10 days) before I decided to > kill it: > > create table taz_and_lancover_10_fast_2 as > SELECT p.lc_class, n.taz > , CASE >WHEN ST_CoveredBy(p.geom, n.the_geom) >THEN p.geom >ELSE > ST_Multi( > ST_Intersection(p.geom,n.geom) > ) END AS geom > FROM lancover_polygons_snap AS p >INNER JOIN tazjan2 AS n > ON (ST_Intersects(p.geom, n.the_geom) > AND NOT ST_Touches(p.geom, n.the_geom) ); > > Explain shows this, it's using the spatial index: > > "Nested Loop (cost=0.00..310492.35 rows=483105 width=9973)" > " Join Filter: (_st_intersects(p.geom, n.the_geom) AND ((NOT (p.geom && > n.the_geom)) OR (NOT _st_touches(p.geom, n.the_geom" > " -> Seq Scan on tazjan2 n (cost=0.00..843.09 rows=2709 width=9493)" > " -> Index Scan using lancover_polygons_snap_geom_idx on > lancover_polygons_snap p (cost=0.00..21.67 rows=5 width=480)" > "Index Cond: (geom && n.the_geom)" > > There are 2709 rows in tazjan2 and 998031 rows in lancover_polygons_snap, so > I appreciate that it's a bit of a large problem. But MapInfo was able to do > it interactively in a few days and Geomedia was also able to do it in about a > day. > > Both MapInfo and Geomedia ran out of memory (8GB machines) until the problem > was broken into two regions (North and South), but postgresql seems to be > chewing on the problem using only 400MB. The interactive approach in MapInfo > was to divide lancover_polygons_snap by lc_class, to further divide each > region into about 10 subproblems. Perha
[postgis-users] ST_Intersection very slow.
So I've was running this query for 866000 s (10 days) before I decided to kill it: create table taz_and_lancover_10_fast_2 as SELECT p.lc_class, n.taz , CASE WHEN ST_CoveredBy(p.geom, n.the_geom) THEN p.geom ELSE ST_Multi( ST_Intersection(p.geom,n.geom) ) END AS geom FROM lancover_polygons_snap AS p INNER JOIN tazjan2 AS n ON (ST_Intersects(p.geom, n.the_geom) AND NOT ST_Touches(p.geom, n.the_geom) ); Explain shows this, it's using the spatial index: "Nested Loop (cost=0.00..310492.35 rows=483105 width=9973)" " Join Filter: (_st_intersects(p.geom, n.the_geom) AND ((NOT (p.geom && n.the_geom)) OR (NOT _st_touches(p.geom, n.the_geom" " -> Seq Scan on tazjan2 n (cost=0.00..843.09 rows=2709 width=9493)" " -> Index Scan using lancover_polygons_snap_geom_idx on lancover_polygons_snap p (cost=0.00..21.67 rows=5 width=480)" "Index Cond: (geom && n.the_geom)" There are 2709 rows in tazjan2 and 998031 rows in lancover_polygons_snap, so I appreciate that it's a bit of a large problem. But MapInfo was able to do it interactively in a few days and Geomedia was also able to do it in about a day. Both MapInfo and Geomedia ran out of memory (8GB machines) until the problem was broken into two regions (North and South), but postgresql seems to be chewing on the problem using only 400MB. The interactive approach in MapInfo was to divide lancover_polygons_snap by lc_class, to further divide each region into about 10 subproblems. Perhaps subprobleming this is the way to go? Can't the query subproblem it based in the indices or would I have to do that manually? One potential thing I've realized is that a few of the geometries in tazjan2 are multipolygons, not single polygons. But it's only a few. There are a few very large and complex polygons in lancover_polygons_snap, but again, most of the 998031 are simple small polygons, about half would be ST_CoveredBy the polygons in tazjan2 and most of the rest would only overlap two or three of the polygons in tazjan2. I must be doing something wrong. Any hints? I have max_connections set to 100 (currently only 11 connections active) work_mem was defaulting to 1MB, I just bumped it to 256MB (machine has 32GB but it's Postgresql9.1 (x86), i.e. 32bit version shared_buffers was 32MB, I am trying 512MB, could go higher effective_cache_size was defaulting to 128MB, I am trying 20GB. random_page_cost was 4.0, this is VMWare virtual server. I am trying 2.0 I'm trying an "explain analyze" with just 5 rows of tazjan2, a 5 x 998031 intersection problem with indices shouldn't take that long, should it? It's been running for a little over an hour. I'm wondering if the few complex polygons in lancover_polygons_snap are causing the slowness? Is there some fast way to divide complex polygons, e.g. apply a 1km grid over them? PostGISfullversion(): "POSTGIS="2.0.3 r11128" GEOS="3.3.8-CAPI-1.7.8" PROJ="Rel. 4.8.0, 6 March 2012" GDAL="GDAL 1.9.2, released 2012/10/08" LIBXML="2.7.8" LIBJSON="UNKNOWN" (core procs from "2.0.3 r11132" need upgrade) TOPOLOGY (topology procs from "2.0.3 r11132" need upgrade) RASTER (raster procs from "2.0.3 r11132" need upgrade)" -- John Abraham ___ postgis-users mailing list postgis-users@lists.osgeo.org http://lists.osgeo.org/cgi-bin/mailman/listinfo/postgis-users
Re: [postgis-users] Getting TopologyExections when trying to node linestrings to create an overlay
Well I could say that using PostGIS ST_Intersects with "messy data" always seems to give me TopologyExceptions. I've had luck with various combinations of ST_SnapToGrid and ST_Buffer(0), but with messy data there always seems to be some weird case that requires manually edits. Sorry to be the bearer of bad news. Frankly, I don't understand why the GEOS library has to throw that error. I would encourage you to isolate particular problems and file bug reports. Improvements in GEOS to eliminate the underlying error(s) would certainly be welcome. -- John Abraham > On Feb 8, 2015, at 1:43 PM, BladeOfLight16 wrote: > > On Tue, Feb 3, 2015 at 11:28 AM, BladeOfLight16 <mailto:bladeofligh...@gmail.com>> wrote: > I'm trying to create a polygon overlay. The basic process is relatively > simple: 1) Get the boundaries 2) Union the boundaries to node the linestrings > 3) Polygonize the noded outlines 4) Filter out holes using a contains or > intersects test. The problem I'm running into is that I'm getting > "GEOSUnaryUnion: TopologyException: found non-noded intersection between > LINESTRING" errors from ST_Union. > > [a lot snipped] > > I haven't seen any response to this. I was just wondering if anyone else had > a chance or intentions to look this over. Granted, it's pretty long and > involved (sorry for that), but I thought all the details I included were > important. I do know it went through the mailing list; someone on IRC helped > me find it in the... I guess it's not the archives; I don't know what it's > called. But the online browsing mechanism. Thanks to anyone who's taking a > look. > ___ > postgis-users mailing list > postgis-users@lists.osgeo.org > http://lists.osgeo.org/cgi-bin/mailman/listinfo/postgis-users ___ 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)
Just FYI in case you’re interested I found it much easier to use the triangle algorithm here https://github.com/tudelft3d/pprepair on an exported shape file It’s called Planar Partition Repair. It got rid of all the slivers and overlaps with ease. Here is the paper http://www.gdmc.nl/ken/files/12_pfg.pdf and thanks to Martijn Meijers and team. Then I rebuilt the topology layer in PostGIS from the fixed up Simple Feature Class multi polygon layer. -- John Abraham j...@hbaspecto.com 403-232-1060 On Nov 3, 2014, at 9:13 AM, Rémi Cura wrote: > Just for completeness sake, > if your input geom is really ambiguous with the precision you choose, > it may be faster in both human and machine time to use grass 7 to import and > clean the topology, > then export from grass 7 to postgis topology with the built-in grass function. > > I ressort sometime to this trick because grass cleaning functionality are > very advanced. > > Cheers, > Rémi-C > > > 2014-11-03 17:04 GMT+01:00 Guillaume Drolet : > 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 > > ___ > postgis-users mailing list > postgis-users@lists.osgeo.org > http://lists.osgeo.org/cgi-bin/mailman/listinfo/postgis-users ___ postgis-users mailing list postgis-users@lists.osgeo.org http://lists.osgeo.org/cgi-bin/mailman/listinfo/postgis-users
Re: [postgis-users] postgis-users Digest, Vol 152, Issue 3
Yes that was basically my algorithm: combinations of buffer, intersection, then difference. I was using st_makevalid too, it didn't help, everything was considered valid, yet I was getting errors. I was unaware of st_snap as a separate function than st_snaptogrid. Maybe that was the missing piece. Thanks. Nonetheless, the triangulation algorithm from Martijn seems so much cleaner and more powerful than any of this. http://tudelft-gist.github.io/pprepair/ . I don't even see a reason to try st_snap. -- John Abraham Sent from my iPhone, please excuse any typos. > On Oct 4, 2014, at 6:18 AM, Mark Wynter wrote: > > Don't have the benefit of seeing your image... > For those that intersect, what about ST_difference to rid overlaps and > ST_Snap to fill voids... > > Or ST_Buffer both, take the intersection, union it to one of the polys and > difference it with the first original > > I recall a project 6 months back where client handed me a bunch of polys that > lacked common edges because they were created almost freehand ... > > I can't pull up my code because I'm on annual leave, but I recall I wrote a > stored procedure that looped through each poly at a time, found the nearest, > and used a combination of functions to produce a valid edge between each pair. > > HTH > > Sent from my iPhone > >> On 4 Oct 2014, at 4:30 am, postgis-users-requ...@lists.osgeo.org wrote: >> >> 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. A little tired of non-noded intersections betweenpolygons -- >>node them all (John Abraham) >> >> >> -- >> >> Message: 1 >> Date: Fri, 3 Oct 2014 12:29:30 -0600 >> From: John Abraham >> To: postgis-users@lists.osgeo.org >> Subject: [postgis-users] A little tired of non-noded intersections >> betweenpolygons -- node them all >> Message-ID: <5471efe5-b80b-45e5-8154-b8c3f691b...@hbaspecto.com> >> Content-Type: text/plain; charset="windows-1252" >> >> I?m trying to clean up a geometry layer, a division of a province into >> smaller areas. Once it?s cleaned up I?d like to start using a topology, so >> that no-one else ever has to deal with these problems? but? meanwhile... >> >> I keep getting non-noded intersections or other topology exceptions when I?m >> searching for (or removing) overlapping polygons or searching for (or adding >> in) missing bits between polygons. >> >> I?ve used all kinds of combinations of st_buffer(0), st_cleangeometry(), >> custom versions of st_cleangeometry(), and st_snaptogrid to try to fix >> things. >> >> I?ve noticed a pattern with the types of problems I?m currently dealing >> with. Basically, they are non-noded intersections between polygon edges. >> Eg, there is a polygon that has a not-quite-exactly-north-south boundary on >> its west side that it shares with its western neighbor. But the neighbor's >> boundary is longer, extending north past the corner of the original, and >> it?s defined by a different pair of points. If the image comes through, you >> can see it below. The Greenish polygon obviously has node at its north-west >> corner, but the brownish polygon, to the west of it, does not have a node at >> that location, its boundary is longer so is defined by more distant end >> points. >> >> I?m wondering if I could use something like st_node to add a node in the >> western polygon, so the corner between them is defined identically on both >> sides? Then, presumably, I wouldn?t be getting these >> ERROR: GEOSUnaryUnion: TopologyException: found non-noded intersection >> between LINESTRING (144297 5.66546e+006, 144201 5.67011e+006) and LINESTRING >> (144203 5.67002e+006, 144204 5.67002e+006) at 144203.19006961776 >> 5670017.7300232062 >> type of errors. But ST_Node only works on lines, not polygons. >> >> >> >> >> >> >> PS I was wat
Re: [postgis-users] A little tired of non-noded intersections between polygons -- node them all
Wow Martijn. Thank you SO much. It seems to have worked. This will get a lot of use, I think. You really should get it into PostGIS as a function, for better visibility and for the benefit of the world. I can’t wait to try it on one of my larger problems (e.g. a cadastral data set with about 2 million parcels). For everyone else, here’s the paper: http://www.gdmc.nl/ken/files/12_pfg.pdf -- John Abraham j...@hbaspecto.com 403-232-1060 On Oct 3, 2014, at 1:23 PM, Martijn Meijers wrote: > Hi John, > > We know the pain you describe. > > For this we've developed at our research group: > > * https://github.com/tudelft-gist/pprepair > > Documentation here: > > * http://tudelft-gist.github.io/pprepair/ > > Not integrated with PostGIS at the moment and only reading shapefiles, but it > should clean exactly the type of data that you describe into properly noded > polygons (with gaps and overlaps removed/filled). > > > Martijn > > -- > Martijn Meijers > > mailto:b.m.meij...@tudelft.nl > http://www.gdmc.nl/martijn > > GIS-technology > Faculty of Architecture and the Built Environment > Delft University of Technology > > P.O. Box 5030 | Julianalaan 134 (building 8), room 01.oost.330 > 2600 GA Delft | 2628 BL Delft > The Netherlands > > tel (+31) 15 2785642 > ___ postgis-users mailing list postgis-users@lists.osgeo.org http://lists.osgeo.org/cgi-bin/mailman/listinfo/postgis-users
[postgis-users] A little tired of non-noded intersections between polygons -- node them all
I’m trying to clean up a geometry layer, a division of a province into smaller areas. Once it’s cleaned up I’d like to start using a topology, so that no-one else ever has to deal with these problems… but… meanwhile... I keep getting non-noded intersections or other topology exceptions when I’m searching for (or removing) overlapping polygons or searching for (or adding in) missing bits between polygons. I’ve used all kinds of combinations of st_buffer(0), st_cleangeometry(), custom versions of st_cleangeometry(), and st_snaptogrid to try to fix things. I’ve noticed a pattern with the types of problems I’m currently dealing with. Basically, they are non-noded intersections between polygon edges. Eg, there is a polygon that has a not-quite-exactly-north-south boundary on its west side that it shares with its western neighbor. But the neighbor's boundary is longer, extending north past the corner of the original, and it’s defined by a different pair of points. If the image comes through, you can see it below. The Greenish polygon obviously has node at its north-west corner, but the brownish polygon, to the west of it, does not have a node at that location, its boundary is longer so is defined by more distant end points. I’m wondering if I could use something like st_node to add a node in the western polygon, so the corner between them is defined identically on both sides? Then, presumably, I wouldn’t be getting these ERROR: GEOSUnaryUnion: TopologyException: found non-noded intersection between LINESTRING (144297 5.66546e+006, 144201 5.67011e+006) and LINESTRING (144203 5.67002e+006, 144204 5.67002e+006) at 144203.19006961776 5670017.7300232062 type of errors. But ST_Node only works on lines, not polygons. PS I was watching Paul Ramsey’s presentation on the web from Foss4G and I agree that PostGIS is fantastic and gaining so much momentum, and perhaps even “asymptotically approaching perfection”. But these types of topology errors always seem to trip me up. Isn’t there something that can be done in the underlying GEOS functions to make these less common? Or is everything an edge case that needs a special algorithm? I’m frankly rather surprised that these errors don’t seem common for more people (based on my limited google search hits): are others not using ST_Union or ST_Intersects very much? Or does everyone else have boundaries that are defined in topologies, or with every point specified? Or is there just some trick, that no-one’s ever told me about, to add in the necessary points on edges so that identical line segments are defined by the same pairs of points? Feeling frustrated, -- John Abraham j...@hbaspecto.com 403-232-1060 ___ postgis-users mailing list postgis-users@lists.osgeo.org http://lists.osgeo.org/cgi-bin/mailman/listinfo/postgis-users
Re: [postgis-users] postgis-users Digest, Vol 143, Issue 15
Hey Paul there are many problems with the clean geometry functions, so this contribution is very useful. I'm trying it now, I notice it returns a GeometryCollection sometimes, so I had to use ST_CollectionExtract(cleangeometry(geom),3) to get my multipolygons back. The old function wasn't returning "null", but it was returning geometries where st_area() = null. Whatever that means. Yours seems to have fixed the problem in general, thank you. However comparing the ST_Area of the original and the cleaned, I see a few of my multi polygons have lost considerable area. In one case, a tiny infinitesimally narrow "finger" stuck out from the polygon, and a tiny section of the finger was retained but the bulk of the polygon was lost. In another case, it looks like the corner was twisted around itself, and so the twisted corner was retained and the bulk of the polygon was lost. I can manually make these edits, I think, and your algorithm is better than the original, but it's not perfect yet :) -- John Abraham j...@hbaspecto.com On Thu, Jan 16, 2014 at 1:00 PM, wrote:Message: 3 Date: Thu, 16 Jan 2014 14:42:29 +1000 From: Paul Pfeiffer To: postgis-users@lists.osgeo.org Subject: [postgis-users] Updated cleangeometry function Message-ID: Content-Type: text/plain; charset="iso-8859-1" I have been having problems with the clean geometry function referenced at http://trac.osgeo.org/postgis/wiki/UsersWikiCleanPolygons returning null geometries. So I have fixed the script and I'm providing it below if anyone else wants to test it and let me know I could update the wiki. (Requires Postgres v9 +) -- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -- -- $Id: cleanGeometry.sql 2014-01-16 Paul Pfeiffer -- -- cleanGeometry - remove self- and ring-selfintersections from -- input Polygon geometries -- -- Copyright 2014 Paul Pfeiffer -- Version 2.0 -- contact: nightdrift at gmail dot com -- -- modified from cleanGeometry.sql 2008-04-24 from http://www.kappasys.ch -- -- This is free software; you can redistribute and/or modify it under -- the terms of the GNU General Public Licence. See the COPYING file. -- This software is without any warrenty and you use it at your own risk -- -- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - CREATE OR REPLACE FUNCTION cleangeometry(geom "public"."geometry") RETURNS "public"."geometry" AS $BODY$ DECLARE inGeom ALIAS for $1; outGeom geometry; tmpLinestring geometry; sqlString text; BEGIN outGeom := NULL; -- Clean Polygons -- IF (ST_GeometryType(inGeom) = 'ST_Polygon' OR ST_GeometryType(inGeom) = 'ST_MultiPolygon') THEN -- Check if it needs fixing IF NOT ST_IsValid(inGeom) THEN sqlString := ' -- separate multipolygon into 1 polygon per row WITH split_multi (geom, poly) AS ( SELECT (ST_Dump($1)).geom, (ST_Dump($1)).path[1] -- polygon number ), -- break each polygon into linestrings split_line (geom, poly, line) AS ( SELECT ST_Boundary((ST_DumpRings(geom)).geom), poly, (ST_DumpRings(geom)).path[1] -- line number FROM split_multi ), -- get the linestrings that make up the exterior of each polygon line_exterior (geom, poly) AS ( SELECT geom, poly FROM split_line WHERE line = 0 ), -- get an array of all the linestrings that make up the interior of each polygon line_interior (geom, poly) AS ( SELECT array_agg(geom ORDER BY line), poly FROM split_line WHERE line > 0 GROUP BY poly ), -- use MakePolygon to rebuild the polygons poly_geom (geom, poly) AS ( SELECT CASE WHEN line_interior.geom IS NULL THEN ST_Buffer(ST_MakePolygon(line_exterior.geom), 0) ELSE ST_Buffer(ST_MakePolygon(line_exterior.geom, line_interior.geom), 0) END, line_exterior.poly FROM line_exterior LEFT JOIN line_interior USING (poly) ) '; IF (ST_GeometryType(inGeom) = 'ST_Polygon') THEN sqlString := sqlString || ' SELECT geom FROM poly_geom
