Search code examples
sqlpostgresqlpostgis

How to create pie wedge (sector) in PostGis by coordinates and azimuth?

I couldn't find correct answer and try to create function in PostGis which returned me polygon with geometry type. This function can be use for visualize cellular network topology on map. Below you can find this function for PostGis which is have next input parameters: lon , lat , azimuth , distance , width. distance - length of pie; width - width of pie.

create or replace function sector_3(lon float, lat float, azimuth float, distance integer, width integer)
returns geometry
language plpgsql
as
$$
declare
   sector geometry;
begin
    sector = ST_MakePolygon(ST_MakeLine(ARRAY[ST_SetSRID(ST_MakePoint(lon,lat),4326),
                ST_Project(ST_SetSRID(ST_MakePoint(lon,lat),4326)::geography, distance, pi()*(azimuth-(width/2))/180.0)::geometry,
                ST_Project(ST_SetSRID(ST_MakePoint(lon,lat),4326)::geography, distance, pi()*(azimuth-(width/2-1*(width/2/5)))/180.0)::geometry,
                ST_Project(ST_SetSRID(ST_MakePoint(lon,lat),4326)::geography, distance, pi()*(azimuth-(width/2-2*(width/2/5)))/180.0)::geometry,
                ST_Project(ST_SetSRID(ST_MakePoint(lon,lat),4326)::geography, distance, pi()*(azimuth-(width/2-3*(width/2/5)))/180.0)::geometry,
                ST_Project(ST_SetSRID(ST_MakePoint(lon,lat),4326)::geography, distance, pi()*(azimuth-(width/2-4*(width/2/5)))/180.0)::geometry,
                ST_Project(ST_SetSRID(ST_MakePoint(lon,lat),4326)::geography, distance, pi()*(azimuth)/180.0)::geometry,
                ST_Project(ST_SetSRID(ST_MakePoint(lon,lat),4326)::geography, distance, pi()*(azimuth+(width/2-4*(width/2/5)))/180.0)::geometry,
                ST_Project(ST_SetSRID(ST_MakePoint(lon,lat),4326)::geography, distance, pi()*(azimuth+(width/2-3*(width/2/5)))/180.0)::geometry,
                ST_Project(ST_SetSRID(ST_MakePoint(lon,lat),4326)::geography, distance, pi()*(azimuth+(width/2-2*(width/2/5)))/180.0)::geometry,
                ST_Project(ST_SetSRID(ST_MakePoint(lon,lat),4326)::geography, distance, pi()*(azimuth+(width/2-1*(width/2/5)))/180.0)::geometry,
                ST_Project(ST_SetSRID(ST_MakePoint(lon,lat),4326)::geography, distance, pi()*(azimuth+(width/2))/180.0)::geometry,                                   
                ST_SetSRID(ST_MakePoint(lon,lat),4326)  
                     ]));
   return sector;
end;
$$;

As a result you will see sector which starts from given point (lon, lat) and with given distance(length) and width.

enter image description here

Solution

  • You could simplify things a bit by first creating a buffer from the parameters distance, lon and lat using ST_Buffer, so that you already have the are you want covered. The second step would be to cut out the piece of the buffer you want, and you can do that using ST_Split. The function ST_Split expects a geometry and a blade to cut the geometry, and this blade can be created as a LineString based on the center point of the buffer (the reference coordinate pair), the azimuth and the distance, using ST_Project. At this point you have two geometries - the area you want and the rest of the buffer - so you can use ST_GeometryN to only return the area that interests you:

    WITH j (geom,azimuth,distance,width) AS (
  VALUES (ST_SetSRID(ST_MakePoint(38.94,45.05),4326),90,120,40)
)
SELECT  
  ST_GeometryN(
    ST_Split(
      ST_Buffer(geom::geography,distance)::geometry,
      ST_MakeLine(ARRAY[
        ST_Project(geom,distance+1,radians(azimuth-(width/2)))::geometry, geom,
        ST_Project(geom,distance+1,radians(azimuth+(width/2)))::geometry]
      )
   ),2)
FROM j;

    enter image description here

    Demo: db<>fiddle

    To better illustrate the proposed solution the next query will depict the buffer, the projected points and the blade created based on them in three different CTEs.

    WITH 
  original_values (geom,azimuth,distance,width) AS (
    VALUES (ST_SetSRID(ST_MakePoint(38.94,45.05),4326),90,120,40)), 
  points (array_points) AS (
    SELECT ARRAY[
      ST_Project(geom,distance+1,radians(azimuth-(width/2)))::geometry, geom,
      ST_Project(geom,distance+1,radians(azimuth+(width/2)))::geometry] 
    FROM original_values),
  line (blade) AS (
    SELECT ST_MakeLine(array_points) 
    FROM points),
  area (buffer) AS (
    SELECT ST_Buffer(geom::geography,distance) 
    FROM original_values   
  ) 
SELECT buffer FROM area
UNION
SELECT unnest(array_points) FROM points 
UNION 
SELECT blade FROM line;

    enter image description here