Azimuth calculation from One to Many Points in PostGIS

Question

I have vector data of 3D City Model in PostgreSQL database. So from my data geometry column, which is a surface geometry. (These are geometries of all the surfaces of buildings). I had extracted the points through ST_DumpPoints query as follows:

SELECT (ST_x(geom)) AS X, ST_y(geom) AS Y FROM (
     SELECT (St_DumpPoints(geometry)).geom from surface_geometry WHERE cityobject_id=95
   ) AS foo

Now I would like to calculate the azimuth from the origin (the below-created point) to all the points extracted with the above query from all the surface polygons. Considering this point created below as position, viewpoint or GPS position of an observer, I am building it with MakePoint as below. (So as observer position, I would be changing this point everytime or is there any way to automatise this?)

SELECT (ST_x(StandPoint)) AS X, ST_y(StandPoint) AS Y FROM (
     SELECT ST_SetSRID(ST_MakePoint(391139.27, 5821816.69), 25833) AS StandPoint) AS 
 foo

I have seen in PostGIS that there is a way to calculate azimuth between two points. The function ST_Azimuth calculate between two points as in PostGIS manual. But how to calculate azimuth from one to many? Any tips on how to proceed?

Later I want to store that calculate the azimuth in a variable in python using PsyCopg 2 (This part is not difficult) for my further calculations.

Answer 1

If I understand the question correctly, the solution is the following

1) Run the script:

create table azimuth_1_m as
select a.gid as gid, b.gid as id, 
degrees(st_azimuth(a.geom,b.geom)) as azimuth
from
your_source_table as a, your_source_table as b
where a.gid=1

You can choose any point...a.gid=2,3,4, etc.

2) To get rid of empty lines, execute the following script

create table azimuth_1_many as
select * from azimuth_1_m 
where azimuth::character varying IS NOT NULL  
and azimuth::character varying <> ''
and gid IS NOT NULL
and id IS NOT NULL

:-)...

Answer 2

A join is what you are looking for; ST_DumpPoints returns a record set that, when expanded, can be joined as a table to your StandPoint.

---

In this simple case however, you can even just expand the set and run ST_Azimuth with your explicit point generation within; the basic concept would be:

SELECT  (dmp).path[1],
        ST_Azimuth((dmp).geom, ST_SetSRID(ST_MakePoint(391139.27, 5821816.69), 25833)) AS azm
FROM    (
    SELECT  ST_DumpPoints(geometry) AS dmp
    FROM    surface_geometry
    WHERE   cityobject_id=95
) q

Or, which I find more elegant, with the use of LATERAL to expand the record set:

SELECT  dmp.path[1] AS id,
        ST_Azimuth(dmp.geom, ST_SetSRID(ST_MakePoint(391139.27, 5821816.69), 25833)) AS azm
FROM    surface_geometry,
        LATERAL ST_DumpPoints(geom) AS dmp
WHERE   cityobject_id=95

In both cases, I use the path[1] element to generate an id corresponding to the order of the vertices of the surface.

This is a simple use-case for a single point created at run-time and a filter on your surface table; you could set up this query string dynamically in Python, i.e replace the point coordinates and the filter id.

---

As you hinted at automation: it is also perfectly possible to use this as the core for a two-table join (possibly with a CTE) between all objects and a collection of observer positions; in that case you'd need to provide some sort of id filter between the two tables to limit one observer to one city object. As your question stands it's hard to predict your actual scenario.