Create area that cell tower covers

Question

I have a dataset with cell tower information, as you can see below. The lat and lon fields are the location of the tower.

The objective is to get the area (geometry) of the sector the cell tower covers, from start_angle to end_angle. As you can see in the next image, using the first row of the dataset as example, I can obtain the lines from start_angle 275 and end_angle 35, but I want the rest of the buffer to disappear.

Query used:

This first query is used to create and rotate the lines from start_angle and end_angle and also a line with 0 degrees.

WITH vertices AS
(SELECT 
      id,
 start_angle,
 end_angle,
      (ST_DumpPoints(geom)).path[1] AS v_id, 
      (ST_DumpPoints(geom)).geom AS vertex
FROM celulas
), teste as
(SELECT
    id,
    v_id,
    ST_SetSRID(ST_Translate(ST_Rotate(ST_MakeLine(ST_MakePoint( 1.0,0.0), 
                                                  ST_MakePoint(-1.0,0.0)),
                                      radians(start_angle * -1)), ST_X(vertex), ST_Y(Vertex)),
                                      ST_SRID(vertex)) AS startline,
    ST_SetSRID(ST_Translate(ST_Rotate(ST_MakeLine(ST_MakePoint( 1.0,0.0), 
                                                  ST_MakePoint(-1.0,0.0)),
                                      radians(end_angle * -1)), ST_X(vertex), ST_Y(Vertex)),
                                      ST_SRID(vertex)) AS endline,
    ST_SetSRID(ST_Translate(ST_Rotate(ST_MakeLine(ST_MakePoint( 1.0,0.0), 
                                                  ST_MakePoint(-1.0,0.0)),
                                      radians(0)), ST_X(vertex), ST_Y(Vertex)),
                                      ST_SRID(vertex)) AS midline
    
FROM vertices
)

I also used the next query to union all geometries: the radius buffer and the lines

select St_intersection(st_split(buffer, midline), st_split(buffer, angulo))
from angulo

Hey there. What do you exactly mean by "I want the rest of the buffer to disappear"? — Jim Jones, Sep 23 '20 at 10:31
Hello. As you can see in the image I have the entire radius of the cell tower. But I only want to show the sector from the start_angle to the end_angle, i.e., from 275 degrees to 35. The rest of the circle doesnt matter. — Gonçalo Silva, Sep 23 '20 at 10:50
Okay. Can you also add to your question the query you used to create this geometry/buffer? — Jim Jones, Sep 23 '20 at 11:23

geozelot · Answer 1 · 2020-09-25T13:17:18.023

Let's create us a handy inline function on the composite type of your relation specifically.

Also, and first off, let's create two equally handy support functions for all sorts of use cases in this context; getting the clockwise (CW) (and, since we're at it, also the counterclockwise (CW)) angular difference in degree between two azimuths:

CREATE OR REPLACE FUNCTION CWAngle(
  IN    sdeg FLOAT,
  IN    edeg FLOAT,
  OUT   ddeg  FLOAT
) LANGUAGE SQL AS
  $$
    SELECT CASE (sdeg <= edeg)
             WHEN TRUE THEN
               edeg - sdeg
             ELSE
               360.0 - sdeg + edeg
           END;
  $$
;

CREATE OR REPLACE FUNCTION CCWAngle(
  IN    sdeg FLOAT,
  IN    edeg FLOAT,
  OUT   ddeg  FLOAT
) LANGUAGE SQL AS
  $$
    SELECT 360.0 - CWAngle(sdeg, edeg);
  $$
;

There may or may not be more efficient ways, or more idiomatic code; this works well enough.

Now, the concept of the following function is to create a 'wedge', your sector, of a circle by projecting the center point along the circle segment between start_angle and end_angle at the given "radius"; we will project <degrees_between_angles>/FLOOR(<degrees_between_angles>) points, so every 1 <= step < 2, with step being degrees (which is a lot; ST_Buffer defaults to a circle having 8 vertices per quarter circle):

CREATE OR REPLACE FUNCTION sector(
    rec spatial.celulas
) RETURNS GEOMETRY(POLYGON, 4236) AS
  $$
    DECLARE
      delta FLOAT := CWAngle(rec.start_angle, rec.end_angle);
      step  FLOAT := delta / FLOOR(delta);
      
      wedge GEOMETRY(POINT)[];
      
    BEGIN
      wedge := wedge || ST_SetSRID(ST_MakePoint(rec.lon, rec.lat), 4326);
  
      FOR n IN 0..FLOOR(delta)
        LOOP
          wedge := wedge || ST_Project(wedge[1]::GEOGRAPHY, rec."radius", RADIANS(MOD(rec.start_angle+(n*step)::NUMERIC, 360.0::NUMERIC)))::GEOMETRY;
        END LOOP
      ;
  
      wedge := wedge || wedge[1];
      
      RETURN ST_MakePolygon(ST_SetSRID(ST_MakeLine(wedge), 4326));
    END;
  $$
LANGUAGE 'plpgsql';

This assumes

that the CWAngle function is present; change to CCWAngle in this function if needed
that start_angle -> end_angle is in clockwise direction
that "radius" is in meter
the exact relation and column names of your example, i.e celulas.lat, celulas.lon, celulas."radius", celulas.start_angle, celulas.end_angle; you need to alter all occurences of these identifiers in the function accordingly if they are actually different

This function is thus specific to the relation in your example, and you can call it like a relation qualified column:

SELECT [*,] celulas.sector AS geom
FROM   celulas
;

Alternatively, here is a more general function based on the same concept, which also allows to set the max vertices per quarter circle, and an optional inner radius.

Create area that cell tower covers

1 Answers1