Tuesday, 7 January 2020

Drawing wind rose with QGIS from PostGIS?


I am working on wind rose data. I have some points on lines (streets) in my Postgres/PostGIS DB. The sample data scenario is:


A point on line (street)


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.


expected output


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.


Wind rose using PostGIS



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

arcpy - Changing output name when exporting data driven pages to JPG?

Is there a way to save the output JPG, changing the output file name to the page name, instead of page number? I mean changing the script fo...