# Python library or algorithm to generate arc geometry from three coordinate pairs?

I'm trying to write code to generate fixtures for a django/geodjango project. I need to generate geometry for arcs (line string) given three lat/long pairs; start of arc, end of arc and centre of arc. I read through the geodjango documentation but this functionality doesn't seem to be available. Does anyone know if there is a Python library that provides this functionality? Or does anyone have an algorithm I can port to python? Thanks

• what kind of arc: a LineString, a curve ? – gene Jul 11 '14 at 17:07
• Oh sorry, a line string. – Mike Stoddart Jul 11 '14 at 18:48
• If you're underlying database is postgis, it supports curves, IE ST_GeomFromText('CIRCULARSTRING(220268 150415,220227 150505,220227 150406) from postgis.net/docs/ST_CurveToLine.html – DPierce Jul 24 '14 at 3:14
• Thanks. I was originally using a plpgsql function that generates an arc from the three coordinate pairs. Unfortunately, we are re-designing our installer, which means the data is now written to json files instead of being loaded straight into PostGIS/PostgreSQL. So I can't use the function anymore and I have to generate the arg geometry using Python. – Mike Stoddart Jul 24 '14 at 10:50

## 2 Answers

Shapely PyQGIS, and GeoDjango use the same API based on the GEOS library:

### With Shapely:

with a list of points:

``````from shapely.geometry import Point, LineString, mapping
pt1 = Point(0,0)
pt2 = Point(20,20)
pt3 = Point (50,50)
line = LineString([pt1,pt2,pt3])
#GeoJSON format
print mapping(line)
{'type': 'LineString', 'coordinates': ((0.0, 0.0), (20.0, 20.0), (50.0, 50.0))}
``````

or with a list of coordinates:

``````line = LineString([(0, 0), (20,20),(50,50)])
#GeoJSON format
print mapping(line)
{'type': 'LineString', 'coordinates': ((0.0, 0.0), (20.0, 20.0), (50.0, 50.0))}
``````

### with PyQGIS:

with a list of points:

``````line = QgsGeometry.fromPolyline([QgsPoint(0,0),QgsPoint(20,20),QgsPoint(50,50)])
#GeoJSON format
print line.exportToGeoJSON()
{ "type": "LineString", "coordinates": [ [0, 0], [20, 20], [50, 50] ] }
``````

### with GeoDjango (look at Django: GEOS API)

with a list of coordinates:

``````from django.contrib.gis.geos import LineString
line = LineString((0, 0), (20, 20), (50, 50))
# GeoJSON format
print line.json
{ "type": "LineString", "coordinates": [ [ 0.0, 0.0 ], [ 20.0, 20.0 ], [ 50.0, 50.0 ] ] }
``````

The resulting line is made up of two segments: (0,0) to (20,20) and (20,20) to (50,50) • Thanks but this doesn't give me an arc, it only gives me a line between three points. – Mike Stoddart Jul 12 '14 at 0:52
• ""Oh sorry, a line string" -> what is an arc ? a segment of a LineString ? – gene Jul 12 '14 at 7:47
• I found this image after a quick search that hopefully shows what I'm after. I may be using the wrong terminology so I apologise if I'm confusing people. Refer to item "2" in this image: docs.autodesk.com/CIV3D/2013/ENU/images/…. The arc is comprised of multiple line segments. – Mike Stoddart Jul 21 '14 at 11:34
• In item 2, its is a polyline with three segments (look above) – gene Jul 21 '14 at 14:04
• two segments, sorry – gene Jul 24 '14 at 16:01

Building on @gene 's answer above if you're looking to approximate a curve as a sequence of points using spline interpolation. In Python you can do this through the `scipy.interpolate` library. Particularly 1d interpolation.

For example

``````import numpy as np
import scipy.interpolate

coords = np.array([[0, 0], [25, 10], [50, 50]])

#The curve fits as a quadratic equation on three points
f = scipy.interpolate.interp1d(coords[:, 0], coords[:, 1], kind='quadratic')

#New points will be evenly distributed along x
new_x = np.linspace(np.min(coords[:, 0]), np.max(coords[:, 0]), 10)
new_y = f(new_x)

new_coords = np.vstack([new_x, new_y]).T
`````` With a bit more work (play with `shapely.interpolate`) you can get a curve in segments of equal like so:

``````from shapely.geometry import Point, LineString, mapping

fine_x = np.linspace(np.min(coords[:, 0]), np.max(coords[:, 0]), 1000)
fine_y = f(fine_x)

fine_coords = zip(fine_x, fine_y)
fine_line = LineString(fine_coords)

even_line = LineString([
np.array(fine_line.interpolate(i))
for i in np.arange(0, fine_line.length, 5) #points 5 units apart
])
`````` • you can simplify: `fine_coords = zip(fine_x, fine_y)` – gene Jul 24 '14 at 8:22
• Thanks but I'm looking to generate segments using three coordinates; start of arc, end of arc and centre of arc. – Mike Stoddart Jul 24 '14 at 11:03
• @Stodge In which case take a look at this answer on Stack Overflow - it should have exactly what you're after – om_henners Jul 24 '14 at 15:43
• The result is not a LineString with three points. Otherwise you can also use all the interpolation algorithms. – gene Jul 24 '14 at 15:57
• @gene Yeah. I'm now thinking Stodge is after a circular arc as generalised to a line string. After looking at it a couple of times, in the picture supplied in the comment to your solution I think the centre point of the arc is the green circle to the left of the "2" – om_henners Jul 24 '14 at 16:05