Intersection between an angle and a polygon (postgreSQL + postGIS or Python2.7 + shapely + psycopg2)

Question

I have a point geom feature (in wkb) in a PostgreSQL database table (with postGIS extension). There is also an azimuth value (degrees), as float, in an other field of the same table.

How to perform, as simple as it would be to cut a part of a cake, an intersection between a x[m] buffer around the point and a 20° angle around the azimuth (±10° on both side, the azimuth should be the bisecting line of the angle)?

From now I've calculated the buffer as an other geom (as a polygon also in wkb) which is stored in an other field of the table.

I wonder what would be the faster and more CPU effective calculation to intersect this circle with the angle?

The only postGIS function I have found yet is ST_azimuth: http://www.postgis.org/docs/ST_Azimuth.html but it may be rather "complicated" to implement as it would need to calculate two "virtual secondary points" that are "far enough" from the buffer and apart respectively from -10 and +10° from the existing azimuth and then make the difference between these two new azimuth angles... And after, I really can't figure out how to intersect an angle, which is some kind of "not clearly defined" geometry, with the circular buffer.

An other point; at the end I would need to implement this in a Python script (I'm basically using shapely and psycopg2 modules yet). But for now it could be done in an SQL query directly on the postgreSQL DB.

--

http://toblerity.org/shapely/manual.html

http://initd.org/psycopg/

Answer 1

Answer:

In fact, one don't really have to calculate (unless high precision is needed) the resulting intersection between the buffer around the point and an angle of ± 10° (or whatever angle) on both sides from the azimuth direction.

On can simply build the triangle. Let's have a look...

Code:

import math
from shapely import geometry
OFFSET = 1000 # Set a big enough offset from the initial point to suit your needs
azim   = 148                 # degrees
point0 = geometry.Point(0,0) # center initial point geometry
E1     = OFFSET*math.sin(math.pi*(azim-10)/180.) # east coord. of 1st secondary point
E2     = OFFSET*math.sin(math.pi*(azim+10)/180.) # east coord. of 2nd secondary point
N1     = OFFSET*math.cos(math.pi*(azim-10)/180.) # north coord. of 1st secondary point
N2     = OFFSET*math.cos(math.pi*(azim+10)/180.) # north coord. of 2nd secondary point
point1 = geometry.point.Point((point0.x+E1, point0.y+N1))
point2 = geometry.point.Point((point0.x+E2, point0.y+N2))
coords = ((point0.x, point0.y),(P1.x, P1.y),(P2.x, P2.y),(point0.x, point0.y))
triangle   = geometry.Polygon(coords)
# If you really need the intersection with a circular buffer around the
# initial point:
BUFFSIZE   = 600 # must be lower than OFFSET, otherwise useless.
buff       = point0.buffer(BUFFSIZE)
intersec   = triangle.intersection(buff)

That's it.
Result look nice, here some intersections between circular buffers around some points and some given azimuths, with a ± 20° "angular buffer" around these azimuths: