I am working on wind rose data. I have some points on lines (streets) in my Postgres/PostGIS DB. The sample data scenario is:
First, I computed street orientation of this line using PostGIS commands ST_StartPoint()
, ST_EndPoint()
and ST_Azimuth()
like this:
Select
-- line_dir = line direction (orientation) in degrees
(degrees(ST_Azimuth(line_end_pt, line_start_pt)))::int As line_dir,
-- assign the line direction as the lower bound of sector 1 (of 30 degree)
(degrees(ST_Azimuth(line_end_pt, line_start_pt)))::int As lower_bound_sector1,
-- add 30 degree to lower bound to get the upper bound as each sector is of 30 degree max.
(degrees(ST_Azimuth(line_end_pt, line_start_pt)))::int + 30 As upper_bound
from
(
Select
pt_on_line.geom as pt,
line.geom As line,
ST_StartPoint(ST_LineMerge(line.geom)) As line_start_pt,
ST_EndPoint(ST_LineMerge(line.geom)) As line_end_pt
from pt_on_line
Left join line on ST_DWithin(pt_on_line.geom, line.geom, 0.1)
)
As foo;
Now starting with this street orientation, I need to generate a wind rose of 12 sectors (each of 30 degree having its lower bound and upper bounds) as geometry around this point on line in a radius of 50 meter. In the above code, I have assigned street direction as lower bound of sector_1 and added 30 degree to the lower bound. This approach might be unrealistic though. The expected output could be something similar to this.
How can I draw these 12 wind sectors (geometry) around this point starting from line's orientation?
Using the answer of @MichalZimmermann, here is the screen shot of wind rose generation using PostGIS and visualized in QGIS.
Answer
CREATE OR REPLACE FUNCTION ST_WindRose(
line geometry,
directions int,
radius numeric
)
RETURNS TABLE (
id integer,
geom geometry(LINESTRING)
)
AS $ST_WindRose$
BEGIN
IF directions % 2 <> 0 THEN
RAISE EXCEPTION 'Odd number of directions found, please provide even number of directions instead.';
END IF;
IF radius > ST_Length(line) THEN
RAISE EXCEPTION 'Inner circle radius is bigger than the wind rose diameter, please make it smaller.';
END IF;
RETURN QUERY
WITH rose AS (
SELECT
ST_Rotate(_line, radians(360) / directions * dirs.id, ST_Centroid(_line)) _line
FROM (
SELECT line _line
) a
CROSS JOIN (
SELECT generate_series(1, directions / 2) id
) dirs
)
SELECT
row_number() OVER ()::integer id,
_line geom
FROM (
SELECT _line FROM rose
UNION ALL
SELECT ST_ExteriorRing(ST_Buffer(ST_Centroid(line), radius, 30)) -- inner circle
UNION ALL
SELECT ST_ExteriorRing(ST_Buffer(ST_Centroid(line), ST_Length(line)/2, 30)) -- outer circle
) a;
END
$ST_WindRose$
LANGUAGE PLPGSQL;
Use it with
SELECT * FROM ST_WindRose(ST_MakeLine(ST_MakePoint(0,0), ST_MakePoint(0,1)), 12, 0.01);
Note it's half past 9pm here, it might not be the best solution ever, but it works.
No comments:
Post a Comment