Take the 2-minute tour ×
Geographic Information Systems Stack Exchange is a question and answer site for cartographers, geographers and GIS professionals. It's 100% free, no registration required.

I am looking for the Math and steps involved not using a tool. We have PostGIS but its an old version so we are working with the data in Python.

Given three coordinates:

  • -107.37, 52.46
  • -103.24, 52.62
  • -101.09, 53.97

I then join them to form a line.

I then build buffered polygons around each line segment (20NM) - resulting in two rectangles. If there was another coordinate I would have 3 rectangles and so on.

This is what the output looks like when the rectangles are assembled:

Joined rectangles

However, what I really want is to have that triangular gap - sealed up with a curve.

Like this (did it by hand):

with curve

So I have everything working just fine and dandy, but I have no idea how to properly join these rectangles so that they have a nicely fitted curve to fill in that triangular gap.

Any help would be greatly appreciated.

Edit:

Since I am assembling the polygon data myself I can easily identify the points that need to be connected. Is there a way to join two points with a curve?

Thanks!

share|improve this question
    
The proper name for this exercise is "polygon simplification". It looks like you want a variant of the "bend simplify" algortihm. Identifying the interior bend is easy, but smoothing over the gap is likely to be more difficult. –  Vince Nov 7 '13 at 20:28
3  
@Vince It is hard to conceive of this operation as simplifying a polygon, given it begins with a polyline! Even starting with the rectangles, it creates a much more complex polygon, not a simpler one. All the O.P. wants is an algorithm for buffering a polyline. jreid: Doesn't even your old version of PostGIS offer buffering capabilities? Regardless, the algorithm described at gis.stackexchange.com/a/70827/664 answers a more general version of your question. It still relies on the assumption that the GIS can simplify a union of polygons. –  whuber Nov 7 '13 at 20:29
    
It's the polygon they want smoothed. Using a proper buffer algorithm would change the geometry from that which is "fine and dandy". –  Vince Nov 7 '13 at 20:35
1  
The issue is it does not support the proper capping on end points. Its very easy for me to get a buffered line, but the version we have does not allow changing of the end caps nor the join bevel. Ideally I want something that lets me perform this (without PostGIS) ST_Buffer( ST_GeomFromText( 'LINESTRING(50 50,150 150,150 50)' ), 10, 'endcap=flat join=round'); My version does not allow me to change the endcaps or joins (less then 1.5) –  jreid42 Nov 7 '13 at 20:40
    
Have you considered removing the two terminal line segments, buffering what is left, and then unioning the rectangles around the terminal line segments with that buffer? This workaround will perform nicely provided the terminal segments are at least as long as the buffer radius. Another workaround, which is fine when the termini of the polyline do not curve back almost to touch the polyline itself, is to buffer the polyline and then remove two half-disks at either end. (I used to use that method successfully in an old GIS with no control over end caps.) –  whuber Nov 7 '13 at 20:50

1 Answer 1

up vote 1 down vote accepted

I am not sure if this is the most elegant approach, but it worked for me and maybe it will help somebody get started in the right direction. If you have 1.5+ of PostgreSQL then use ST_Buffer with different endcaps - I am sure it performs this in a much more efficient manner).

Because I have the two rectangles and all their coordinates, the biggest problem was insuring that I am joining the exposed open parts and not touching the other side (which side is the kink on?)

graphic

How do I know which points to join?

I've built a "Rectangle" class that has the points TL, TR, BL, BR. They are always oriented along the line.

That way I always know that a TL will link to a BL and a TR will link to a BR.

How do I know whether to join BL to TL or BR to TR?

Based on the difference in angle between the two rectangles I can determine what side has the kink.

In my Rectangle class I am also keeping track of the angle and normalizing it to 0->360 degrees. I don't think this was the best approach but it made it easier for me to understand the problem in my head. I had to do some basic math when the difference in angle was greater then 180 degrees as it causes the kink to flip to the other side.

How do I build the curve?

Once I have identified which points I wish to join (j1 from rectangle 1, j2 from rectangle 2), I need to find a third point somewhere out in the middle (see example above (pa1)). This is used for quadratic bezier curves which I use to build the arc between the points. This third point is used to "pull" the curve away from the rectangles.

The python function I used to build the curves is:

def get_bezier_curve(j1, j2, pa1, segments=8):
    """
    Given 2 points (j1, and j2) calculate the bezier curve
    between them - using point pa1 to pull  the curve away
    from the points.

    Segments determines the number of points calculated for
    the curve.
    """
    # j1 is always first item
    coords = [j1]

    step = 1.0 / segments

    # Quadratic bezier curve
    for t in drange(step, 1.0-step, step):
        x = ((1-t)*(1-t)*j1[0] + 2*(1-t)*t*pa1[0]+t*t*j2[0])
        y = ((1-t)*(1-t)*j1[1] + 2*(1-t)*t*pa1[1]+t*t*j2[1])
        coords.append((x,y))

    # j2 is always last item
    coords.append(j2)
    return coords

To find this, I find the intersection between the inner lines (lines on the side of a rectangle that doesnt have a kink) (i1). Then I create a vector from this towards the original point (p1) that both rectangles share (the original point before buffering). I then increase the length of this to find my point to pull the curve (pa1).

I then build the curve using the above algorithm between j1, and j2 using pa1 pull out the curve.

I then create a polygon out of i1 -> j1 ... bezier curve points ... j2 -> i1.

I then take all the kink polygons, and the rectangle polygons I have and perform a cascaded union (using GEOS from GeoDjango - although you could use Shaply) to collapse these polygons into a single polygon.

Note: To reasons I do not quite fully understand when taking a MultiPolygon of all this data and performing a cascaded union I am not always given a Polygon. Sometimes I am given a GeometryCollection of a useless LineString (no idea where this is coming from) and the combined Polygon. I am guessing this is some side-effect to collapsing certain geometries and stuff not fitting correctly - maybe rounding errors - not sure. If anybody knows why please let me know.

share|improve this answer

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

Not the answer you're looking for? Browse other questions tagged or ask your own question.