2

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.

  • @ Cyril, That doesn't help me. I know how to calculate the azimuth if only two points. – FRQ Feb 12 at 17:12
  • It is not about a central point. I Want to calculate azimuth from one point coordinate to each coordinate of the surfaces. Plus I have seen I guess almost all the questions related to ST_Azimuth here. – FRQ Feb 12 at 17:21
  • but you can create one column in your table, where are the coordinates of one point (in accordance with the number of observed points), and another column where are the other points of type point?... – Cyril Feb 12 at 17:23
  • 1
    I think what you want to do is use ST_DumpPoints in a sub-query or CTE and then join this against the point you want to calculate the azimuths for. Your question is a little bit hard to understand, though. – John Powell Feb 13 at 8:19
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

:-)...

  • I have re-edited my question once again. Hope you understand it now. Because the above query is not relevant to my question. – FRQ Feb 13 at 3:57
  • Yeah but that's not I am trying to do. I don't want to duplicate the table either I want to store those duplicate points in a new table. I am trying to figure it out till someone has better expertise and thought or idea to guide me to proceed. – FRQ Feb 13 at 15:07
1

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.

  • It seems that it may be something like that I want, but I am confused about how the query is working. Are you calculating the azimuth between 391139.27, 5821816.69? If so that is not my two different points instead those are the x and y coordinates in 25833. Nevertheless could you please explain your query about how it is working in steps, please. – FRQ Feb 13 at 21:08
  • @user24610 let's take the first query: ST_DumpPoints creates an own record for every point in the surface in the sub-select, quite but not exactly like a table. the outer query takes each row of that sub-select and expands the row type ((dmp).path & (dmp).geom, check the link and further docs) and uses each points geometry in ST_Azimuth, together with your observer point. both are geometry data type (go read about PostGIS standard geometry types before doing anything else...if that is new to you ,) ) that stores both coordinates with it. – ThingumaBob Feb 13 at 21:27
  • @user24610 or have I misunderstood you? anyway, ST_Area returns radians, use degrees(ST_Area(...)) AS azm to convert if you need. – ThingumaBob Feb 13 at 21:36
  • No, I don't need to know the area. I was dumping points to extract x and y as in Spatial reference ID 25833 unit is in meter, and I want all in the same Spatial reference. So at the end is it calculating the azimuth between all surface points and one observer point? That's what I want to do, setting up this query string dynamically in Python. Plus I don't know how with CTE this will work. I am a beginner on my way to learn more and new stuff. – FRQ Feb 13 at 22:03
  • 1
    oh, you meant probably ST_Azimuth and degrees(ST_Azimuth(...)) AS azm instead of ST_Area and degrees(ST_Area(...)) AS azm respectively... Got it – FRQ Feb 14 at 1:34

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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