Thursday 4 May 2017

sql - Create polygon from lines in PostGIS


I'm trying to divide a polygon in smaller polygons from a line, but I think I can't use the st_split function. What I need is to create small polygons inside a big one, using linestring grid.


I've tried some ways, but I can't get the result. What I've tried:


Divide one polygon from a LineString using st_split()



From a boundary polygon.


enter image description here


And Linestring table:


enter image description here


I'd need the following polygons:


enter image description here


Problem: I can't Split a polygon from several lines, neither a polygon from a multilinestring.


The other method I'm trying is to create a polygons from the lines with st_polygonize() The SQL I was trying is:


SELECT 
g.path[1] as gid,

g.geom::geometry(polygon, 22033) as geom
FROM
(SELECT
(ST_Dump(ST_Polygonize(geom))).*
FROM linestable
) as g;

Extracted from Creating polygons from line segments using PostgreSQL and PostGIS


Problem: I can only get one polygon (the boundary).


Can someone tell me which would be the best way to get the polygons from the linestring, or if am I missing something?



Note: Tables are in the same SRID, and geometries are snapped into a grid. In QGIS I can run the polygonize process from lines to polygon perfectly.


As John's demand, here is the linestring table. https://drive.google.com/file/d/0B603y_m735jfS014S0EyVnpMUEU/view?usp=sharing



Answer



I got this working by using ST_Node first, in conjunction with ST_Collect, to convert the lines into a set of noded linestrings within a MultiLinestring.


As it says in the docs for ST_Node:



Fully node a set of linestrings using the least possible number of nodes while preserving all of the input ones.



What this means, is that all of the linestrings are combined in all possible combinations, so as to make up the equivalent to the exterior ring of a polygon. Whereas, if you attempt to ST_Polygonize a set of LineStrings, none of which on it's own describes a polygon, you simply get the LineStrings back. So, this works:


WITH multi(geom) AS (

SELECT ST_Node(ST_Collect(geom))
FROM leyenda_digitalizar00
)
SELECT ST_AsText( (ST_Dump(ST_Polygonize(geom))).geom )
FROM multi;

If you just run the first part of this, ie, the CTE multi, the output looks like:



MULTILINESTRING((204.5 69.9000000000004,204.5 69.9000000000004),(204.5 68.9,205.4 68.9),(204.5 68.9,204.5 69,204.5 69.1,204.5 69.2,204.5 69.3,204.5 69.4,204.5 69.5,204.5 69.6,204.5 69.7,204.5 69.8,204.5 69.9,204.5 69.9000000000004),(209.5 68.9,209.5 68.8,209.5 68.7,209.5 68.6,209.5 68.5,209.5 68.4,209.5 68.3,209.5 68.2,209.5 68.1,209.5 68,209.5 67.9,209.5 67.8,209.5 67.7,209.5 67.6,209.5 67.5,209.5 67.4,209.5 .......




Now, when you now feed this MultiLinestring to ST_Polygonize it works as expected, eg,



POLYGON((205.4 68.9,204.5 68.9,204.5 69,204.5 69.1,204.5 69.2,204.5 69.3,204.5 69.4,204.5 69.5,204.5 69.6,204.5 69.7,204.5 69.8,204.5 69.9,204.5 69.9000000000004,205.4 69.9,205.4 69.3,205.4 68.9))


POLYGON((204.5 69.9000000000004,204.5 70,204.5 70.1,204.5 70.2,204.5 70.3,204.5 70.4,206.8 70.4,209.5 70.4,209.5 70.3,209.5 70.2,209.5 70.1,209.5 70,209.5 69.9,205.4 69.9,204.5 69.9000000000004))


POLYGON((206.8 70.4,204.5 70.4,204.5 70.5,204.5 70.6,204.5 70.7,204.5 70.8,204.5 70.9,204.5 71,204.5 71.1,204.5 71.2,204.5 71.3,204.5 71.4,206.8 71.4,206.8 70.4))



Obviously, the ST_AsText is just for illustration, and you will have to tweak, if you want the path ID too.


The key takeaway is that ST_Polygonize expects linestrings that already describe the outline of a polygon, which is what ST_Node(ST_Collect(.... does in the above.


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