[postgis-users] Turn a Polygon into Lines

[Randall, Eric]

Lee Hachadoorian wrote:
> I'm looking through the PostGIS reference, and I can't seem to find a 
> way to take a geometry of polygons and turn it into lines.  What I'm 
> looking for is something like the ArcGIS Feature to Line geoprocessor, 
> which will create a line shapefile where each feature is an arc 
> representing the boundary between neighboring polygons with a field 
> indicating the ids of the polygon on either side.
>
> Functions like ST_MakeLine require point geometries, and I don't see 
> anything else that seems to be what I'm looking for.  Any ideas would 
> be welcome.
>
> Thanks,
> Lee Hachadoorian
> PhD Student in Geography
> Program in Earth & Environmental Sciences
> CUNY Graduate Center


A custom function is the fastest and most efficient way to do this, as others have indicated,  
but it's often fun and instructive to try to implement a task in a single statement.
Here's my attempt at the "ArcGIS Feature to Line geoprocessor" mentioned.
It needs to be faster and smarter, but seems to work.  If anyone has suggestions for this or even a different
approach,  please chime in.  I'm not always good at seeing other paths once I get a ways into the woods.
Thanks.

-Eric





/* 
    turn_poly_into_lines.sql

   I used parcel polygons in my experiment.

   Replace "polytable" with your poly table.  Replace gid value with your poly's gid.

 */ 

select p4.poly_id,p4.adjacent_poly_id,p4.pos_order,p4.geom 
from
(

/*  Get edges shared with other polys  */

select t1.gid as poly_id,t2.gid as adjacent_poly_id, 
ST_line_locate_point(ST_linemerge(ST_boundary(t1.geom)),ST_startpoint(ST_linemerge(ST_intersection(ST_boundary(t1.geom),ST_boundary(t2.geom))))) as pos_order,
ST_linemerge(ST_intersection(ST_linemerge(ST_boundary(t1.geom)),ST_linemerge(ST_boundary(t2.geom)))) as geom
from polytable t1, polytable t2
where t1.gid = 309
and ST_touches(t1.geom, t2.geom)

union

/* 
   Create edges not shared with other polys.  Uses -1 for pseudo external universe polygon.  
   Produces extra line necessitating the self-join and the ST_covers filter at the end.
   Need help, more work on this.
*/

select p1.poly_id, p1.adjacent_poly_id, p1.position as pos_order,ST_line_substring(p3.geom,p1.position,p2.position)  as geom
from
(
select t1.gid as poly_id,-1::integer as adjacent_poly_id,
ST_line_locate_point(ST_linemerge(ST_boundary(t1.geom)),ST_endpoint(ST_linemerge(ST_intersection(ST_boundary(t1.geom),ST_boundary(t2.geom)))))
as position, ST_endpoint(ST_linemerge(ST_intersection(ST_boundary(t1.geom),ST_boundary(t2.geom)))) as geom
from polytable t1, polytable t2
where t1.gid = 309
and ST_touches(t1.geom, t2.geom)
and ST_endpoint(ST_linemerge(ST_intersection(ST_boundary(t1.geom),ST_boundary(t2.geom)))) not in
(
select ST_startpoint(ST_linemerge(ST_intersection(ST_boundary(t1.geom),ST_boundary(t2.geom)))) as geom
from polytable t1, polytable t2
where t1.gid = 309
and ST_touches(t1.geom, t2.geom)
)
)
as p1,
(
select t1.gid as poly_id,-1::integer as adjacent_poly_id, 
ST_line_locate_point(ST_linemerge(ST_boundary(t1.geom)),ST_startpoint(ST_linemerge(ST_intersection(ST_boundary(t1.geom),ST_boundary(t2.geom))))) 
as position, ST_startpoint(ST_linemerge(ST_intersection(ST_boundary(t1.geom),ST_boundary(t2.geom)))) as geom
from polytable t1, polytable t2
where t1.gid = 309
and ST_touches(t1.geom, t2.geom)
and ST_startpoint(ST_linemerge(ST_intersection(ST_boundary(t1.geom),ST_boundary(t2.geom)))) not in
(
select ST_endpoint(ST_linemerge(ST_intersection(ST_boundary(t1.geom),ST_boundary(t2.geom)))) as geom
from polytable t1, polytable t2
where t1.gid = 309
and ST_touches(t1.geom, t2.geom)
)
)
as p2,
(
select ST_linemerge(ST_boundary(geom)) as geom
from polytable 
where gid = 309
)
as p3
where p1.position < p2.position
)
as p4,

/* Repeat for needed self-join */

(

/* Get edges shared with other polys */

select t1.gid as poly_id,t2.gid as adjacent_poly_id, 
ST_line_locate_point(ST_linemerge(ST_boundary(t1.geom)),ST_startpoint(ST_linemerge(ST_intersection(ST_boundary(t1.geom),ST_boundary(t2.geom))))) as pos_order,
ST_linemerge(ST_intersection(ST_linemerge(ST_boundary(t1.geom)),ST_linemerge(ST_boundary(t2.geom)))) as geom
from polytable t1, polytable t2
where t1.gid = 309
and ST_touches(t1.geom, t2.geom)

union

/* 
   Create edges not shared with other polys.  Uses -1 for pseudo external universe polygon.  
   Produces extra line necessitating the self-join and the ST_covers filter at the end.
   Need help, more work on this.
*/

select p1.poly_id, p1.adjacent_poly_id, p1.position as pos_order, ST_line_substring(p3.geom,p1.position,p2.position)  as geom
from
(
select t1.gid as poly_id,-1::integer as adjacent_poly_id,
ST_line_locate_point(ST_linemerge(ST_boundary(t1.geom)),ST_endpoint(ST_linemerge(ST_intersection(ST_boundary(t1.geom),ST_boundary(t2.geom)))))
as position, ST_endpoint(ST_linemerge(ST_intersection(ST_boundary(t1.geom),ST_boundary(t2.geom)))) as geom
from polytable t1, polytable t2
where t1.gid = 309
and ST_touches(t1.geom, t2.geom)
and ST_endpoint(ST_linemerge(ST_intersection(ST_boundary(t1.geom),ST_boundary(t2.geom)))) not in
(
select ST_startpoint(ST_linemerge(ST_intersection(ST_boundary(t1.geom),ST_boundary(t2.geom)))) as geom
from polytable t1, polytable t2
where t1.gid = 309
and ST_touches(t1.geom, t2.geom)
)
)
as p1,
(
select t1.gid as poly_id,-1::integer as adjacent_poly_id, 
ST_line_locate_point(ST_linemerge(ST_boundary(t1.geom)),ST_startpoint(ST_linemerge(ST_intersection(ST_boundary(t1.geom),ST_boundary(t2.geom))))) 
as position, ST_startpoint(ST_linemerge(ST_intersection(ST_boundary(t1.geom),ST_boundary(t2.geom)))) as geom
from polytable t1, polytable t2
where t1.gid = 309
and ST_touches(t1.geom, t2.geom)
and ST_startpoint(ST_linemerge(ST_intersection(ST_boundary(t1.geom),ST_boundary(t2.geom)))) not in
(
select ST_endpoint(ST_linemerge(ST_intersection(ST_boundary(t1.geom),ST_boundary(t2.geom)))) as geom
from polytable t1, polytable t2
where t1.gid = 309
and ST_touches(t1.geom, t2.geom)
)
)
as p2,
(
select ST_linemerge(ST_boundary(geom)) as geom
from polytable 
where gid = 309
)
as p3
where p1.position < p2.position
)
as p5

where ST_covers(p4.geom,p5.geom)
group by p4.poly_id,p4.adjacent_poly_id,p4.pos_order,p4.geom
having count(ST_covers(p4.geom,p5.geom)) = 1
order by p4.pos_order




---------
 
Eric Randall
County of Erie
140 W 6th St
Room 119
Erie, PA 16501

ph. 814-451-6063
fx. 814-451-7000

-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.osgeo.org/pipermail/postgis-users/attachments/20080331/5c5af262/attachment.html>