Geodesic Convolution Buffers

[mclaeysb]

I had some free time and dove into the sea of buffer ellipses, took on the winding number snake and sword-fought with edge-intersections. It was hard, but I endured, and today I can finally present to you the result of my quest:

This pull implements geodetic buffers for GeoJSON Points, LineSegments and Polygons (possibly having inner loops) and their Multi- counterparts, ... without JSTS dependency!
I included a basic test, but please test your own files and let me know if you encounter errors or problems.

Here's how I went about it:

• The CGAL documentation explains two different approaches for buffering general polygons (convex polygons are trivial): decomposition and convolution. It first seemed like a viable option to use the decomposition approach in combination with some triangulation algorithm (like earcut), but union-ing all constituent polygons is probably computationally too intensive (although cascading buffers could be used to speed up things). The convolution approach seemed a better idea.
• Basically, one should first compute the offset of a polygon (or lineSegment) by going around it right-handedly (and walk around convex angles right-handedly by going around the vertex), and then select only the area's with positive winding number. GeoJSON, however, doesn't know about winding numbers and only a handful of packages exist that use GeoJSON + winding numbers in some way.
• I looked in literature, and found an interesting master thesis entitled 'Partition of a Non-simple Polygon Into Simple Polygons' [1] describing an algorithm that decomposes a self-intersecting polygon into it's constituent simple polygons (+ yields their winding number). It was developed in the field of image rendering techniques, printers and page description languages, where the problem occurs of rendering a polygon using either the even-odd parity rule or the non-zero winding rule. It is (was) a novel object-precision algorithm in a world where discrete image-precision scan-line algorithms are most common.
• Based on the article I wrote the simplepolygon package. It allows for doing the necessary selection: give it an offset polygon, and it will return all constituent simple polygons with their parent simple polygon, winding number and net winding number.
• I also wrote the geojson-polygon-self-intersections package which computes the intersection vertices and their properties, and is used in simplepolygon.
• I ended up using a trick from the 'Polygon Offsetting by Computing Winding Numbers' [2] paper some people mentioned before: the convex angles can be shortcut since their other side is usually also part of the buffer polygon. This makes the offset polygon less heavy and prevents it of having non-unique vertices or zig-zagging double arcs in case of lineSegments.
• Finally, I (re-)wrote turf-buffer (I branched from [not ready] artisinal #4 'artisinal' but ended up only using it's point buffer). It reads the feature and preprocesses (rewinding polygons if necessary, and removing duplicate vertices since simplepolygon doesnt like them). It creates offset buffer polygons by making arcs around every vertex of a line, polygon, ... (shortcutted if convex) which use the turf-destination package to ensure geodesicness. It then passes the offset polygon to simplepolygon, and filters for netWinding = 1 (outer line) and netWinding = 0 (holes). It then unions all outer lines and holes, and subtracts them. The result is a geodesic buffer!

Known issues

• Polygons which aren't conform the Simple Feature standard can give errors. Most often, they are due to non-unique vertices (e.g. inner ander outer rings of polygons touching) which simplepolygon can't handle. Other cases include polygons or lineString with spikes (..., [0,0], [1,0], [1,1], [1,0], [2,0], ...) which also create offset polygons with non-unique vertices.
• Polygons with a very high amount of points and a buffer radius in the range of the feature-diameter can produce artifacts which could be due to floating point errors or an error with a euclidean (non-geodetic) algorithm such as 'turf-whithin' (note that this yields a super super complex offset polygon which very very high density of self-intersection points).

How geodesic is this approach?

• It uses turf-destination which uses the haversine to compute destinations etc.
• Issues arise when crossing the pole with a buffer, as with many polygon algorithms.
• simplepolygon is in se a euclidean algorithm, but should behave the same on a sphere. One important note is the computation of intersection points. If the input data is fine enough, this should be ok. Of not, we could use a geodesic line intersection algorithm. Ref. the note about the geodetic-ness of simplepolygons.

I really like geodesic-ness (I hope that's a word) discussions, and I hope to see more of them in Turf / live / ... discussions (since the Turf approach to not "project it + do a euclidean algorithm" but instead "do it on a sphere directly" is pretty fundamental to it) (Also also ref many discussion in turf-square-grid etc.). Let's talk!

Further development

• test & bench
• improve speed by optimising code / spatial index in simplepolygon / simplifying polygon prior to treatment if radius is large compared to inter-vertex-distance (turf-simplify?)
• Inside buffer for polygons

Hopefully soon...

Citations

[1] Subramaniam, Lavanya. Partition of a Non-simple Polygon Into Simple Polygons. Diss. University of South Alabama, 2003
[2] Chen, Xiaorui, and Sara McMains. "Polygon offsetting by computing winding numbers." ASME 2005 International Design Engineering Technical Conferences and Computers and Information in Engineering Conference. American Society of Mechanical Engineers, 2005.

[mclaeysb]

I just notices I should have better submitted by pull-request to the artisinal branch. Can I undo this? If you can (and agree), please do!

[tmcw]

Per the JavaScript norms currently, extending built-in objects can be dangerous, since the extensions may break other code or native behavior. Could this be implemented instead as a function that takes two arrays, and the same with the mod extension?

[mclaeysb]

Thanks for the feedback on my previous comment!

As you can see, I’ve continued to work on this turf-buffer PR and on its dependencies simplepolygon and geojson-polygon-self-intersections - albeit bit by bit, since this is spare time work for now (and it's lovely summertime outside and stuff).

I tested this implementation against all fixtures in ‘tile-cover’, as suggested by @morganherlocker. I found three minor issues, and was able to fix all three of them quickly. Now all 179 features from this test-set produce a correct buffer (I visually checked all of them). And look: the buffer radius nicely scales with latitude when projected in webmercator: :)

Next, I implemented a spatial index in simplepolygon’s two double for-loops, which are now one for-loop with one tree search. The resulting new benchmark measurements give:

	JSTS ops/sec	geodesicconvolution ops/sec	fastest	by
point	8,784	52,964	geodesicconvolution	6x
lineString	4,382	1,321	JSTS	3.3x
lineString large	249	21.19	JSTS	11.8x
polygon	6,455	2,346	JSTS	2.8x
polygon large	710	101	JSTS	7x

I’m very happy with these speed gains here (and loving rbush). There's no clear slowest step in turf-buffer or simplepolygon anymore. Next (and probably final) attempt to gain speed will be to turf-simplify large features. Comming up! Ideas/opinions welcome!

[mclaeysb]

Here we are!

I implemented line and polygon simplification. Super small edit, then quite some testing and debugging.

Remarkably, although turf-simplify doesn't guarantee to preserve topology when polygons have holes (and breaks it in some of the test cases), this doesn't prevent the calculation of a correct buffer. :)

What's the effect, you ask...

Performance

I set the tolerance of turf-simplify to be the bufferradius/20, which was the best trade-off between speed-gain and detail-loss. With this value, there is quasi no difference with the buffer of the non-simplified feature. As a result, there is a huge speed gain for features where the inter-coordinate-distance is much smaller then the buffer radius (for example when buffering super-detailed features, or when taking very large buffers). In those cases, this implementation beats JSTS in speed by up to 10x.

In other cases, such as for features with very few coordinates (which can't be simplified much more) or small buffer radii compared to the inter-coordinate-distance (such as in the tests I did earlier), there is quasi no speed gain. Therefore, the previous metrics still apply.

Known issues

• Features not conforming to the Simple Feature standard (see earlier)
• LineString that have very large (> 45°) jumps between successive points (inter-coordinate-distance) and are self-intersecting. This represents a new category of features I did find issues with. These are probably caused by the fact that simplepolygon and geojson-polygon-self-intersections are not perfectly geodetic (some sub-routines compute intersection points, within-ness, ... in euclidean space). I do think it’s fair to require a sufficiently oversampled input feature though.
  To fully solve this problem would include: 1) making a geodesic line intersection variant of geojson-polygon-self-intersections (not so hard), 2) making simplepolygon fully geodesic (computing angles locally and using bearings, ...) (quite some work) and 3) making union and difference geodesic (super hard).
  Polygons may not self-intersect according to the Simple Feature standard, and I (tried to but) didn’t find any polygons with this issue. That being said, it’s noteworthy to mention that the buffer of a polygon with large inter-coordinate-distance can appear to be wrong when projected on webmercator, but isn’t when displayed on a sphere (although it’s hard to find software that doesn’t resample these large steps - most often is wgs84 - before drawing on a sphere).
• There's one strange case I can't get my head around: the 'degenring.geojson' file from the tile-cover test fixture can be buffered succesfully (and quickly) with radius 1km and 100km, but 10 and 20 km fail. It appears simplepolygon makes a mistake when computing the 'parents'. I don't yet understand what goes wrong...

Furthermore, this commit fixes the second category of issues I mentioned in the first post.

And now...?

This pull request kindly yet impatiently awaits further testing, comments and/or approval. @tmcw @mourner @morganherlocker @tcql or anyone else: let's here it from you!

[morganherlocker]

Wow, amazing work!

Features not conforming to the Simple Feature standard (see earlier)

Are you saying this implementation fails when it receives non-simple features, or that it sometimes outputs non-simple features?

[mclaeysb]

I mean that this implementation fails on some input features that aren't conform the Simple Feature standard. I talked a bit about some cases in my first comment on this PR, and just wanted to mention it again here for completeness.

[aaronlidman]

Hey @mclaeysb, great meeting you in Brussels.

Super interesting work here. Are the known issues listed above something you're working on? What else does this branch need?

[mclaeysb]

Hey @aaronlidman, thanks for checking this out!

I'm not currently working on these issues, since most of them are OK:

• The first is solved by requiring Simple Feature, which seems reasonable to me and @morganherlocker.
• The second is solved by requiring sufficiently spatially oversampled input, which is an assumption made by almost any geo-tool.
• The third was looked into for quite some time by me, but I couldn't figure it out and it's the only case I had an issue with it.

In my opinion, this branch mostly needs:

• A more detailed time analysis. Either by having a decent profiler telling how much time is spent in which function (which could for example tell us much speed could be gained by replacing turf.destintion() by it's cheap-ruler counterpart). Or by manually adding a timer around each little function (but I'm not in for that + gets tricky with nesting).
• If there's more explicit need for this functionality, one could try to rewrite simplepolygon for speedup and better topological support.
• Testing in a full-scale workflow (tile-reduce / fancy features / ...)

Also, this is my first JS project and my first PR, so any comments on my code, or tips & helpful tools, are much appreciated!

In the meantime, I found out about @tmcw's mom's cookies and am hoping this is somewhat eligible...?

[DenisCarriere]

@mclaeysb Reference your profiling comment, the debug library can help you out with that, it can show you the time difference in milliseconds between the each debug call. To include it in the source code would be optional, but it's great when your developing, it acts like a more robust console.log.

https://www.npmjs.com/package/debug#millisecond-diff

[mclaeysb]

@DenisCarriere debug looks cool, I'll consider it for my next debugging session (just printing a tuned version of process.hrtime() while debugging now).

But more than knowing the time spent at each successive line in this (asynchronous) code, it'd be good to know, say, the total (summed) time spent in the many turf-destintion() or turf-destintion() or turf-helpers.point() calls, which are spread throughout the code. From what I understand from debug, it can't do that, or can it?

[DenisCarriere]

@mclaeysb I think what you're looking for goes beyond the scope of debug, sounds like your debugging solution will get really complex really fast. Best of luck! Let me know what you're final outcome was.

[morganherlocker]

@mclaeysb v8 has some decent profiling flags for this type of info that you can access from the node binary. I have never been able to found particularly good docs for them, so I'll demonstrate a simple case here.

example

Say we have a benchmark we want to run called bench.js that calls some complex algorithm code. Normally we would run this using node bench.js. Instead, use the --prof flag when calling node:

node --prof bench.js

This will run your benchmark in a special v8 debug mode that samples the v8 process at a set interval (default 1 ms), and log some info about what the process was working on each time it was sampled. You can set additional flags to log more or less data (ie: garbage collection events, memory usage, etc.). Once the process exits, you will see a .log file in your current directory containing a bunch of giberish-y looking info like this:

[...]
code-creation,LazyCompile,0,0xd17a6c4bf20,956,"tickDone node.js:299:22",0x1cc29dc0b370,~
code-creation,LazyCompile,0,0xd17a6c4c2e0,888,"emitPendingUnhandledRejections node.js:494:44",0x1cc29dc0bb50,~
tick,0x7fff9594dc22,323695,0,0x0,5
tick,0x7fff9594dc22,324694,0,0x0,5
code-creation,LazyCompile,0,0xd17a6c4c660,1556,"afterWrite net.js:773:20",0x363b3009dac0,~
code-creation,Stub,12,0xd17a6c4cd60,221,"CompareICStub"
code-creation,Handler,3,0xd17a6c4ce40,145,symbol("nonexistent_symbol" hash 4fd0399)
code-creation,Stub,12,0xd17a6c4cfc0,221,"CompareICStub"
code-creation,LazyCompile,0,0xd17a6c4d0a0,236,"WritableState.onwrite _stream_writable.js:88:26",0x363b300accd8,~
[...]

Node has another flag that allows for v8 to parse and aggregate these dumps, --prof-process. You can run the aggregation like this:

# note that the name of the log is generated on the fly; 
# change the name based on your actual log file
node --prof-process v8.log

This will output a helpful summary of how much time your process spent on each function, split out by naitive, C++ functions and JavaScript functions.

[JavaScript]:
   ticks  total  nonlib   name
     10    3.8%    4.2%  LazyCompile: *Voronoi.attachCircleEvent /Users/morgan/Documents/projects/islands/node_modules/voronoi/rhill-voronoi-core.js:1063:47
      9    3.4%    3.8%  LazyCompile: *Voronoi.leftBreakPoint /Users/morgan/Documents/projects/islands/node_modules/voronoi/rhill-voronoi-core.js:714:44
      4    1.5%    1.7%  LazyCompile: ~Voronoi.RBTree.rbInsertSuccessor /Users/morgan/Documents/projects/islands/node_modules/voronoi/rhill-voronoi-core.js:140:64
      3    1.1%    1.3%  LazyCompile: *Voronoi.RBTree.rbInsertSuccessor /Users/morgan/Documents/projects/islands/node_modules/voronoi/rhill-voronoi-core.js:140:64
      2    0.8%    0.8%  Stub: FastNewClosureStub

[...]

 [Summary]:
   ticks  total  nonlib   name
     64   24.3%   26.9%  JavaScript
    171   65.0%   71.8%  C++
     17    6.5%    7.1%  GC
     25    9.5%          Shared libraries
      3    1.1%          Unaccounted

[...]

I usually start at the top of this list until every function is as fast as I can possibly make it. Functions with an * in front of their name have been optimized by the v8 engine, so pay particularly close attention to functions without that notation. Most of the time, you'll find some obvious functions doing work that can be cached or approximated in some way.

PS: There are many ways to visualize these logs -- some built into v8 and some not. Here's some samples I often rely on, but generating them is probably more appropriate for an extended blog post. 😄

--log_timer_events with v8's profviz tool

flamegraph visualizing callstack depth

[mourner]

Another approach which should be simpler and more convenient in this case is running Chrome DevTools profiler in the browser (preferably on a very big polygon to get statistically significant data).

[flipper44]

I am so looking for this, I noticed Turf was sending out inaccurate Geodesic buffers so I dropped it.
The one thing I am curious about, was that it looks the Arc around the feature points is segmented as expected but doesn't visually seem to prioritize that the a segment end directly off the point(vertex), but looks as if creates segments from the start to the end..I was thinking it should Calculate the Center segment off the point, then calculate the segments of the arc to the right and left, forcing at least one between even if it is too small. This is hard to word so I can try to clarify.
Is there any chance to Have an option of how much segmentation per 360 degree arc like the circle tool has...so we could have smoother edges. I started playing around with the turf.destination tool and the circle tool and I loved how they worked even as they went past the poles.

[flipper44]

I notice it has been since September since the last comment, is there something holding this up? Would I be able to grab the non-Minified version and replace the buffercode with the code in index.js
from https://github.com/Turfjs/turf-buffer/pull/33/files to test this out. Or would this not work?

Edit: This doesn't look Easy to do after looking at the non-minified code. I could try to figure out how to build this on windows with these changes I guess. Any shortcuts on how I can research making a build with your changes?

[mclaeysb]

Hey @flipper44. Thanks for your intrest in this pull!

I don't understand you first question/remark regarding the 'center segment' etc. Could you try to explain it again?

Regarding your other remarks, there points might help you out:

• The resolutionparameter in the main function is there to let you tune the amount of segments in each arc.
• This pull is currently quite quiet because, as I explained in my post on Sep 27, 2016, we mainly need investigate how we can speed this code up, since it's still a bit slow as compared to JSTS in some cases. And the main reason such time-testing has not happend yet is because I have a full time job now.. :) I'll see if I can find some time soon to continue this work..
• The easiest way to test this out is to download all the files from this pull on my geodesicconvolution branch, run npm install to get the dependancies and then create a script that loads this code and runs it, e.g.:

var buffer = require('your/location/of/turf-buffer');
var fs = require('fs');
var input = fs.readFileSync('your/location/of/input.geojson','utf-8');
input = JSON.parse(input);
var output = buffer(input, 100, 'kilometers');

Cheers!
Manuel

[flipper44]

• My first comment was based on the arc segments I saw, but knowing that the resolution can be changed helps and makes that discussion Obsolete, it would be nice if it was an optional parameter however.
• I did see your comment about speed concerns, but thought since the current buffer attempts geodesic but does it inaccurate that any fix, even a slower fix, would be better than no fix. See image below, wouldn't it be better to get accurate Line, Polygon and Point buffers slower than Getting something incorrect fast?
• I may be missing something but version 3https://npmcdn.com/@turf/turf/turf.min.js attempts geodesic buffers and fails.

Circle vs Buffer Above

[mclaeysb]

Hi @flipper44

I'll answer your points here:

• OK!
• You could indeed argue that it's better to provide the right solution in more time than the wrong solution in a faster way, but a decision to include (the current state of) this pull request would be up to the turf maintainers. Since turf-buffer is most often used for small polygons, the error is not as case-dependent as with larger polygons - it is more systematic: the current turf-buffer outputs are not made by 2D-buffering a polygon 'in WebMercator projection' (which would result in a circle in WebMercator as in your image), but they are made by 2D-Buffering 'in WGS84 projection' (resulting in a ellipse in WebMercator).
• This pull refers to a branch on my fork of turf-buffer, and hence is not part of turf's npm upload. The branch 'artisinal' that you'll find was an earlier attempt of geodesic buffers by @morganherlocker and others.

Regards!

[flipper44]

In terms of speed, JSTS isn't using Haversine calculations and distances, there is no way that A geodesic buffer can even come close in speed to a planar buffer in speed. And the Current buffer only seems to do Haversine North and south and nothing in between, When I look at the sheer area of error I still am at a loss for why More GIS people aren't also pushing for this. Perhaps It would make sense to have three functions, one called PlanarBuffer which would be best for local coordinates systems, One called FastGeodesicBuffer, and one called TrueGeodesicBuffer (yours)

[flipper44]

Got past the warnings on npm install
When I try to load the script the requires error out.
var simplepolygon = require('simplepolygon');
var destination = require('@turf/destination');
var bearing = require('@turf/bearing');
var helpers = require('@turf/helpers');
var union = require('@turf/union');
var difference = require('@turf/difference');
var simplify = require('@turf/simplify');

is there a javascript library that is missing?

[tmcw]

When I try to load the script the requires error out.

General rule of thumb: if you see an error, specify it. What is the error? What is the exact text of the error message?

[flipper44]

I was missing the browserify command

…

On Feb 16, 2017 3:06 PM, "Tom MacWright" ***@***.***> wrote: When I try to load the script the requires error out. General rule of thumb: if you see an error, specify it. What is the error? What is the exact text of the error message? — You are receiving this because you were mentioned. Reply to this email directly, view it on GitHub <#33 (comment)>, or mute the thread <https://github.com/notifications/unsubscribe-auth/APQv9_em9SEsqJcJhFEoExQip9xODf3Qks5rdLpbgaJpZM4Iz2pa> .

[mclaeysb]

An update: Turf has been modularised, and turf-buffer has been replaced by @turf/buffer. I pushed my geodesicconvolutionbranch to master in my fork. Hope this makes it a bit more visible and easy to try out for those who want!

[ZacLiveEarth]

I'm very interested in seeing a geodetic buffer in turf. Did this discussion die off or move somewhere else?