Wednesday, 2 October 2019

Filling a polygon with lines using Postgis


I need to calculate a specific distance in an area. The distance is calculated from parallel lines in a polygon, And I need to find out how many parallel lines (with x meters distance to each other) can a fill in a polygon. Is there some ready function for this operation in postgis or mapscript? Or do I need to implement it myself?



Here is sample image for what I want to do:


enter image description here


I want to calculate how many orange lines can I put in polygon, and total length of these lines.



Answer



This can be done quite easy in SQL


All the below examples can be tested directly on http://postgisonline.org/map.php. Just paste the query and press Map1


SELECT GENERATE_SERIES(FLOOR(ST_YMin(the_polygon))::int , CEILING(ST_YMax(the_polygon))::int,200) y_value, ST_XMin(the_polygon) x_min, ST_XMax(the_polygon) x_max from 
(SELECT the_geom AS the_polygon FROM lakes) l

Next step is to build lines from those coordinates and and use ST_Intersection to only get the parts of the lines intersecting the polygon:



SELECT ST_Intersection(the_geom, the_polygon)  AS the_geom FROM
(SELECT the_polygon, ST_Setsrid(ST_MakeLine(ST_MakePoint(x_min, y_value),ST_MakePoint(x_max, y_value) ), ST_Srid(the_polygon)) AS the_geom FROM
(SELECT the_polygon, GENERATE_SERIES(FLOOR(ST_YMin(the_polygon))::int , CEILING(ST_YMax(the_polygon))::int,200) y_value, ST_XMin(the_polygon) x_min, ST_XMax(the_polygon) x_max from
(SELECT the_geom AS the_polygon FROM lakes) l
)c
) lines

Then, at last we can just sum all the lengths. So what you have got is a query looking like this doing the whole operation:


SELECT SUM(ST_Length(the_geom)) FROM 
(SELECT ST_Intersection(the_geom, the_polygon) AS the_geom FROM

(SELECT the_polygon, ST_Setsrid(ST_MakeLine(ST_MakePoint(x_min, y_value),ST_MakePoint(x_max, y_value) ), ST_Srid(the_polygon)) AS the_geom FROM
(SELECT the_polygon, GENERATE_SERIES(FLOOR(ST_YMin(the_polygon))::int , CEILING(ST_YMax(the_polygon))::int,200) y_value, ST_XMin(the_polygon) x_min, ST_XMax(the_polygon) x_max from
(SELECT the_geom AS the_polygon FROM lakes) l
)c
) lines
) intersection_lines

edit


This was the length part of your problem mentioned in the end of the question. To count the number of lines is no problem of course. Then just use count(*) on the top row.


HTH



Nicklas


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...